Common Variants in the Glycerol Kinase Gene Reduce Tuberculosis Drug Efficacy

1The currently used multidrug chemotherapy regimen for tuberculosis (TB) was developed in a series of clinical trials in the 1980s and remains the standard of care for this disease. Infections with drug-sensitive strains of Mycobacterium tuberculosis are treated with a 6-month regimen that includes four drugs, isoniazid (INH), rifampin (RIF), pyrazinamide (PZA), and ethambutol (EMB). While this regimen cures 90% of drug-sensitive cases, the long period over which antibiotics must be administered represents a major limitation. Not only is a complete regimen difficult to deliver, but even in clinical trial settings, incomplete sterilization leads to relapse in a significant fraction of patients. The emergence of drug-resistant strains of M. tuberculosis further confounds therapy and necessitates even longer regimens with less effective drugs.

2The factors that necessitate this extended drug regimen for TB remain difficult to dissect because the in vitro efficacy of individual drugs does not predict their effect during infection. For example, PZA is critical for sterilizing an infected host, but it has very modest activity in vitro, where it may act via different mechanisms. Conversely, both INH and RIF cause relatively rapid cell death in vitro but kill bacteria much more slowly during infection. As a result, virtually all TB drug regimens kill bacteria at a lower rate during infection than they do in axenic culture. Two general mechanisms have been proposed to explain the generally drug-tolerant phenotype that is observed during infection. Growth in mammalian tissue triggers changes in mycobacterial gene expression and metabolism that can reduce drug efficacy. In addition, a number of distinct stochastically generated subpopulations have been observed, which arise either via asymmetric cell division or nonheritable regulatory events. Some of these subpopulations are relatively insensitive to antibiotics in vitro and could prolong the treatment period necessary for sterilization. While none of these mechanisms involve heritable genetic changes, many other bacteria rely on high-frequency reversible genetic variation to produce subpopulations that are tolerant to environmental insults. This process of phase variation generally relies on specific DNA sequences, such as homopolymeric regions, that are subject to frequent mutation. While phase variation has been observed in several mycobacterial species, it has not been specifically characterized in M. tuberculosis, and its potential role in determining drug efficacy is unknown.

3To understand the processes that determine drug efficacy during infection, we employed two complementary approaches. A forward genetic study identified bacterial functions that alter drug efficacy in mice. In parallel, whole-genome sequence analysis of M. tuberculosis clinical isolates identified genetic variants in candidate genes that are associated with resistance. Together, these approaches defined a variable homopolymeric region in the glpK gene that controls glycerol metabolism and drug efficacy. Heritable genetic variation at this site produces a drug-tolerant phenotype that reduces treatment efficacy and is associated with the emergence of resistant clones.

4To specifically define bacterial functions that limit efficacy during infection, we used transposon sequencing (TNseq) to identify mutations that alter bacterial killing. Groups of mice were infected with a nearly saturated library of M. tuberculosis transposon mutants via the intravenous route. After allowing 2 weeks for bacterial growth and the establishment of adaptive immunity, animals were treated with a regimen based on first-line TB chemotherapy, a mixture of INH, RIF, PZA, and EMB (HRZE). Each drug in this regimen was shown to be effective individually at the dose given, and the four-drug mixture reduced organ burden by more than 100-fold after 14 days of treatment (Fig. 1A). To identify mutations with relatively rapid effects on bacterial killing by antibiotics, mutant pools were recovered from the spleen after 1 week of therapy by plating organ homogenates. This analysis was performed in the spleen because this organ contained an adequate bacterial population size to ensure that complexity of the library was maintained throughout the infection. In vivo-selected libraries were compared to each other using TNseq, which quantifies the relative abundance of each mutant in a given pool by sequencing all the transposon-chromosome junctions that are present.

5Mutant pools were collected from three groups of animals. A pretreatment pool was collected immediately before drug administration. One week later, pools from antibiotic-treated or untreated groups were collected. This study design allowed the relative fitness of each bacterial mutant to be assessed in the presence and absence of drug therapy. Pairwise analyses of the treated and untreated pools with the pretreatment library identified distinct sets of genes that altered bacterial representation under each condition (see Table S1 in the supplemental material). As many antibiotics act in a growth rate-dependent manner, we investigated whether bacterial fitness in vivo was an important determinant of antibiotic efficacy. When the relative abundance of each mutant in untreated animals was compared with their abundance in time-matched drug-treated mice, we found no global correlation between bacterial fitness in the presence and absence of antibiotic (Fig. 1B). However, several individual mutations were observed that reduced fitness under both conditions. Thus, to more formally focus our study on drug-related phenotypes, we performed a three-way analysis to identify those mutations that alter fitness preferentially in antibiotic-treated animals. This analysis defined 61 mutants that increased, and 8 that reduced, the effect of therapy (Table S2). These mutants corresponded to a variety of functional pathways (Fig. 1C). In several cases mutation resulted in dramatic alterations in fitness after drug exposure, as the relative representation of these mutants under the pre- and posttreatment conditions varied by more than 100-fold.

6A number of functions found to alter bacterial fitness in drug-treated mice were already known to impact drug efficacy. For example, mutants lacking PncA, which converts PZA into its active pyrazinoic acid form, were highly overrepresented in the treated mice, highlighting the singular importance of PZA in the activity of this regimen (Fig. 1D). In addition, we found that mutations in mmaA1, mmaA2, and cmaA2 sensitized the bacterium to drug treatment. All of these genes encode functions necessary for mycolate modification, and chemical inhibition of these partially redundant activities has been shown to increase cellular permeability to antibiotics. Additional protein families associated with cell wall structure altered drug efficacy. PE/PPE family members have been implicated in cell envelope integrity, and we found that mutations in the ppe50-ppe51 pair increased killing (Fig. 1E). Conversely, the loss of enzymes (ppsA, ppsC, and drrA) necessary for the synthesis of the major cell envelope lipid, phthiocerol dimycocerosate, decreased clearance. Drug access appeared to be similarly limited by multiple classes of efflux pumps, as mutations in members of the ATP-binding cassette (rv1747), major facilitator superfamily (rv3728), and MmpL (mmpL8 and mmpL10) families were found to increase bacterial clearance (Table S2).

7In addition to these known mechanisms, we identified a number of novel functions that altered bacterial killing. Prominent among these were pathways involved in carbon metabolism. For example, mutation of both assayable steps of the glycogen synthetic pathway (glgA and glgC) increased antibiotic activity (Fig. 1F). This pathway promotes carbon storage through carbohydrate anabolism, a general process that has been previously implicated in drug tolerance. In addition, we identified the glpK gene, which encodes the sn-glycerol-3 kinase of M. tuberculosis. Twenty-six of the 29 insertional mutants in this gene showed a similar decrease in clearance rate upon drug treatment (Fig. 1G). When pools isolated before and after treatment were directly compared, only pncA mutations produced a statistically significant reduction that was greater than those in glpK (Fig. 1H). This relatively dramatic phenotype was explored in more detail.

8GlpK is responsible for phosphorylating the 3-position of glycerol, which is necessary for its catabolism via the lower glycolytic pathway. To determine if glycerol catabolism per se was capable of enhancing the activity of TB drugs, the effects of INH, RIF, and moxifloxacin (MOX) were compared in media containing glycerol or other carbon sources known to be used during infection, fatty acid and cholesterol. These carbon sources supported different growth rates, which can confound endpoint-based determinations of antibiotic activity, such as standard MIC measurements. As a result, we quantified the growth rate (GR) of bacteria over a time course and determined the concentration of each drug that was necessary to decrease this rate by 50%, which is expressed as GR50. Using this approach, we found that glycerol catabolism produced a modest but reproducible decrease in the GR50 for several drugs. Growth in glycerol significantly increased the efficacy of RIF and MOX compared to growth in valerate and enhanced the efficacy of INH and MOX compared to growth in cholesterol (Fig. 2A).

9To further investigate the role of glycerol metabolism in drug efficacy, a glpK deletion mutant of M. tuberculosis was constructed. The ΔglpK mutant was unable to grow in media containing glycerol as the sole carbon source (Fig. 2B), indicating that the deleted gene encodes the sole glycerol-3 kinase activity. The effect of antibiotics on the growth rates of glpK-sufficient and glpK-deficient strains was then compared in media containing different carbon sources. When glycerol was present in the medium, the ΔglpK mutant was significantly less sensitive to INH and RIF than the wild type or the complemented mutant (Fig. 2C). This difference largely disappeared when glycerol was replaced with either the nonglycolytic substrate, butyrate, or a glycolytic product that bypasses the triose phosphate pool, pyruvate. The differential effects of these carbon sources indicated that the assimilation of exogenous glycerol was primarily responsible for glpK’s influence on drug sensitivity. PZA sensitivity was assessed at pH 5.8 to maximize the in vitro efficacy of the drug. However, under these conditions, glpK deletion did not alter PZA sensitivity. As INH, RIF, and MOX have distinct mechanisms of action, the effect of glycerol catabolism on antibiotic activity in vitro did not appear to be specific to a particular drug or target pathway.

10The mouse model was then used to explore the role of glpK during infection. Consistent with both our TNseq data and previous work, deletion of the glpK gene did not affect the growth or persistence of M. tuberculosis in the lungs of mice after aerosol infection (Fig. 3A). To quantify the effect of glpK deletion on the efficacy of individual drugs, mice were inoculated via the intravenous route with a mixture of wild-type and ΔglpK bacteria and treated with antibiotics, as was done for the initial TNseq screen. Another mutant lacking the ppe51 gene, which TNseq predicted to be hypersensitive to multidrug treatment (Fig. 1E), was included as an additional control. Since this study did not require maintaining the complex mutant mixture needed in the TNseq study, more prolonged treatment regimens could be used. As we observed previously, all drug regimens reduced the bacterial burden, and PZA or combination therapy had the greatest effect (Fig. 3B). Surviving bacteria were recovered by plating at the indicated time points, and the relative abundance of the three M. tuberculosis strains was determined by quantitative PCR. Both mutants demonstrated the predicted phenotypes in animals treated with the four-drug combination therapy for 2 weeks, as the ΔglpK mutant was cleared significantly more slowly and the Δppe51 mutant significantly more rapidly than the wild type (Fig. 3C). These phenotypes were even more pronounced in animals treated with PZA alone. In contrast to the broadly sensitizing effect of glycerol catabolism in vitro, the ΔglpK mutant behaved similarly to the wild type in mice treated with RIF, INH, or EMB. The decreased efficacy of PZA against the ΔglpK mutant was also found in the lungs of mice infected via aerosol. In this model, treatment with PZA between 21 and 35 days postinfection reduced the bacterial burden of wild-type and ΔglpK complemented strains by at least 1,000-fold but had a significantly reduced effect on the glpK-deficient mutant (Fig. 3D). Thus, while glycerol catabolism can nonspecifically alter antibiotic susceptibility in vitro, glpK deletion preferentially reduced the effect of PZA and a PZA-containing multidrug regimen in this animal model.

11As the glpK deletion did not alter bacterial fitness during infection and conferred a benefit upon drug treatment, we hypothesized that mutations altering glycerol catabolism are positively selected during the evolution of drug resistance in natural populations. As a first test of this hypothesis, we characterized a panel of Korean M. tuberculosis isolates that varied in drug sensitivity profiles, from fully sensitive strains to extensively evolved clones that were phenotypically resistant to more than ten different antibiotics (see Table S3 for strain characteristics). To investigate whether glycerol catabolic defects were selected during the evolution of resistance in these strains, we subcultured a random subset of drug-sensitive or extensively resistant isolates in media containing glycerol as the sole carbon source. Drug-sensitive strains grew at a rate similar to that of a standard laboratory strain (H37Rv). However, while the extensively drug-resistant isolates could grow in butyrate, none of the tested isolates could grow in the glycerol-containing media (Fig. 4A and Table S3).

12The whole-genome sequences (WGS) of these isolates were determined. Based on WGS, this collection was predominantly comprised of a Korean sublineage of East Asian strains, and the multidrug-resistant (MDR) phenotypes could generally be attributed to known high-level resistance-conferring mutations (Table S4). Inspection of the WGS data revealed that the glycerol catabolic defect in 9 of the 11 tested strains could be attributed to loss-of-function mutations in the glpK gene. These strains all harbored a one-base expansion of the same homopolymeric sequence (GGGGGGG) in the 5′ half of the glpK open reading frame. The sequence of the glpK homopolymer was verified in the entire panel by targeted Sanger sequencing (Fig. 4B). This mutation is predicted to eliminate GlpK enzymatic activity, as it introduces a premature termination codon that eliminates the majority of the open reading frame, and the same homopolymer expansion has been previously observed in M. bovis strains lacking glycerol kinase activity. An additional missense mutation altering amino acid 169 was identified in an otherwise glpK wild-type allele, but the functional significance of this mutation is unclear. Two of the phenotypically glycerol-deficient strains carried no obvious mutations in glpK or other glycerol catabolic genes. Thus, while glpK frameshifts appear to be the most common lesion associated with glycerol catabolic defects in this collection, other mechanisms contribute in a fraction of isolates.

13All ten of the glpK frameshifts identified in this panel were found in multidrug-resistant strains, particularly the highly evolved strains that were resistant to more than eight different drugs (Fig. 4C). Based on the phylogenetic relationship between these strains, the identified glpK frameshifts represent at least three independent mutational events. Alternative tree topologies necessary to accommodate fewer mutational events were significantly less likely (P < 10−4). Thus, in this relatively small collection of strains, inactivating mutations in glpK were frequent and associated with defective glycerol utilization and drug resistance.

14Frameshift mutations in homopolymeric DNA sequences can represent high-frequency and reversible events. To assess the frequency of this mutation in a larger population and to further explore its association with drug resistance, we analyzed the whole-genome sequences of a larger collection of isolates from Peru. Of 1,031 sequenced strains, 68 isolates harbored nonsynonymous variants in the glpK gene. Of these, 45 contained a single-base expansion of the glpK homopolymer and 2 contained a two-base expansion. These frameshifts were found in all phylogenetic clades of M. tuberculosis, indicating that this mutation has arisen through multiple distinct mutational events in different lineages of the pathogen (Fig. 5A and B). In total, homopolymer expansion accounted for 66% of the nonsynonymous glpK variants, and 4.6% of all isolates harbored frameshift mutations disrupting the glpK open reading frame.

15As we found in the smaller set of Korean strains, glpK frameshifts were significantly associated with drug resistance. No instances of glpK frameshifts were found in the 90 phenotypically drug-sensitive strains. In contrast, 44 of the 739 isolates that met the WHO criteria for MDR carried these mutations (P = 3 × 10−5). In this large collection, we were also able to test the association between glpK genotype and resistance to individual drugs and found significant associations with RIF, INH, and ethionamide (ETH) (Fig. 4C). The lack of observed association with PZA resistance could have been due to the unreliability of this phenotypic assay. Indeed, using a genotypic assay we found a significant association between glpK frameshifts and nonsynonymous pncA variants (P = 0.001). Thus, the glpK homopolymer is the site of the majority of variation in this gene, and expansion of this hypervariable region is associated with the evolution of drug resistance.

16Drug tolerance has been proposed to contribute to both relapsing TB disease and the emergence of drug-resistant clones. Most current models to explain tolerance in mycobacteria are largely restricted either to nonheritable processes, such as changes in gene expression or functionally asymmetric cell division, or stably heritable mutations. Our data provide a new mechanism by which alterations in a hypervariable region in the glpK gene produces bacteria that persist during antibiotic treatment and could contribute to the emergence of drug-resistant clones.

17Insertions and deletions in a homopolymeric region of an open reading frame is a common mechanism to produce high-frequency reversible phenotypic variation in bacteria. These mutations are generally thought to result from slipped-strand mispairing during DNA replication. However, additional DNA repair mechanisms, such as mismatch repair or base excision repair, can alter the frequency and directional bias of the process. The exclusive bias for +1 and +2 frameshifts in the glpK gene of clinical isolates argues for a more complex process than simple replicative error, and the ultimate frequency of mutants may also be influenced by the specific fitness effect of each frameshift. As a result, it is difficult to anticipate the rate of variation that occurs during infection. Regardless, we identified +1 and +2 frameshifts in 6% of the Peruvian MDR isolates. Since clinical samples are routinely cultured in glycerol-containing media, which would be expected to select for reversion to the wild-type glpK coding sequence, it is likely our data underestimate the frequency of glpK mutants. The observed prevalence of this mutation is clearly high enough to produce a significant population of glpK-deficient clones that alter drug efficacy.

18Our studies in the mouse model demonstrate that glpK-deficient bacteria are drug tolerant during infection. The mechanism(s) that underlies the drug tolerance of glpK-deficient bacteria is likely to be complex. Common fates for glycerol-3-phosphate are catabolism via the lower glycolytic pathway, incorporation into anabolic pathways, and spontaneous degradation to methylglyoxal. As a result, glycerol assimilation can alter growth rate, metabolism, and cellular structure. While not yet conclusive, our in vitro studies argue against some of these mechanisms. Both the in vitro effects of glycerol supplementation and the in vivo effects of glpK expression were independent of growth rate (Fig. 1 and 2), and the differential effect of glpK in glycerol versus pyruvate growth media indicates that glycolytic flux per se is not the major determinant of drug efficacy. Thus, we speculate that the abundance of the triose phosphate pool or some derivative of this pool is primarily responsible for the general enhancement of antibiotic efficacy that we observed upon glycerol assimilation in vitro.

19The nonspecific effect of glpK on multiple drugs that we observed in vitro is consistent with previous studies that identified frameshifts in the glpK homopolymer in mutants selected for spontaneous resistance to investigational antimycobacterial compounds in vitro. Similarly, glpK mutations have also been found in conjunction with additional mutations in strains selected to be drug tolerant or PZA resistant in vitro. Despite these relatively general effects on drug activity in vitro, glpK deletion preferentially reduced the efficacy of PZA-containing regimens in the mouse model used in this study. This apparent discrepancy could reflect differences in drug exposure, bacterial physiology, or GlpK functions in these two settings. The poor activity of PZA in vitro, where we observed no effect of glpK on PZA activity, makes it difficult to dissect these mechanisms in a more controlled system. Regardless, the identification of mutations that affect PZA efficacy only during infection highlights the importance of performing the original TNseq screen in an animal model.

20The prevalence of glpK-deficient strains in natural populations and the preferential survival of these bacteria in drug-treated animals suggested that glpK-deficient clones contribute to the persistence of M. tuberculosis during therapy and provide precursors for the emergence of clones with high-level resistance-conferring mutations. It is unlikely that the effects of glpK variation would be noted in standard phenotypic drug susceptibility testing (DST). This situation is similar to common variation in the prpR gene, which specifically influences drug tolerance but not DST results. These observations raise the possibility that genotypic tests for common drug tolerance-inducing variants could predict treatment failure and eventually be used to tailor therapy. We note that our TNseq study identified a number of additional loss-of-function mutations that alter drug efficacy, and the genome contains more than 100 genes with homopolymeric regions that are at least as long as the one found in glpK. Together, these observations suggest that many phenotypically distinct subpopulations arise via similar mechanisms and influence antibiotic efficacy.

21BALB/cJ (stock no. 000651) mice were purchased from the Jackson Laboratory (Bar Harbor, ME, USA). Housing and experimentation were in accordance with the guidelines set forth by the Department of Animal Medicine of University of Massachusetts Medical School and Institutional Animal Care and Use Committee and adhered to the laws of the United States and regulations of the Department of Agriculture. Eight- to 12-week-old female animals were infected with 106 CFU of a himar1 transposon library via the intravenous route. Groups of mice were treated with antibiotics starting at 14 days postinfection. Antibiotics were administered via drinking water at the following concentrations: 0.1 g/liter isoniazid (Sigma), 0.6 g/liter ethambutol (Sigma), 0.1 g/liter rifampin (Sigma), and 15 g/liter pyrazinamide (Sigma). At the indicated time points, mice were sacrificed, spleens and lungs were isolated and homogenized, and CFU numbers were determined by plating dilutions on 7H10 agar with 10 μg/ml kanamycin. For library recovery, approximately one million CFU per mouse were plated on 7H10 agar with kanamycin (10 μg/ml). Genomic DNA was extracted, and the relative abundance of each mutant was estimated as described previously. Statistical analysis of log2 fold change (log2FC) in counts between conditions (two-way analysis) was performed by resampling. The three-way analysis measures the difference in log2FC (Δlog2FC) measured under two selective conditions relative to a common starting condition: Δlog2FC = log2FC(condition 1) − log2FC(condition 2).

22In the present case, condition 1 was 14 days postinfection plus 7 days of antibiotic treatment, condition 2 was 21 days postinfection, and the starting condition was 14 days postinfection (the start of drug treatment). Statistical significance was assessed by resampling. For each gene, the sampling distribution of Δlog2FC was obtained by resampling with replacement of the insertion counts at each TA within the gene (after normalization across all libraries). Counts for replicates were pooled prior to resampling. For each of 10,000 resamples, Δlog2FC was calculated. The P value was taken as the fraction of the cumulative frequency distribution of Δlog2FC falling outside Δlog2FC = 0, on the negative side for values measured as Δlog2FC > 0, or on the positive side for values measured as Δlog2FC < 0 (equivalent to a 1-tailed test). The resulting P values were adjusted for multiple testing by Benjamini-Hochberg false discovery rate.

23M. tuberculosis H37Rv was maintained in Middlebrook 7H9 medium containing oleic acid-albumin-dextrose-catalase (OADC), 0.2% glycerol, and 0.05% Tween 80 and grown with shaking (200 rpm) at 37°C. Hygromycin (50 μg/ml) or kanamycin (20 μg/ml) was added when necessary. All work with M. tuberculosis adhered to the CDC-NIH Guide for Biosafety in Microbiological and Biomedical Laboratories. glpK and ppe51 were deleted by allelic exchange as described previously, and this work adhered to NIH guidelines for research involving recombinant DNA molecules. Nucleotides 4138237 to 4139720 (glpK) or 3501829 to 3502901 (ppe51) were replaced by the vector pKM464, including qTag-22 or qTag-26 for glpK and ppe51, respectively. The ΔglpK strain was cultured in glycerol-free 7H9. Glycerol-dependent growth was assessed in minimal medium containing asparagine (0.5 g/liter), KH2PO4 (1 g/liter), Na2HPO4 (2.5 g/liter), ferric ammonium citrate (50 mg/liter), MgSO4·7 H2O (0.5 g/liter), CaCl2 (0.5mg/liter), ZnSO4 (0.1mg/liter), 0.1% tyloxapol, and either 0.1% glycerol or 0.1% dextrose. For in vitro antibiotic susceptibility testing, isoniazid and rifampin were used at 2 and 1 μg/ml, respectively, and serially diluted 2-fold. Bacteria were inoculated to a starting optical density at 600 nm (OD600) of 0.05 in 96-well plates with 7H9 medium containing OADC, 0.05% Tween 80, and 0.2% glycerol, butyrate, or pyruvate. Pyrazinamide was used at 400 μg/ml and serially diluted 2-fold. Bacteria were inoculated to a starting OD600 of 0.01 in inkwells containing 7H9 medium supplemented with OADC, 0.2% glycerol, and 0.05% tyloxapol at pH 5.8 and grown with shaking. Growth was monitored by OD600. Conditions were assessed in triplicate. Antibiotic efficacy was determined by comparing growth rate under increasing drug concentrations. OD600 was plotted and the rate constant (k) value was determined for all conditions using an exponential growth model. Rate constants posttreatment were normalized to levels for no-antibiotic controls.

24ΔglpK mutant fitness in vivo was determined by inoculating mice with a 1:1 mixture of ΔglpK (hygromycin resistant) and H37Rv (harboring pJEB402 chromosomally integrated plasmid encoding kanamycin resistance) strains via the aerosol route. At the indicated time points, mice were sacrificed and CFU numbers in spleen and lung homogenate were determined by plating on 7H10 agar. Fitness in the presence of antibiotic was assessed by a similar competitive assay. Mice were infected with a pool of strains at equal ratios via the intravenous route (106 total CFU/mouse). Groups of mice were treated with antibiotics starting at 14 days postinfection, as described for the TNseq study. At the indicated time points, approximately 10,000 CFU from the spleen homogenate of each mouse were plated on 7H10 agar. Genomic DNA was extracted for quantitative real-time PCR analysis. Briefly, the abundance of the constant and variable regions of the q-Tag present in each mutant was determined by TaqMan PCR assay, as described previously, and used to calculate a variable/constant region ratio for each strain. The abundance of each mutant strain was then plotted relative to that of wild-type H37Rv (mutant/wild type). Values were normalized to initial day 0 ratios.

25Bacteria were grown in minimal medium with 0.1% glycerol, 0.1% valeric acid, or 0.1% cholesterol on 96-well plates. Isoniazid, rifampin, and moxifloxacin were used at 1, 0.062, and 1 μg/ml, respectively, and serially diluted 2-fold. A no-antibiotic control was included in each experiment. Bacteria were inoculated to a starting OD600 of 0.05, and growth was monitored by OD600 and fluorescence. Conditions were prepared in triplicates. Antibiotic efficacy was determined by growth rate inhibition. The exponential growth constant (k) value was determined for all conditions. The k value of each antibiotic concentration was normalized to the k value of the no-drug control. The GR50 value was determined as the concentration of antibiotic that resulted in a 50% decrease in growth rate, as previously described.

26Strains were collected from the National Culture Collection for Pathogens, which is maintained by the Korea Centers for Disease Control and Prevention. Phenotypic DST testing for each strain was conducted by an absolute concentration method using Löwenstein-Jensen agar with critical concentrations of TB drugs (in μg/ml): isoniazid (0.2), rifampin (40), ethambutol (2), streptomycin (10), kanamycin (40), prothionamide (40), cycloserine (30), para-aminosalicylic acid (1), ofloxacin (2), pyrazinamide (50; pH 4.65), capreomycin (40), moxifloxacin (2), amikacin (40), levofloxacin (2), p-nitrobenzoic acid (500), and rifabutin (40). Strains were classified as resistant if drug-containing media produced more than 1% of the CFU observed in control cultures. To test growth on glycerol, M. tuberculosis was grown in Middlebrook 7H9 broth with 0.5% glycerol, 0.05% tyloxapol, catalase, and fatty acid-free bovine serum albumin (Sigma). Inocula were cultivated in 7H9-OADC-Tween 80 to an OD of 0.1 to 0.2, washed with phosphate-buffered saline plus tyloxapol (0.05%), and diluted to ∼106 CFU/ml.

27Genomic DNA was sequenced either by Ion Torrent (yielding an average read length of 170 bases) or Illumina (300-base paired-end reads) platforms. In both cases, reads were aligned using bwa mem (version 0.7.12) against Mycobacterium tuberculosis H37Rv reference GenBank accession no. NC_018143.2. Variants were called using GATK 3.3-0 by following the developer’s best practices: (i) picard 1.96 MarkDuplicates, (ii) GATK Realigner Target Creator, (iii) GATK IndelRealigner, (iv) GATK BaseRecalibrator, (v) GATK UnifiedGenotyper, and (vi) GATK GenotypeGVCFs. Base recalibration was performed iteratively using the initial Ion variant calls, obtained without recalibration, to obtain a set of polymorphic sites for use in step 4. Final filtering was performed separately for single-nucleotide polymorphism (SNP) and insertion/deletion (indel) calls: for SNPs, FS of >60.0, MQ of <40.0, MQRankSum of <−12.5, and ReadPosRankSum of <−8.0; for indels, FS of >200.0 and ReadPosRankSum of <−20.0. Call-passing filters were combined and a final filter QD of <20.0 was applied. Only calls passing all filters were combined into the final vcf file (a total of 7,418 variants). For phylogenetic analysis, only variant SNPs were used. The final alignment (which included the corresponding NC_181843.2 reference bases) consisted of 51 taxa × 7,247 positions. A maximum likelihood phylogeny was obtained using PHYML version 20120412 with the generalized time-reversible model. Trees were visualized with iTOL. Alternative tree topologies were generated using PHYLIP retree, and SH tests were performed using PAML baseml, version 4.8.

28Drug-resistant phenotypes were determined by measuring MICs to 12 antituberculosis drugs. For strains found to be sensitive at the critical concentration recommended by the WHO for each drug, we tested two MIC levels below the critical concentration, and for those resistant at the critical concentration, we tested six levels above it. The testing concentrations deviated from the traditional doubling in order to better detect intermediate-level MICs that are theoretically achievable levels in patient sera.

29Strains were sequenced on the Illumina platform to produce 100 to 150 paired-end reads and coverage of at least 50-fold. The paired-end raw sequence data were mapped to the H37Rv reference genome using the BWA mem algorithm. We used SAMtools (default settings) and pilon to identify single-nucleotide variants and insertions and deletions up to approximately 100 bp using a coverage-based approach. We assigned a variant call as missing if the valid depth of coverage at a specific site was less than 10 reads, if the mean read-mapping quality at the site did not reach 7, or if none of the alternative alleles accounted for at least 90% of the valid coverage. The genotype of glpK was determined by a binary burden score that represented the presence of any nonsynonymous single-nucleotide variant, insertion, or deletion observed in that gene. M. tuberculosis genetic lineages were determined using a previously published SNP barcode. A neighbor-joining tree was derived using https://cran.r-project.org/web/packages/ape/index.html.

30A linear mixed model was used to examine the associations between glpK genotype and the rank-transformed MICs phenotypes of 12 drugs and pncA genotype. We adjusted for the population structure using a genetic relatedness matrix (GRM), calculated from a pairwise distance matrix using synonymous single-nucleotide variants of the complete genome using the software GEMMA. The type I error rate was set at 0.01 after a Bonferroni correction accounting for the multiple comparisons. The linear mixed model was also performed using GEMMA.