Whole Genome Sequencing Results Associated with Minimum Inhibitory Concentrations of 14 Anti-Tuberculosis Drugs among Rifampicin-Resistant Isolates of Mycobacterium Tuberculosis from Iran

Accurate and timely detection of drug resistance can minimize the risk of further resistance development and lead to effective treatment. The aim of this study was to determine the resistance to first/second-line anti-tuberculosis drugs in rifampicin/multidrug-resistant Mycobacterium tuberculosis (RR/MDR-MTB) isolates. Molecular epidemiology of strains was determined using whole genome sequencing (WGS)-based genotyping. A total of 35 RR/MDR-MTB isolates were subjected to drug susceptibility testing against first/second-line drugs using 7H9 Middlebrook in broth microdilution method. Illumina technology was used for paired-end WGS applying a Maxwell 16 Cell DNA Purification kit and the NextSeq platform. Data analysis and single nucleotide polymorphism calling were performed using MTBseq pipeline. The genome-based resistance to each drug among the resistant phenotypes was as follows: rifampicin (97.1%), isoniazid (96.6%), ethambutol (100%), levofloxacin (83.3%), moxifloxacin (83.3%), amikacin (100%), kanamycin (100%), capreomycin (100%), prothionamide (100%), D-cycloserine (11.1%), clofazimine (20%), bedaquiline (0.0%), and delamanid (44.4%). There was no linezolid-resistant phenotype, and a bedaquiline-resistant strain was wild type for related genes. The Beijing, Euro-American, and Delhi-CAS were the most populated lineage/sublineages. Drug resistance-associated mutations were mostly linked to minimum inhibitory concentration results. However, the role of well-known drug-resistant genes for D-cycloserine, clofazimine, bedaquiline, and delamanid was found to be more controversial.


Introduction
Multidrug-resistant tuberculosis (MDR-TB), with a wide range of prevalence around the world is considered as a global public health concern [1]. In 2017, approximately 558,000 (483,000-639,000) rifampicin-resistant tuberculosis (RR-TB) cases were estimated, of which 82% were MDR. Overall, 3.5% of new cases and 18% of previously treated cases were RR/MDR-TB and reported to be 1.3% and 8.3%, respectively, in Iran [2]. However, drug resistance in the neighboring countries of Iran is remarkably higher, as in Afghanistan, Armenia, and Azerbaijan, MDR incidence was 4.1%, 11%, and 13% among new cases, and it was 16%, 47%, and 39% among previously treated patients, respectively [2].
Treatment of patients with MDR-TB is more time-consuming and requires expensive drugs with more toxicity, compared to drug-susceptible TB. Additionally, inappropriate and inadequate treatment increases the risk of resistance to other drugs, possible transmission, and prolonged treatment [3]. In other words, accurate drug susceptibility testing (DST) is required to achieve effective treatment and minimize the risk of further resistance development [4].
Spontaneous mutations in the genes encoding drug targets or drug-converting enzymes are the main cause of acquired resistance in Mycobacterium tuberculosis complex (MTBC) [3]. Other mechanisms can also be involved, such as drug efflux pumps [5]. Several phenotypic DST (pDST) methods have become available and are done on bacterial cultures after a positive primary culture of the specimen, and thus detecting resistance is technically complex and time-consuming [6]. Furthermore, these techniques have their own problems with test standardization and stability of various antibiotics in different media [7]. Detection of mutations in drug-resistant conferring genes has been proposed as a fast tool for identifying resistance against anti-tuberculosis drugs. However, the sensitivity of molecular methods to detect the resistance of different drugs is debatable [8]. Molecular techniques like Xpert MTB/RIF (Cepheid, Sunnyvale, CA, USA) and line probe assays did not cover all the drugs (e.g., linezolid, bedaquiline, delamanid, etc.) and have high sensitivities for drugs which only a few mutations within a gene are responsible for resistance (e.g., rifampicin) [9]. Additionally, these methods show low sensitivity for heteroresistant strains with low frequencies of mutant variants [10].
Unfortunately, quantitative phenotypic DST (MICs) data are often limited and genetic background of resistance for some antibiotics were not clearly understood [11]. Accurate identification of chromosomal mutations is accessible using whole genome sequencing (WGS) and linking with the minimum inhibitory concentrations are necessarily needed to better infer the susceptible strains from the resistant ones [12].

Study Design and Mycobacterial Strains
Rifampicin-resistant M. tuberculosis strains isolated from patients referred to Tehran Regional Reference Laboratory for Tuberculosis (TRL-Tehran) between 2014 and 2018 were included in the present study. M. tuberculosis was identified by biochemical tests after positive primary culture on Lowenstein Jensen (LJ) medium [13]. Resistance to the first-line drugs was determined by the proportion method [14]. In brief, a 1.0 McFarland suspension of bacterial colonies were prepared and dilutions of 10 −2 and 10 −4 were inoculated onto LJ medium with 40 µg/mL RIF (Sigma Aldrich, St. Louis, MO, USA) and without drug as test control. After 21 days of incubation at 37 • C, the results were considered sensitive when the proportion of up to 1% of colony growth was seen in the tubes containing the drug. TRL-Tehran undergoes annual external quality control by the World Health Organization (WHO) TB Supranational Reference Laboratory Network in Stockholm, Sweden. Selected isolates were transported to the TB Supranational Reference Laboratory (SRL-Milan) to undergo WGS.

Broth Microdilution
All drugs were purchased from Sigma Aldrich except for BDQ (BOC Sciences, New York, NY, USA) and DLM (Toronto Research Chemicals, North York, ON, Canada). The MICs of first-and second-line drugs were determined by broth micro dilution method using 96-well microplates with OADC (Oleic Albumin Dextrose Catalase) supplemented 7H9 Middlebrook broth [4]. In summary, the final bacterial concentrations were a 10 −2 dilution of prepared 1.0 MacFarland bacterial suspensions. Medium without any antibiotics inoculated with bacterial suspensions used as growth control wells and negative controls were not inoculated. The plates were sealed in plastic bags and incubated 14 days at 37 • C. Growth and negative-control wells were checked to ensure adequate growth and no contamination, respectively. The first well without any visible growth was as the minimum drug concentration could inhibit M. tuberculosis growth. The 2-fold serially diluted drug concentration ranges were as following; 0.12-16 µg/mL for RIF, 0.02-3.2 for INH, 0.25-32 µg/mL for EMB, 0.12-16 µg/mL for LFX, 0.03-4 µg/mL for MFX, 0.06-32 µg/mL for AMK, 0.15-20 µg/mL for KAN, 0.3-40 µg/mL for CAP, 0.3-40 µg/mL for PTO, 0.5-64 µg/mL for DCS, 0.12-16 µg/mL for CFZ, 0.03-8 µg/mL for LZD, 0.015-4 µg/mL for BDQ, and 0.015-3.84 µg/mL for DLM. Since there were no well-defined critical concentrations (CCs) of all anti-TB drugs for drug susceptibility testing on 7H9 broth, drug concentration ranges and the assumed CCs for tested drugs were based on the previous literature [15][16][17]. This experiment was performed in duplicate and M. tuberculosis strain H37Rv ATCC 27294 was used for quality control.
SNP positions, with a reliable base call (10-fold coverage and 75% allele frequency) in at least 95% of the isolates, covered in all samples were concatenated to a sequence alignment, excluding SNPs within a window of 12 bp from each other in the same isolate. From the aligned sequences of concatenated SNPs of all 35 isolates, maximum likelihood trees were calculated using RAxML version 8 [25], with a general time reversible (GTR) substitution model, 1000 resamples, and Gamma20 likelihood optimization to account for rate heterogeneity among sites. The consensus tree was rooted with the "midpoint root" option in FigTree and annotated using the Evolview software. All raw sequencing reads have been deposited at Bioproject at NCBI (project accession nr: PRJNA566379).

Statistical Analysis
Statistical analysis was performed using MedCalc for Windows, version 18.10 (MedCalc Software, Ostend, Belgium), to define the sensitivity, specificity, and accuracy of WGS in predicting each tested drug resistance, compared to the pDST.

Mycobacterial Strains
In total, 37 RR-MTB strains were collected among the approximately 845 positive cultures for M. tuberculosis examined for RIF resistance using proportion method. This includes all the new cases and pretreated patients' specimens. Based on the MICs results, 35 out of 37 examined strains were confirmed as rifampicin-resistant, and the remaining two isolates, which were susceptible to RIF, were excluded from further analysis. Out of the 35 RR-MTB strains included, one strain was isolated from an extrapulmonary abscess specimen and all others were from pulmonary TB patients.   Overall, 10 pre-extensively drug-resistant (pre-XDR) (MDR-TB also resistant to any FLQ or one of the second-line injectable drugs (SLIDs)) strains and three XDR-TB were found whose resistance pattern is displayed in Table 1. All the pre-XDR/XDR strains were susceptible to LZD, BDQ, and DLM, except two pre-XDR strains (KAN mono-resistant), which were resistant to DLM. Overall, 10 pre-extensively drug-resistant (pre-XDR) (MDR-TB also resistant to any FLQ or one of the second-line injectable drugs (SLIDs)) strains and three XDR-TB were found whose resistance pattern is displayed in Table 1. All the pre-XDR/XDR strains were susceptible to LZD, BDQ, and DLM, except two pre-XDR strains (KAN mono-resistant), which were resistant to DLM.

Genotypic Resistance to Fluoroquinolones
Five (83.3%) out of six LFX/MFX-resistant strains had mutations in the gyrA gene, including Gly88Cys (two strains), Asp94Gly (two strains), and Ser91Pro (one strain). However, an LFX/MFX susceptible strain (MICs on critical concentration for both drugs) showed a gyrA Asp94Ala mutation. There was no mutation in quinolone-resistant determining region (QRDR) of the gyrB gene (Table 3).

Genotypic Resistance to Second-Line Injectable Drugs
All mutations found in rrs, eis, and tlyA genes associated SLIDs resistance are shown in Table 4. The A1401G (eight strains) mutation was the dominant mutation in rrs gene, followed by G1484T (one strain). In addition to the rrs gene, four different mutations were found in eis promoter (C-10T, G-12A, G-14A, and G-15C), indicating low-level resistance to KAN. However, strain with a G-14A mutation, was phenotypically susceptible to KAN. In the case of CAP, a C insertion at position 210 of tlyA gene was found in a CAP-resistant strain along with eight strains having rrs mutations. Hence, 8/8 (100%), 11/11 (100%), and 9/9 (100%) strains showed resistant genotype for AMK, KAN, and CAP, respectively. Nevertheless, Gly196Glu mutation of tlyA gene was observed in a CAP-susceptible strain.

Genotypic Resistance to Other Second-Line Drugs (WHO Group C)
The mutations found in the genes associated with resistance to PTO, DCS, CFZ, and LZD are detailed in Table 5. Among the 19 PTO-resistant strains, mutations in ethA, ethA-ups, inhA, and fabG1-ups were detected in all (19/19, 100%) strains, of which 10 showed a fabG1-ups C-15T substitution. Furthermore, two susceptible strains had a ethA Gln254Pro mutation. All the nine DCS-resistant strains were wild type, except for one strain (1/9, 11.1%) that showed cycA Pro188Ala mutation. Additionally, 10 DCS-susceptible strains were found to have different mutations in targeted genes, of which seven showed mixed population. Moreover, six strains had MICs = 16 µg/mL, identical to the H37Rv strain ( Table 5). All the genes involved in resistance to CFZ were wild type in CFZ-resistant strains, except Rv1979c, which showed adenine deletion at position 1014 (1/5, 20%). However, the MICs of three strains with various substitutions were on breakpoint (1 µg/mL). Two strains with rrl mutations (C1537T and C1331T) and two with a rplC Ala72Thr mutation were susceptible to LZD (Table 5).

Genotypic Resistance to Add-On Agents
The only mutation detected in BDQ resistance-related genes was Rv0678 Gly87Arg, which was found in a BDQ-susceptible strain. The atpE, atpE-ups (upstream), and Rv0678-ups genes did not harbor any mutations in any strain. A BDQ-resistant strain (MICs = 0.5 µg/mL) was found as wild type in studied genes (0/1, 0%). Three different mutations (G-deletion at position 91 of ddn, fbiA Arg321Ser (two strains) and fbiC Trp678Gly) were found in 4/9 (44.4%) of DLM-resistant strains. However, fbiC Trp678Gly mutation was also identified in a DLM-susceptible strain.

Sensitivity, Specificity, and Accuracy of WGS
WGS showed the highest sensitivity in detecting resistance to EMB, PTO, and SLIDs (100%). Additionally, the lowest sensitivity was found in the case of BDQ, DCS, and CFZ by 0.0%, 11.1%, and 20.0%, respectively. Table 6 summarizes the sensitivity, specificity, and overall accuracy of WGS for all tested drugs, compared with pDST.

WGS-Based MTB Genotyping
WGS confirmed the heterogeneous TB epidemiology among our isolates with the four major MTB complex lineages. Lineage 1; East-African Indian (includes EAI and EAI Manila), Lineage 2; East-Asian (includes Beijing), Lineage 3; Indo-Oceanic (includes Delhi-CAS), and Lineage 4; Euro-American (includes Cameroon, Haarlem, H37Rv-like, LAM, mainly T, S-type, Ural and X-type). Most represented sublineages were Beijing, Euro-American (sublineage 4.5), and Delhi-CAS (Figure 2). Cluster analysis using the threshold of five SNPs as maximum distance between group members identified three distinct clusters (A, B, and C). Four pre-XDR strains  from three provinces from west to the capital of Iran classified in cluster A (Beijing). Cluster B (Ural) includes two MDR-MTB strains (154-19, 173-19) recovered from two patients in the capital. Strains in cluster A and B were isolated during several years and might distributed due to referring drug-resistant patients from other cities to the capital for treatment and following. Two pre-XDR/XDR strains (160-19, 176-19) isolated from members of a family during a year were placed in cluster C (Figure 2). The genotype of five remaining pre-XDR strains were Beijing (two strains from the west), Euro-American (one strain from the capital and one strain from the Afghani patient at south), and LAM (1 strain, Capital). An XDR-TB patient was from the north (Beijing) and another XDR strain isolated from the Iraqi patient at south-west of Iran (LAM). Among the RR-MTB strains, three Euro-American, one EAI, and one Delhi-CAS lineages/sublineages were detected in patients from the capital.

Discussion
To determine the MIC distributions of first-/second-line anti-TB drugs among RR/MDR-MTB strains, phenotypic DST was performed using the broth microdilution method. Different methods are available to determine the MICs and based on the results and some other factors including, cost, labor requirements, published MIC data, and experiences from other pathogens, Middlebrook 7H9 was also proposed as the EUCAST reference method for MTBC [34]. Additionally, WGS analysis was performed on 35 collected RR-MTB strains to screen genetic mutations in drug-resistance-related genes.
In the case of rifampicin, accuracy of the WGS to detect rifampicin-resistant strains was 97.1%, and regardless of strains with mixed mutations, all the mutations were located at 81-bp rifampicinresistant determinant region (RRDR) of the rpoB gene, except for one (Ile572Phe). Thus, screening only the 81-bp region could detect 94.3% of rifampicin-resistant strains that some resistant strains with no mutation in this region have previously shown [35,36]. In addition, there were no strains having mixed mutations within the RRDR. In accordance with the previous studies from different geographical regions [37][38][39], codon 531 of the rpoB gene was the most abundant mutation in this study. Moreover, 85.7% of the RR-MTB strains were found to be resistant to INH (MDR-TB), and testing only the katG gene could detect 83.3% of the strains. However, 83.3% of INH-resistant strains had katG Ser315Thr and fabG1-inhA C-15T mutations which could have been increased to 96.6% if the sequence analysis included the entire katG gene. These data were supported by other studies [40][41][42][43][44].

Discussion
To determine the MIC distributions of first-/second-line anti-TB drugs among RR/MDR-MTB strains, phenotypic DST was performed using the broth microdilution method. Different methods are available to determine the MICs and based on the results and some other factors including, cost, labor requirements, published MIC data, and experiences from other pathogens, Middlebrook 7H9 was also proposed as the EUCAST reference method for MTBC [34]. Additionally, WGS analysis was performed on 35 collected RR-MTB strains to screen genetic mutations in drug-resistance-related genes.
In the case of rifampicin, accuracy of the WGS to detect rifampicin-resistant strains was 97.1%, and regardless of strains with mixed mutations, all the mutations were located at 81-bp rifampicin-resistant determinant region (RRDR) of the rpoB gene, except for one (Ile572Phe). Thus, screening only the 81-bp region could detect 94.3% of rifampicin-resistant strains that some resistant strains with no mutation in this region have previously shown [35,36]. In addition, there were no strains having mixed mutations within the RRDR. In accordance with the previous studies from different geographical regions [37][38][39], codon 531 of the rpoB gene was the most abundant mutation in this study. Moreover, 85.7% of the RR-MTB strains were found to be resistant to INH (MDR-TB), and testing only the katG gene could detect 83.3% of the strains. However, 83.3% of INH-resistant strains had katG Ser315Thr and fabG1-inhA C-15T mutations which could have been increased to 96.6% if the sequence analysis included the entire katG gene. These data were supported by other studies [40][41][42][43][44]. The ahpC, and furA structural genes were wild type in strains examined in this study and had a less role in resistance to INH, as limited data and mutations with lower incidence of these genes have been reported in previous studies [42,43,45,46]. Ethambutol resistance was mostly linked to the mutations within the embB gene (codon 306 as the dominant position), but not embA. Khosravi et al. reported embA mutation from EMB-susceptible strains [47]. Additionally, Sun et al revealed that 82.6% of EMB-resistant isolates from China had at least a mutation in the embABC locus, but contrary to this study, they were located mostly in the embA and embB upstream region [48].
Sequencing results of gyrA gene covered five FLQ resistance (three XDR and two pre-XDR strains) with the accuracy of 94.3%, and the QRDR of the gyrB gene was wild type in all strains. The gyrA Asp94Ala, detected in an LFX/MFX susceptible strain with MICs equal to 1 and 0.25 µg/mL. This mutation may lead to only ciprofloxacin resistance or is associated with lower MICs to FLQ than other substitutions within this codon, as was reported in the literature [26,49]. In addition, gyrA Gly88Cys was more associated with higher MICs for both LFX and MFX than other detected mutations. A strain with low-level resistance to LFX (MIC = 2 µg/mL) and MFX (MIC = 0.5 µg/mL) did not harbor any mutation in the related genes, suggesting additional resistance mechanisms, such as efflux pumps [50].
Testing for only the rrs gene could detect all eight AMK/KAN/CAP-R phenotypes, including three XDR and five SLIDs resistant (pre-XDR) strains and resistance to KAN, and CAP was correctly diagnosed in all strains using rrs, eis-ups, and tlyA genes. In concordance with previous studies, two KAN mono-resistant strains (low-level resistance with MIC = 10 µg/mL) were found to harbor C-10T and G-15C substitutions at the eis promoter [51][52][53]. However, Pholwat et al. indicated that the G-15C mutation does not confer KAN resistance in M. tuberculosis [54]. A KAN/CAP-resistant strain, with MIC = 5 µg/mL for both drugs, possesses eis G-12A and tlyA C210 insertion. The only strain with a G-14A mutation at the eis promoter was phenotypically susceptible to KAN (MIC = 0.3 µg/mL). This mutation was investigated as KAN-resistance conferring mutation in literature [51,52]. Furthermore, tlyA Gly196Glu was found in a CAP-susceptible strain that may not determine CAP resistance, while Gly196Glu was previously reported in resistant isolates [30].
Approximately 54% of rifampicin-resistant strains were identified as PTO-resistant in the present study. Additionally, the sensitivity and specificity of PTO resistance detection were 100% and 87.5%, respectively, using WGS. Cross resistance among PTO and INH was noted in 94.7% of PTO-resistant strains. Similar to the previous studies, in this study, more than half of the total mutations in PTO-resistant phenotypes were found in the inhA/fabG1 regulatory regions [55,56]. The ethA Gln254Pro mutation was found in two susceptible phenotypes, while the later mutation was also identified in two resistant strains (MIC = 5, 10 µg/mL). Nevertheless, the Gln254Pro substitution was previously reported from ethionamide-resistant M. tuberculosis strains [57]. Surprisingly, DCS resistance-related genes targeted in this study showed a single cycA Pro188Ala mutation in one of nine strains recognized as DCS-resistant (the low sensitivity by 11.1%). Despite this, 10 out of 26 DCS-susceptible strains harbored different mutations and most of them, with MICs equal to or less than what was observed in the H37Rv, showed mixed wild type/mutant population. Results of this study emphasized the controversial role of alr, ald, ddlA, and cycA genes in molecular detection of DCS resistance. Previous attempts indicated complex mechanisms in resistance to this drug [58,59]. Only the alr D344N substitution, identified by Chen et al., and novel mutations in genes other than known DCS resistance-associated genes have been reported [59]. WGS, with low sensitivity (20%), traced only A1014 deletion of Rv1979c among five CFZ-resistant strains. The reason for this low sensitivity is that currently, BDQ and CFZ do not yet have well-defined mechanism of resistance or well-documented genomic targets [60]. However, three strains, with MICs on breakpoint (MIC = 1 µg/mL), had various mutations in three candidate genes, including Rv0678, Rv1979c, and pepQ. These results were inconsistent with those of the previous studies. Zhang et al. showed that 97% of the CFZ-resistant mutants had a mutation in the Rv0678 gene, encoding the MmpR repressor [61]. Additionally, mutations in these three genes and their upstream region were reported in the CFZ-resistant strains [32,61].
In the present study, LZD and BDQ were the most effective drugs against MDR/pre-XDR and XDR strains of M. tuberculosis. However, unlike the present data, 10% LZD resistance in MDR/XDR isolates was previously reported from Iran [62]. Different mutations in rrl and rplC genes had a role in LZD resistance [32], while the four mutant strains in the current study, with C1537T and C1331T substitutions in the rrl gene and rplC Ala72Thr mutation (two strains), were LZD-susceptible. In addition, Yang et al. reported several rrl mutations in LZD-susceptible strains with MIC below 0.5 mg/L [63]. There was no cross resistance between CFZ and BDQ in the tested strains in this study, and the only Rv0678 Gly87Arg mutant was also susceptible to BDQ (MIC = 0.01 µg/mL). The atpE gene was wild type in all strains. Meanwhile, the strain with MIC greater than the epidemiological cut-off (0.12 µg/mL, reported by the European Committee on Antimicrobial Susceptibility Testing-EUCAST) were genetically wild type. Thus, the role of the Rv0678 and atpE in resistance is controversial, as mutations of these genes in susceptible and resistant strains have been reported. In other words, both susceptible and BDQ-resistant strains may carry mutations in these genes and the presence of mutations does not necessarily indicate resistance [32,63,64].
Among the rifampicin-resistant strains, 25.7% were resistant to DLM. In other words, three rifampicin-resistant, four MDR, and two pre-XDR strains were resistant to DLM and all three XDR strains were susceptible to DLM. In addition, DLM-resistant strains were belonged to different lineages/sublineages including Beijing (three strains), Euro-American (three strains), Delhi-Cas (three strains). Additionally, six of them were isolated from same city and the remaining three strains obtained from two different cities. As the first investigation, the high rate of DLM resistance was amazingly considerable, while this new compound has not yet been used as a treatment choice of TB patients in Iran. Moreover, the sensitivity of WGS in DLM resistance detection was almost low (44.4%). The fbiA Arg321Ser was the common mutation found in two (22.2%) DLM-resistant strains. However, fbiC Trp678Gly mutant was found in an DLM-susceptible phenotype. Molecular basis of resistance to DLM is still poorly understood, and the present results around the DLM MICs and genetic background of its resistance-associated genes are in accordance with previous studies [33,63].

Conclusions
The highest resistance rate among RR/MDR strains of M. tuberculosis was found in PTO. However, LZD and BDQ were found as the most effective drugs against all kinds of drug-resistant phenotypes (RR/MDR/pre-XDR/XDR). It was surprising to find more than 25% DLM resistance among the drug-resistant phenotypes without exposure to this agent. WGS showed the lowest sensitivity for detecting resistance to BDQ (0.0%), DCS (11.1%), and CFZ (20.0%), followed by DLM (44.4%) which is related to the current poor knowledge of genetic variants related to the resistant phenotype. Additionally, the highest sensitivity among the second-line drugs belonged to PTO and SLIDs (100.0%). Finally, mutations in drug resistance-related genes were mostly linked to the MICs results, especially for key anti-TB drugs; hence, the RR/MDR/pre-XDR/XDR strains can be identified with high confidence using WGS. However, if fully susceptible strains were present in this study, a better interpretation of the mutations could be achieved. In addition, limited data around the molecular background and resistance mechanisms of some drugs, especially new compounds, was another restriction.