Discordance Between Phenotypic and WGS-Based Drug Susceptibility Testing Results for Some Anti-Tuberculosis Drugs: A Snapshot Study of Paired Mycobacterium tuberculosis Isolates with Small Genetic Distance

Introduction

Tuberculosis (TB) is a devastating public health threat worldwide, and drug resistance of the causative agent, Mycobacterium tuberculosis (Mtb), is a significant challenge for successful disease control. Drug-resistant TB treatment courses are longer (up to 20 months), regimens require a higher pill burden, patients have poorer treatment outcomes, and second-line medications are associated with higher toxicity and more serious adverse events than drug-susceptible TB treatments.¹ Notably, Mtb lineages vary in geographical distribution, virulence, immunogenicity,^2,3 and propensity to acquire drug resistance,^4–6 highlighting the importance of local studies targeting Mtb strain-level genetic diversity and drug resistance for proper population risk assessment.

According to the World Health Organization (WHO), the global estimated total TB incidence reached 10.6 million people by 2022, and 410 thousand TB cases were multidrug-resistant (MDR) or rifampicin-resistant (RR).⁷ Latvia is a WHO-priority country in Europe with a low-to-moderate TB incidence.⁸ The estimated total incidence was 19 TB cases per 100,000 population in 2022, and among them, 2.1 cases were MDR/RR-TB.⁷ Multiple Mtb genotypes belonging to Lineages 2 and 4, which are prevalent in Latvia and neighboring countries, have been associated with drug resistance,⁹ whereas widespread MDR and pre-extensively drug-resistant Mtb strains have been an important healthcare issue in Latvia for decades.¹⁰

The adjustment of the TB treatment regimen relies on the determination of the Mtb drug resistance pattern, and this step is crucial not only for rational and effective therapy prescription but also for the prevention of further transmission of drug-resistant Mtb strains.¹¹ Phenotypic drug susceptibility testing (pDST) is a widely used culture-based approach in drug-resistant TB diagnostics; however, owing to the slow growth rate of Mtb, it typically takes at least two weeks on liquid media and approximately a month on solid media before the results can be interpreted.¹² Furthermore, for some anti-TB drugs, pDST is unreliable; therefore, routine testing is not recommended.¹³ Rapid molecular tests for the identification of Mtb resistance-conferring variants, including automated nucleic acid amplification tests such as Xpert MTB/RIF Ultra (Cepheid, Sunnyvale, USA), and line probe assays such as GenoType MTBDRplus (Hain Lifesciences-Bruker, Nehren, Germany), can be performed directly on sputum specimens and produce results within two–five hours.^14,15 Rapid molecular tests to identify rifampicin (RIF), isoniazid (INH), pyrazinamide (PZA), fluoroquinolone (FQ), ethionamide (ETO), and amikacin (AMK) resistance-associated variants are available.¹⁶ Although these tests have improved drug-resistant TB diagnostics¹⁷ and demonstrated high accuracy in drug resistance detection, especially for the first-line drugs RIF and INH,^18,19 there is still a risk of bias as these assays target only a limited number of genomic variants in the selected Mtb genes.²⁰

Recently, whole-genome sequencing (WGS) has become the most common approach in Mtb studies, as it not only allows genome-wide analysis of resistance-related variants, including single nucleotide variants (SNVs), insertions, and deletions, but is also widely used for Mtb genotyping, transmission chain investigation, and recurrent TB cause determination.^21,22 Currently, WGS-based prediction of Mtb drug resistance patterns mostly uses publicly available automatic bioinformatic tools^23,24 and has demonstrated promising results and great potential for implementation in routine clinical laboratory practice. However, discordant results between phenotypic and WGS-based DST are still being reported, even when the association between the genomic variant and Mtb drug resistance is well known.^25–31 Furthermore, the association of many variants previously detected in drug-resistant Mtb strains with the resistant phenotype remains uncertain.³²

Herein, we used the TB-Profiler tool³³ and the current comprehensive catalog of Mtb variants and their association with drug resistance developed by the WHO³² to retrospectively investigate drug-resistance-conferring variants in the WGS data of Mtb isolate pairs characterized by a small genetic distance of ≤10 SNVs. Isolate pairs were obtained from Latvian pulmonary TB patients who either suffered from recurrent episodes of active infection or had difficulties in achieving successful treatment outcomes during a single TB episode. Subsequently, WGS-predicted drug resistance profiles were compared with the pDST data, and the consistency of the DST results for closely related paired Mtb isolates was assessed. In addition, we analyzed the distribution of drug-resistance-conferring variants among the identified Mtb genotypes.

Materials and Methods

Study Population, Mycobacterial DNA Sample Collection, and Spoligotyping

This study aimed to retrospectively assess the agreement between phenotypic and WGS-based DST results for longitudinal single-episode and recurrent TB patient Mtb isolate pairs with small genetic distances of ≤10 SNVs, and to characterize the drug-resistance-conferring variant distribution among identified Mtb genotypes. This study included patients with pulmonary TB admitted to the Centre of Tuberculosis and Lung Diseases, Riga East University Hospital in 2002–2019, which is the nationwide center of TB diagnostics and treatment. Two Mtb isolates were obtained from all patients and either represented each of the active TB episodes or both were acquired during a single disease episode within at least two-month intervals.

Cultures were grown on Löwenstein–Jensen (LJ) media at the clinical laboratory of the Centre of Tuberculosis and Lung Diseases during the patients’ stay at the hospital, and the harvested bacteria were used for DNA extraction according to the cetyltrimethylammonium bromide protocol.³⁴ Spoligotyping was conducted using commercially available reagent kits (Isogen Life Science, Netherlands; later – Ocimum Biosolutions, India) following a previously published protocol.³⁵ Spoligotype international types (SIT) were inferred using the SITVIT2 database (available at http://www.pasteur-guadeloupe.fr:8081/SITVIT2/). Mtb isolate pairs with identical spoligotypes were divided into two cohorts (recurrent TB episode isolates and longitudinal single TB episode isolates), and genetic distances between paired isolates and causes of recurrent TB episodes were determined as described previously.³⁶

Phenotypic Drug Susceptibility Testing

The phenotypic drug resistance profiles of all Mtb isolates were compiled based on the available pDST reports from hospital medical records (Supplementary Table 1). pDST was conducted using either the Bactec Mycobacterial Growth Indicator Tube (MGIT) 960 system (Becton Dickinson, Heidelberg, Germany) or LJ media according to the updated WHO guidelines, which were effective at the time of test performance. Data on capreomycin and kanamycin were excluded from the analysis because they were no longer approved for TB treatment,¹ whereas ofloxacin (OFX), moxifloxacin (MFX), and levofloxacin (LFX) were analyzed as a single FQ group.

The critical concentrations (CC) for Bactec MGIT 960 system (liquid media) were 1 µg/mL (years 2005–2019) and 2 µg/mL (years 2002–2004) for RIF, 0.1 µg/mL for INH, 5 µg/mL (years 2005–2019) and 7.5 µg/mL (years 2002–2004) for ethambutol (EMB), 100 µg/mL for PZA, 2 µg/mL (years 2002–2016) for OFX, 0.25 µg/mL (years 2009–2019), 0.5 µg/mL (years 2015–2016), 1 µg/mL (years 2018–2019), and 2 µg/mL (years 2014–2016) for MFX, 1 µg/mL (years 2018–2019), and 2 µg/mL (years 2014–2016) for LFX, 1 µg/mL (years 2005–2014) and 6 µg/mL (years 2002–2004) for streptomycin (STR), 1 µg/mL for AMK (years 2008–2019), 5 µg/mL (years 2004–2005) for ETO, 4 µg/mL (year 2004) for para-aminosalicylic acid (PAS), 1 µg/mL (years 2014–2019) for linezolid (LZD), 1 µg/mL (years 2018–2019) for bedaquiline (BDQ), 1 µg/mL (year 2019) for clofazimine (CFZ), and 0.06 µg/mL (year 2019) for delamanid (DLM).

The CCs for LJ media were 40 µg/mL (years 2002–2016) for RIF, 0.2 and 1 µg/mL (years 2003–2016) for INH, 2 µg/mL (years 2003–2016) for EMB, 4 µg/mL (years 2008–2016) for OFX, 4 µg/mL (years 2003–2014) for STR, 30 µg/mL (years 2008–2018) for AMK, 30 µg/mL (years 2003–2008), and 40 µg/mL (years 2008–2016) for ETO, 0.5 µg/mL (years 2003–2008), and 1 µg/mL (years 2008–2019) for PAS, and 30 µg/mL (years 2012–2015) for D-cycloserine (DCS).

For most Mtb isolates, pDST was performed using both the methods. In these cases, the results acquired from testing on both solid and liquid media were compared and combined to determine the phenotypic drug resistance profiles. If the results were mismatched and the Mtb isolate demonstrated resistance to any medication in only one of the tests performed, the isolate was considered drug-resistant.

WGS and Bioinformatic Data Processing

The WGS library preparation process and bioinformatics pipeline have been reported previously.³⁶ Sequence alignment files (BAM) generated by the snippy tool (v3.6.0) were used as the TB-Profiler software inputs (v4.1.1) to identify Mtb sub-lineages and detect genetic variants in mycobacterial genes associated with Mtb drug resistance. Variants with allele frequencies ≥10% and minimum coverage depths of four sequencing reads per base were called. In addition, the snippy tool outputs containing annotated genetic variants detected in studied Mtb isolates were used to identify variants in those mycobacterial genes that are mentioned in the catalog of Mtb variants and their association with drug resistance developed by WHO but are not yet included in the TB-Profiler tool algorithm. The significance of all detected variants was assessed using the WHO catalog of Mtb variants and their association with drug resistance.

Amplification of the pncA Gene

In the case of the three isolate pairs, amplification of the Mtb pncA gene was necessary for complete drug resistance profile acquisition. The polymerase chain reaction (PCR) was performed according to a previously published protocol.³⁷ Primers pncA-F 5′-AACAGTTCATCCCGGTTC-3′ and pncA-R 5′-GCGTCATGGACCCTATATC-3′, targeting the entire pncA gene (product length: 668 bp), were used. The PCR mixture (26 µL) was prepared as follows: 13 µL of 2x DreamTaq PCR Master Mix (Thermo Fisher Scientific, US), 2 µL of DNA template, and 0.2 µM of each primer. PCR was performed under the following thermal cycling conditions: denaturation at 95°C for 3 min, followed by 35 amplification cycles (30s at 95°C, 30s at 58°C, and 1 min at 72°C) and an elongation step at 72°C for 5 min. The PCR products were visualized on a 1.5% agarose gel (TopVision Agarose, Thermo Fisher Scientific, US). Sanger sequencing was performed using a BrilliantDye Terminator (v3.1) Cycle Sequencing Kit (NimaGen, Netherlands) according to the manufacturer’s protocol using an ABI Prism 3100 Genetic Analyzer (PerkinElmer, USA).

Results

In total, this study included 46 Mtb isolate pairs (92 Mtb isolates altogether) that had pairwise SNV-distances ≤10; in 30 cases, isolate pairs were obtained from patients who experienced recurrent TB episodes, and in 16 cases, from patients with single TB episodes. Cases were described, and pairwise SNV distances between paired Mtb isolates were calculated as previously described.³⁶ Here, we performed a comparative analysis of the phenotypic and WGS-predicted drug resistance profiles of 46 Mtb isolate pairs. WGS-predicted drug resistance profiles for paired Mtb isolates were identical, except for four cases of FQ resistance development in the subsequent isolates and one case of a heterozygous PZA-resistant allele identified only in the initial isolate of the phenotypically PZA-resistant Mtb isolate pair. All drug-resistance-conferring variants among the studied Mtb isolate pairs were either detected by TB-Profiler software or identified during visual inspection of pncA gene PCR amplicon sequences and included in Supplementary Table 2.

According to WGS, no resistance-associated variants were detected in 21 cases (45.7%), and the pDST data confirmed the drug susceptibility of these isolate pairs. Sub-lineages 4.3.3 (including SIT254 (n=6), and SIT42 (n=2) spoligotypes), 4.8 (including SIT53 (n=4), SIT40 (n=1), and SIT156 (n=1) spoligotypes), and 4.1.2.1 (SIT47; n=1) comprised by the Lineage 4, and sub-lineage 2.2.1 (SIT1; n=6) belonging to the Lineage 2 were identified among drug-susceptible Mtb isolate pairs. In 25 (54.3%) cases, WGS analysis results indicated the presence of resistance-conferring variants. All Mtb isolate pairs resistant to at least one anti-TB medication belonged to Lineages 2 and 4 and represented sub-lineages 2.2.1 (including SIT1 (n=12), and SIT190 (n=1) spoligotypes), 4.1.2.1 (SIT283; n=1), 4.2.1 (SIT1117; n=1), 4.3.3 (SIT42; n=7), and 4.8 (SIT53; n=3). As in most cases, multiple spoligotypes representing the same TB-Profiler-assigned sub-lineage were determined among studied Mtb isolate pairs, spoligotyping demonstrated greater resolution in genotype determination than TB-Profiler software for this dataset. Thus, we used spoligotyping-assigned genotypes for further resistance-conferring variant distribution characteristics.

Resistance-conferring variants of RIF, INH, EMB, PZA, FQ, AMK, STR, ETO, and PAS were detected in drug-resistant Mtb isolate pairs (Figure 1). All the WGS-detected variants were homozygous (allele frequency 90–100%). Further detailed analysis of the detected drug-resistance-conferring variants and comparison of the WGS-predicted and phenotypic drug resistance data are provided below.

Figure 1 Comparison of phenotypic and WGS-based drug susceptibility testing results of drug-resistant M. tuberculosis isolates. Abbreviations: AMK, amikacin; BDQ, bedaquiline; CFZ, clofazimine; DCS, D-cycloserine; DLM, delamanid; DST, drug susceptibility testing; EMB, ethambutol; ETO, ethionamide; FQ, fluoroquinolones; INH, isoniazid; LZD, linezolid; N, drug susceptibility testing results were not available; PAS, para-aminosalicylic acid; pDST, phenotypic drug susceptibility testing; PZA, pyrazinamide; R, resistant; RIF, rifampicin; S, sensitive; STR, streptomycin.

Rifampicin

The phenotypic and WGS-based RIF DST results were identical for 95.7% (44/46) of Mtb isolate pairs and 97.8% (90/92) of the studied Mtb isolates (Table 1). According to WGS-based DST results, 18 Mtb isolate pairs were resistant to RIF. Variants in the RIF resistance-determining region³⁸ of locus Rv0667 (rpoB gene) were detected in all cases. Among them, the Ser450Leu variant was the most prevalent; it was found in 10 Mtb isolate pairs of three different spoligotypes (SIT1, SIT42, and SIT190), and in two more isolate pairs belonging to the SIT1 spoligotype it coincided the Gly332Arg variant at locus Rv0668 (rpoC gene). His445Tyr and Asp435Val were detected in two and three Mtb isolate pairs of the SIT1 spoligotype, respectively.

Table 1 Rifampicin Resistance-Conferring Variants, and Comparison with Phenotypic Drug Susceptibility Testing Results

Two inconsistencies were found between the phenotypic and WGS-based DST. In the SIT1117 isolate pair, only the initial Mtb isolate was phenotypically resistant, although both isolates harbored Asp435Tyr in rpoB gene. While no RIF-resistance-conferring variants were detected in the SIT283 isolate pair, the initial isolate demonstrated phenotypic resistance.

Isoniazid

The phenotypic and WGS-based INH DST results were concordant for 97.8% (45/46) of Mtb isolate pairs and 98.9% (91/92) of the studied Mtb isolates (Table 2). Only three resistance-associated variants were detected among the 24 Mtb isolate pairs. In 11 isolate pairs belonging to SIT1, SIT190, SIT283, and SIT1117 spoligotypes, INH resistance occurred because of the Ser315Thr variant at the locus Rv1908c (katG gene). The C-15T variant at locus Rv1483 (fabG1 gene), which is also associated with ETO cross-resistance, was the only INH resistance-related variant in all three SIT53 isolate pairs; however, in one case, pDST results matched WGS-based DST findings only for the initial isolate in the pair. Furthermore, three SIT1 and six INH-resistant SIT42 isolate pairs harbored both Ser315Thr and C-15T variants. In one more phenotypically INH-resistant isolate pair belonging to the SIT1 spoligotype, the Lys143Glu variant in the katG gene was detected. This variant was not included in the TB-Profiler database as resistance-conferring; however, it was described as of uncertain significance (grade group 3) by the WHO.

Table 2 Isoniazid Resistance-Conferring Variants, and Comparison with Phenotypic Drug Susceptibility Testing Results

Ethambutol

For 93.5% (43/46) of the studied Mtb isolate pairs and 95.7% (88/92) of isolates the phenotypic and WGS-based EMB DST results were concordant (Table 3). EMB-resistance-conferring variants were detected in 17 Mtb isolate pairs, and variants at the locus Rv3795 (embB gene) were detected in all but two isolate pairs. Met306Val was found in six SIT1 isolate pairs; however, in one isolate pair, the initial isolate was phenotypically sensitive. Met306Ile was identified in one SIT190 isolate pair and in two SIT1 Mtb isolate pairs this variant coincided with Asn296His. According to the WHO variant database, the association between the Asn296His variant and EMB resistance remains controversial. Similarly, G-43C at locus Rv3794 (embA gene) is also considered a variant of uncertain significance in causing EMB resistance, and it has not been included in the TB-Profiler database. However, in our study, G-43C was only identified in one phenotypically resistant SIT1 isolate pair. EMB resistance in SIT42 isolates was caused by Tyr319Ser and Gln497Arg variants in the embB gene. In one case, Gln497Arg was found in combination with the C-12T variant at the locus Rv3794 (embA gene). Another SIT42 isolate pair harbored the Val7Gly variant at the locus Rv1267c (embR gene); however, the pDST results did not reveal EMB resistance in either isolate. The association between this variant and EMB resistance remains to be clarified. In addition, the initial isolate of the SIT1117 pair was phenotypically EMB-resistant and no resistance-conferring variants were detected.

Table 3 Ethambutol Resistance-Conferring Variants, and Comparison with Phenotypic Drug Susceptibility Testing Results

Pyrazinamide

The PZA pDST data were available for 17 pairs of Mtb isolates and 40 Mtb isolates together, while for 52 Mtb isolates from 29 pairs, pDST results were missing. Comparison of WGS-based and pDST results demonstrated 81.2% (13/16) concordance for Mtb isolate pairs and 85% (34/40) concordance for the studied isolates (Table 4). According to the WGS data, 10 isolate pairs harbored PZA-resistance-conferring variants at the locus Rv2043c (pncA gene). In the case of three SIT1 Mtb isolate pairs, TB-Profiler detected large deletions, including either the whole pncA gene or a part of the coding sequence; specifically, 2285819–2291116del was identified in two Mtb isolate pairs and 2289105–2289321del in one isolate pair. However, visual exploration of the WGS data showed that this genomic region lacked sequencing coverage in these isolates (Supplementary Figure); therefore, PCR-based analysis was performed.

Table 4 Pyrazinamide Resistance-Conferring Variants, and Comparison with Phenotypic Drug Susceptibility Testing Results

PCR amplification of the pncA gene in these Mtb isolates was successful (Figure 2). The Thr76Pro variant was identified in the first isolate pair, however, in the subsequent isolate this variant was heterozygous as the presence of both reference and alternate allele in this position was detected (Figure 2B); unfortunately, pDST data were not available for this isolate. In the second Mtb isolate pair, the heterozygous Leu35Pro variant was detected only in the initial isolate (Figure 2C), although both the isolates were phenotypically resistant. No resistance-conferring variants were identified in the third Mtb isolate pair that was phenotypically PZA-resistant.

Figure 2 Pyrazinamide-resistance-conferring variants detected by PCR-based analysis in two M. tuberculosis isolate pairs. The picture shows a fragment of sequenced pncA gene PCR amplicons for M. tuberculosis isolate pairs, whose WGS data lacked pncA gene covering sequencing reads. The initial isolate from patient P1SIT1.1 (A) harbored the homozygous resistance-conferring variant, whereas the subsequent isolate (B) harbored both resistant and susceptible gene alleles. For patient P27SIT1.12, the heterozygous resistance-conferring variant was detected only in the initial isolate (C), and the subsequent isolate did not harbor any variants in the pncA gene (D).

PZA-resistance-conferring variants were detected in five SIT1 isolate pairs: Leu35Pro was detected in one more isolate pair, but pDST data were available only for the initial isolate; Cys72Trp was detected in one isolate pair with no pDST data available; and Tyr103His and Leu85Arg were identified in one and two phenotypically resistant Mtb isolate pairs, respectively. Both Leu35Pro and Cys72Trp are listed in the TB-Profiler database as resistance-conferring variants; however, according to WHO, their significance remains uncertain. The Ser59Pro variant was detected in the SIT190 isolate pair, but the pDST data were available only for the initial isolate. Asp63Gly was identified in four isolate pairs belonging to SIT42. However, in one isolate pair, only the initial isolate was phenotypically resistant. Pro62Thr was detected in the phenotypically susceptible SIT42 isolate pair. According to WHO, this variant is currently considered to be associated with PZA resistance; however, it was not included in the TB-Profiler database.

Fluoroquinolones

The results of FQ pDST were available for 22 pairs of Mtb isolates and 60 Mtb isolates together, while for 32 Mtb isolates from 24 pairs, pDST results were missing. The WGS-based DST findings corresponded with the pDST results (Table 5). All resistance-associated variants were found at the locus Rv0006 (gyrA gene). FQ-resistance-conferring variants were detected in five Mtb isolate pairs of SIT1, SIT190, and SIT42 spoligotypes. Furthermore, in four isolate pairs of SIT1 and SIT42, resistance to FQ developed in the subsequent isolates. In three Mtb isolate pairs belonging to SIT1 and SIT42 as well as in three subsequent isolates of the SIT1 spoligotype, resistance was caused by the Asp94Gly variant. Asp94Ala and Ala90Val were detected in the SIT1 and SIT190 isolate pairs, respectively, whereas Gly88Cys was identified in the subsequent isolate of the SIT42 isolate pair.

Table 5 Fluoroquinolone Resistance-Conferring Variants, and Comparison with Phenotypic Drug Susceptibility Testing Results

Streptomycin

The results of STR pDST were available for 24 pairs of Mtb isolates and 68 isolates together, whereas for 24 Mtb isolates from 22 pairs, pDST results were missing. The WGS-based DST findings matched the pDST results for the phenotypically tested Mtb isolates (Table 6). Variants associated with STR resistance were detected in the 24 Mtb isolate pairs. The Lys43Arg variant at locus Rv0682 (rpsL gene) was detected in the SIT190 isolate pair and in nine SIT1 isolate pairs. Two variants were identified in the rrs region of nine isolate pairs: C517T was detected in two SIT1 isolate pairs, whereas A514C was detected in one SIT1 and six SIT42 isolate pairs. All variants harboring Mtb isolates were phenotypically resistant, except for one Mtb isolate pair and subsequent isolates of six more isolate pairs when pDST data were not available. Three isolate pairs belonging to SIT53, SIT1117, and SIT283 harbored Gly34Glu, Leu59Arg, and Ala138Val variants at locus Rv3919c (gid gene), respectively. The pDST results supported the WGS findings for all but the subsequent SIT283 isolate when pDST data were missing. None of the detected gid gene variants were included in the TB-Profiler database; however, according to WHO, all of them had been detected in the STR-resistant Mtb isolates previously, but these variants are yet of uncertain significance.

Table 6 Streptomycin Resistance-Conferring Variants, and Comparison with Phenotypic Drug Susceptibility Testing Results

Amikacin

In the case of AMK, pDST results were available for 13 pairs of Mtb isolates and 32 isolates together, whereas for 60 Mtb isolates from 33 pairs, pDST results were missing. Phenotypic and WGS-based DST results were identical in 92.3% (12/13) of Mtb isolate pairs and 96.9% (31/32) of the studied isolates (Table 7). A1401G in the rrs region was the only AMK-resistance-conferring variant detected in seven Mtb isolate pairs belonging to the SIT1 and SIT42 spoligotypes. Phenotypic DST results were available for all but four of these isolates, and they matched the WGS-based DST findings. No resistance-related variants were detected in the initial isolate of one pair of SIT1 isolates that exhibited phenotypic resistance to AMK.

Table 7 Amikacin Resistance-Conferring Variants, and Comparison with Phenotypic Drug Susceptibility Testing Results

Ethionamide

The results of ETO pDST were available for 16 pairs of Mtb isolates and 35 Mtb isolates together, while for 57 Mtb isolates from 60 pairs, pDST results were missing. Concordance between phenotypic and WGS-based DST was only 50% (8/16) for Mtb isolate pairs and 65.7% (23/35) for Mtb isolates (Table 8). Variants associated with ETO resistance were detected in 19 isolate pairs. Twelve isolate pairs belonging to the SIT1, SIT42, and SIT53 spoligotypes were ETO-resistant because of the aforementioned C-15T variant at locus Rv1483 (fabG1 gene). Unfortunately, pDST data were not available for one subsequent SIT1 isolate, one SIT42 and all three SIT53 Mtb isolate pairs. Two SIT42 isolate pairs and two subsequent SIT1 and SIT42 isolates harboring this variant were phenotypically resistant. In two more phenotypically resistant isolate pairs belonging to the SIT42, the Ile338Ser variant at locus Rv3854c (ethA gene) was simultaneously detected, and its association with Mtb ETO resistance has not yet been clarified.

Table 8 Ethionamide Resistance-Conferring Variants, and Comparison with Phenotypic Drug Susceptibility Testing Results

In seven Mtb isolate pairs belonging to SIT1 and SIT190, ETO resistance was associated with single nucleotide deletions at locus Rv3854c (ethA gene); however, these findings poorly matched the pDST results. Four different deletions were identified in the SIT1 isolate pairs: 110del was detected in two isolate pairs, but only one demonstrated phenotypic resistance; in an isolate pair harboring 768del, only the subsequent isolate exhibited phenotypic resistance; 1029del was detected in one phenotypically susceptible isolate pair, and for two isolate pairs harboring 1152del, pDST data were not available. 1290del was identified in the SIT190 isolate pair; however, only the subsequent isolate was phenotypically resistant and it was the only variant that could not be found in the WHO catalog. According to the pDST results, the initial SIT283 isolate was ETO-resistant, but no resistance-related variants were identified.

Para-Aminosalicylic Acid

In the case of PAS, pDST results were available for 17 pairs of Mtb isolates and 37 isolates together, while for 55 Mtb isolates from 29 pairs, pDST results were missing. The pDST and WGS-based DST results were identical for 70.6% (12/17) of Mtb isolate pairs and 86.5% (32/37) of the studied isolates (Table 9); however, PAS-resistance-conferring variants have not yet been included in the WHO Mtb variant catalog. PAS-resistance-conferring variants were detected in six isolate pairs belonging to the SIT1 and SIT190 spoligotypes, and all variants were identified at locus Rv2447c (folC gene). Two variants were detected in Mtb isolate pairs belonging to the SIT1 spoligotype: Ile43Thr was detected in one phenotypically resistant isolate pair, and Ser150Gly was identified in four isolate pairs; however, for two isolate pairs, pDST results were not available, and in two isolate pairs, either initial or subsequent isolates demonstrated phenotypic resistance to PAS. For SIT190 isolate pair resistance was caused by the Glu153Gly variant; however, only the subsequent isolate was phenotypically resistant. The initial and subsequent isolates of the SIT42 and SIT1 isolate pairs, respectively, exhibited phenotypic resistance to PAS, with no resistance-related variants identified.

Table 9 Para-Aminosalicylic Acid Resistance-Conferring Variants, and Comparison with Phenotypic Drug Susceptibility Testing Results

Bedaquiline, Clofazimine, Linezolid, D-Cycloserine, and Delamanid

No resistance-conferring variants of BDQ, CFZ, LZD, DCS, or DLM were detected. The results of BDQ, CFZ, and DLM pDST were available for only two subsequent isolates of SIT1 and SIT42 Mtb isolate pairs. The SIT1 isolate demonstrated phenotypic resistance to DLM. The SIT42 isolate pair and subsequent isolates of four Mtb isolate pairs belonging to the SIT1 (n=2), SIT42, and SIT283 spoligotypes were phenotypically tested for LZD resistance, and all isolates were LZD-susceptible. DCS pDST data were available for 18 isolate pairs belonging to SIT1 (n=9), SIT42 (n=6), SIT254, SIT190, and SIT283 spoligotypes, as well as for the initial and subsequent isolates of SIT42 and SIT1117 isolate pairs, respectively. Two subsequent isolates of the SIT1 and SIT42 isolate pairs, and the initial isolate of the SIT42 isolate pair, exhibited phenotypic resistance to DCS.

Discussion

In this study, we compared the phenotypic and WGS-predicted drug resistance profiles of Mtb isolate pairs with low genetic distances of ≤10 SNVs acquired from both single-episode and recurrent TB-experienced patients. This strategy allowed us to evaluate phenotypic and WGS-based DST performance from a novel perspective, comparing the results for nearly identical isolates rather than for sample duplicates, and following any possible drug resistance profile changes. In addition, the pDST data for the paired isolates were obtained at different time points, in the majority of cases - more than a year apart, reflecting a wider period.

However, this study has several limitations. The established sample inclusion criteria led to the foremost limitation of the study: a small number of drug-resistant Mtb isolates (n=50). In addition, because the medical records were studied retrospectively, some pDST results were unavailable for some medications (PZA, STR, AMK, ETO, and PAS), and acquiring missing pDST results for all studied Mtb isolates (n=92) would require substantial human and financial resources which the research team did not possess while conducting this study. These aspects significantly limited our choice of statistical data analysis approach. Furthermore, although studied Mtb isolate pairs were acquired within 17 years, the analyzed sample set could not be considered representative of the total drug-resistant Mtb strain population in Latvia, and thus no accurate Mtb drug resistance trends within this timeframe could be determined. Nevertheless, despite these limitations, the obtained WGS data provided insights into the variability of resistance-conferring variants for different Mtb genotypes in the studied geographical area and detected several issues concerning phenotypic and WGS-based DST methods. The results of this study also highlighted specific drugs, Mtb genotypes, and resistance-conferring variants, in which case a thorough DST result analysis is required to accurately determine Mtb isolate’s drug resistance and ensure effective treatment for each patient.

Resistance-conferring variants for nine anti-TB medications were identified, and for most of them (RIF, INH, EMB, FQ, STR, and AMK), the concordance between phenotypic and WGS-based DST was more than 90% when calculated for Mtb isolates or isolate pairs. However, we also found some inconsistencies; notably, in the majority of cases, they occurred in only one Mtb isolate from the pair. Two subsequent isolates exhibited phenotypic susceptibility to RIF and INH, although these isolate pairs harbored Asp435Tyr and C-15T variants in rpoB and fabG1 genes, respectively. Indeed, Asp435Tyr has been reported in RIF-susceptible Mtb isolates before,^39–41 moreover, it tends to cause low to moderate-level RIF resistance,^39,42 and mutant mycobacteria have demonstrated variable minimal inhibitory concentrations (MICs) on different media (0.25–2 µg/mL on Bactec MGIT 960 system, 0.5–64 µg/mL on Middlebrook 7H10 agar plate),^40,41,43,44 which justifies the discordance between WGS-based and pDST results. Similarly, the C-15T in the fabG1 gene alone is associated with low to moderate INH resistance, whereas in combination with Ser315Thr in the katG gene, it causes high-level resistance.^45,46 In the studied dataset, only three isolate pairs belonging to the SIT53 spoligotype, including one isolate with an identified disagreement between the phenotypic and WGS-based DST, harbored only the C-15T variant. All isolates were drug-resistant on liquid media (CC 0.1 µg/mL) if the test was conducted, while results acquired on LJ media were heterogeneous, mostly demonstrating either low-level resistance (CC 0.2 µg/mL) or drug susceptibility (Supplementary Table 1). The phenotypically susceptible isolate was tested for INH resistance only on the LJ medium, which provided inaccurate pDST results. In addition, in two isolate pairs belonging to SIT1, the Gly332Arg variant in the rpoC gene was detected, which is a compensatory non-synonymous variant that is mostly detected in combination with the common RIF-resistance-causing Ser450Leu variant (rpoB gene) in Mtb isolates belonging to the Beijing genotype, which enhances the transmission of RIF-resistant Mtb strains.⁴⁷

A few more mismatches between the phenotypic and WGS-based DST results were found among the EMB-resistant Mtb samples. The initial isolate of the SIT1 Mtb isolate pair harboring the Met306Val variant demonstrated phenotypic susceptibility to EMB. Such cases are frequently reported^26,28,46 because this variant has been associated with a highly variable MIC^46,48 that can be lowered to 2 µg/mL (Middlebrook 7H9 broth), which is below the CC of 5 µg/mL used in pDST; therefore, mycobacteria appear to be EMB-susceptible. This case is particularly notable because without WGS-based DST data, the pDST results indicated the development of EMB resistance in the subsequent isolate. Furthermore, we detected multiple drug-resistance-conferring variants, both present and absent in the TB-Profiler database, in phenotypically resistant Mtb isolates, whose significance in resistance development is still under evaluation according to the currently available WHO data. The Val7Gly variant in the embR gene was included in the TB-Profiler database; however, it was detected in a drug-susceptible SIT42 isolate pair. Commonly, EMB resistance is caused by genomic variants in the embCAB operon, and more rarely, it is associated with variants in the embR gene that encodes arabinosyltransferase activity level modulating transcriptional regulator.^49,50 Little is known about this drug resistance mechanism, highlighting the need for further studies to assess the association of genetic variants in the embR gene with Mtb EMB resistance.

Although PZA has been the first-line TB medication for decades, susceptibility testing remains a challenge in clinical laboratories worldwide. Phenotypic testing requires an acidic environment, which is crucial for PZA activity, and makes the drug most active at a pH of 5.5, and almost inactive when the pH is neutral.⁵¹ Furthermore, to obtain reproducible pDST results for any anti-TB medication, the inoculum preparation of strict standardized size is crucial, as because of increased bacterial cell density drug activity can be reduced leading to increased MIC.⁵² In the case of PZA, if the bacterial inoculum is too rich and non-homogenous, it can also increase the pH level and reduce PZA activity,⁵³ which makes achieving reliable and reproducible pDST results very difficult, leading to false PZA resistance interpretation and complicating the process of identifying potential resistance-conferring variants.⁵⁴

In our study, six mismatches between phenotypic and WGS-based DST were identified; however, we faced a different challenge in acquiring a precise WGS-predicted drug resistance profile, which, to our knowledge, has not been previously reported. For the three Mtb isolate pairs belonging to the SIT1 spoligotype (Beijing family), WGS data for the pncA gene could not be obtained, and this reference sequence region had zero coverage. Although deletions of the whole pncA gene have been identified as a cause of PZA resistance,^55,56 amplification of the target gene and Sanger sequencing were successful; therefore, WGS data were interpreted inaccurately using the TB-Profiler tool. Further investigation is needed to understand the reasons for unsuccessful pncA gene sequencing from fragmented mycobacterial genome samples and how to avoid this, since this genomic region is crucial for WGS-based DST. Moreover, Sanger sequencing revealed heteroresistance in two Mtb isolates from these pairs. The initial isolate of the first isolate pair had a homozygous resistant allele, and the subsequent isolate had a mixed state of PZA-resistant and PZA-sensitive alleles, whereas the heterozygous resistance-associated variant was found in the initial isolate of the second isolate pair, but the subsequent Mtb isolate did not harbor this variant. All isolates were phenotypically PZA-resistant, except for the subsequent isolate of the first isolate pair when pDST results were not available. Considering the reduced frequency of resistant alleles in subsequent Mtb isolates,⁵⁷ we hypothesized that both patients were infected with a polyclonal mycobacterial population rather than experiencing PZA resistance development during treatment, and Mtb sub-clones could have resided in varying proportions in different lung lesions, which led to sampling bias during the sputum collection procedures⁵⁸ following the acquisition of incomplete sequencing data of Mtb samples representing only the drug-susceptible portion of the Mtb population. These findings highlight the importance of molecular studies targeting active TB-causing mycobacterial populations as well as heteroresistance and its triggering factors.

The highest number of discordances between the phenotypic and WGS-based DST results were identified for ETO and PAS. Both are Group C medications for TB treatment;¹ therefore, pDST was conducted relatively less often than for first-line drugs, which is the reason for the missing pDST data for half of the genotypically resistant Mtb samples. ETO is a thermolabile compound, making it difficult to obtain reliable pDST results.¹³ Furthermore, many previously determined ETO-resistance-associated variants in fabG1 and ethA genes, including both SNVs and indels, cause a minor increase in MIC, resulting in an overlap between MICs of mutant and wild-type strains.^46,59,60 Therefore, it is uncertain whether a variant is associated with phenotypic ETO resistance. The data acquired in this study reflect this issue because, based on the WHO evaluation, the majority of detected variants were graded as being associated with ETO resistance; however, pDST results among the studied isolate pairs were highly variable, and in some cases, incorrectly indicated ETO resistance development in the subsequent isolates.

In contrast, PAS resistance-associated variants harboring Mtb strains previously demonstrated significantly higher MIC values than those of CC (1 µg/mL on MYCOTB plates), particularly if variations were detected in the folC gene.⁵⁹ Nevertheless, the variants detected in the folC gene in this study were identified in both PAS-susceptible (Ser150Gly) and PAS-resistant isolates.^25,26,61 Although this medication is widely used worldwide for drug-resistant TB treatment, the WHO does not recommend performing a pDST for PAS because no reproducible and reliable methods have been established.¹³ The CC to distinguish PAS-susceptible and PAS-resistant phenotypes correctly is yet to be determined,⁶² which could have interfered with the acquisition of an accurate Mtb pDST profile.

Unfortunately, drug-resistance-associated variant databases are still incomplete, especially for drugs that were recently introduced into drug-resistant TB treatment algorithms such as BDQ, LZD, CFZ, and DLM. Mtb resistance to these compounds is not as widespread as that to first-line drugs, and although some resistance markers have already been identified,^32,63 the proper variant validation process will take time. The TB-Profiler software used in this study demonstrated the highest sensitivity in comparison with other Mtb drug resistance prediction tools; however, it lacks specificity²³ because its database includes some variants with a low level of confidence. For instance, variants in the ethA gene related to ETO resistance, and the generated drug resistance profile should be analyzed thoroughly.⁵⁹ At the same time, our results suggest that the TB-Profiler database needs to be updated, and in some cases, the software could not identify potential resistance-conferring variants that have been included in the WHO Mtb variant catalog as variants of uncertain significance mostly identified in drug-resistant strains and were detected only in phenotypically resistant isolates in our dataset. Moreover, the lack of reliable pDST methods for PAS also applies to DCS,^64,65 significantly complicating the identification of potentially resistance-related variants. Finally, there are wild-type Mtb strains that demonstrate higher MIC values than widely used CCs,^59,60 which highlights the existence of another drug resistance development mechanism that is not necessarily associated with the mechanism of action of the drug; for instance, it could be efflux-mediated⁶⁶ or provoked by an epigenetic pathway.⁶⁷

We strongly suggest that performing pDST simultaneously with WGS-based DST will provide the most comprehensive drug resistance profile because the most important limitations of both methods can be eliminated.⁶⁸ In this study, cases of some Mtb isolate pairs opposing pDST results were obtained when the test was technically complicated, and the results were difficult to interpret accurately (PZA, ETO, PAS), or mycobacteria harbored drug-resistance-associated variants, which might lower the MICs below used CCs (RIF, INH, EMB, ETO). This implies that certain variants can affect mycobacterial MICs differently even when Mtb strains are closely related. On the other hand, pDST can provide reliable results in cases when the drug-resistance-conferring genetic variants are yet unknown or when the resistance is triggered by a different mechanism, as discussed previously. The pDST has been considered a golden standard for Mtb isolate drug resistance testing for decades, while the WGS-based DST strategy has mostly been applied for research purposes. At this point, we suggest that the WGS-based method might be introduced into routine practice in clinical laboratories as a first-choice Mtb DST strategy, and pDST could be applied only when necessary, for instance, when resistance-conferring variants to first-line drugs have been identified.

We also noticed that drug-resistance-conferring variants and their combinations varied between Mtb genotypes. In this study, drug-resistant Mtb isolate pairs belonged to Lineages 2 (SIT1 and SIT190) and 4 (SIT42, SIT53, SIT283, SIT1117), which corresponds to the previously acquired results.⁹ FQ, STR, and AMK-resistance-related variant diversity between Mtb genotypes was not evident; however, variants varied between lineages and spoligotypes for EMB and PZA, respectively. Only a combination of two identified INH resistance-conferring variants was detected in all SIT42 Mtb isolate pairs, whereas deletions in the ethA gene associated with ETO resistance were identified only in Lineage 2 isolate pairs. Furthermore, all PAS-resistant isolate pairs belonged to Lineage 2, and RIF resistance-associated variant diversity was detected among the SIT1 isolate pairs. Similar observations have been made previously, and the association between Mtb lineages and drug resistance patterns,⁶ lineage-specific drug-resistance-conferring variants,⁵ and variants causing different resistance levels in different lineages⁶⁹ has been identified. Detailed data on drug resistance variability among different Mtb genotypes are important for TB surveillance, and this needs to be investigated in larger sample cohorts.

Conclusion

Our study supports the common opinion that phenotypic and WGS-based DST should be conducted simultaneously, if possible, as this approach would minimize the risk of inaccurate Mtb drug resistance prediction. Further studies are required to better understand the distribution of drug-resistance-conferring variants among different Mtb genotypes.