Genomic and functional analyses of Mycobacterium tuberculosis strains implicate ald in D-cycloserine resistance.

A more complete understanding of the genetic basis of drug resistance in Mycobacterium tuberculosis is critical for prompt diagnosis and optimal treatment, particularly for toxic second-line drugs such as D-cycloserine. Here we used the whole-genome sequences from 498 strains of M. tuberculosis to identify new resistance-conferring genotypes. By combining association and correlated evolution tests with strategies for amplifying signal from rare variants, we found that loss-of-function mutations in ald (Rv2780), encoding L-alanine dehydrogenase, were associated with unexplained drug resistance. Convergent evolution of this loss of function was observed exclusively among multidrug-resistant strains. Drug susceptibility testing established that ald loss of function conferred resistance to D-cycloserine, and susceptibility to the drug was partially restored by complementation of ald. Clinical strains with mutations in ald and alr exhibited increased resistance to D-cycloserine when cultured in vitro. Incorporation of D-cycloserine resistance in novel molecular diagnostics could allow for targeted use of this toxic drug among patients with susceptible infections.

Introduction

Molecular diagnostic testing has transformed pathogen identification and subsequent drug susceptibility testing (DST) in the field of clinical microbiology[1]. In particular, the recent adaptation of rapid genotypic diagnostics for Mycobacterium tuberculosis[2,3] has revolutionized the detection of M. tuberculosis and determination of drug resistance by assaying for the presence of specific resistance-conferring mutations. To minimize morbidity and mortality of drug-resistant M. tuberculosis, it is critical to treat patients with antitubercular drugs to which their infecting strains are susceptible[4,5], and also to prescribe the minimally toxic regimen that is capable of curing the disease[6]. However, current genotypic diagnostics for M. tuberculosis are limited in that they identify only a subset of previously known resistance mutations, and the catalog of mutations that confer resistance is incomplete. For example, the Hain Genotype MTBDRsl test[3] (see URLs) can provide insight into second-line drug resistance in M. tuberculosis strains, but only interrogates a limited number of resistance mutations and fails to diagnose one in three to four cases of extensively drug-resistant tuberculosis (XDR TB)[7]. As many second-line drugs included in MDR regimens are highly toxic[8-12] incomplete data provided by current molecular tests may result in needless exposure to toxic agents, for example in the case of rifampicin mono-resistance detected up by GeneXpert[2] that is treated as MDR.

D-cycloserine is one such highly toxic second-line drug that is not evaluated by current molecular diagnostics. The World Health Organization (WHO) has classified D-cycloserine as a category 4 drug for the treatment of TB, indicating that it is an oral bacteriostatic second-line agent. D-cycloserine remains a cornerstone of MDR and XDR-TB treatment despite severe psychiatric and central nervous system toxicities[9,10,12]. With these associated toxicities, it would be ideal to further restrict its use to patients harboring strains that are susceptible to D-cycloserine. Phenotypic DST for D-cycloserine has only been validated on Löwenstein-Jensen media; however, routine phenotypic DST is performed infrequently and is not recommended by the WHO due to technical difficulties with the assay[13-15].

Given these difficulties with phenotypic resistance characterization and toxicities, it would be ideal to detect D-cycloserine resistance via rapid molecular diagnostics. However, the molecular basis of D-cycloserine resistance among clinical isolates of M. tuberculosis is poorly understood. D-cycloserine functions by competitively inhibiting two enzymes in the alanine metabolism pathway—alanine racemase (encoded by alr, Rv3423c)[16], and D-alanine D-alanine ligase (encoded by ddlA, Rv2981c)[17,18]—which are essential for peptidoglycan synthesis. Overexpression of alr and ddlA have been previously shown to confer resistance to D-cycloserine in Mycobacterium smegmatis[19,20], and acquisition of a nonsynonymous mutation in alr in M. tuberculosis was shown to be correlated with the acquisition of phenotypic resistance to D-cycloserine[21,22]. In addition, Mycobacterium bovis BCG, a live attenuated vaccine strain derived from M. bovis, is intrinsically resistant to cycloserine[23], a phenotype that is only partially explained by a mutation in cycA (an alanine transporter)[24]. However to our knowledge no mutations in M. tuberculosis have been experimentally validated as conferring resistance to D-cycloserine and there is uncertainty over the critical concentration that defines resistance.

Recently, whole-genome sequencing has provided more complete information about the genetic content of drug-resistant M. tuberculosis[21,25-28], including investigations that have used signatures of convergent evolution and positive selection to search for novel genes involved in drug resistance[25,26]. Here, we combined two large studies of M. tuberculosis[26,28] with comparable whole-genome sequencing and phenotypic drug resistance data to identify mutations that associate with drug resistance. These two studies, from South Africa and China, included a total of 498 strains with varying degrees of drug resistance, including large numbers of MDR and XDR strains. By utilizing strategies for amplification of signal from rare variants, and by combining an association test with a correlated evolution test, we identified the gene L-alanine dehydrogenase, ald (Rv2780), as a candidate novel drug resistance-conferring gene. We then demonstrated that loss-of-function of ald confers resistance to D-cycloserine in vitro, and that clinical strains with ald loss-of-function mutations were resistant to this drug.

Results

Identification of candidate resistance-conferring variants

We first compiled a dataset of 498 strains of M. tuberculosis from South Africa and China, including whole-genome sequencing and phenotypic drug resistance data[26,28]. We next sought to assess how much unexplained drug resistance was present among these strains by comparing genomic variants detected in each strain to a list of resistance-associated variants for eight drugs drawn from current diagnostics and literature (Supplementary Note, and Supplementary Tables 1–2). The amount of unexplained drug resistance varied by drug, with anywhere between 8 and 21 strains without known resistance conferring mutations per drug (Table 1). Given the limited numbers of strains that were phenotypically resistant but lacked known genotypic resistance mechanisms to individual drugs, we reasoned that novel mechanisms of resistance might be rare and difficult to identify statistically based on individual variants. We therefore collapsed the genotypic data using two approaches: first by frequency in the population, and second based on functional impact (Methods). Specifically, in the second approach, we collapsed different loss-of-function mutations (defined as a nonsense mutation, frameshift, or large in-frame insertion or deletion (indel)) within the same gene, as loss-of-function of select genes are known to confer resistance to select drugs in M. tuberculosis[29-33]. Synonymous single-nucleotide polymorphisms (SNPs) were excluded in both approaches. Using these criteria, a combined set of 11,704 genotypic features was generated for downstream analysis.

Table 1

Unexplained drug resistance in 498 strains of M. tuberculosis from South Africa and China. The total number of strains assayed for phenotypic resistance to each drug is shown, along with the number of phenotypically resistant (R) and phenotypically susceptible (S) strains. The number of strains with phenotypic drug resistance but without a match to known resistance-associated mutations is also given.

	# of phenotyped strains			# strains w/unexplainedR
Drug	Total	R	S
Rifampicin	498	337	161	8
Isoniazid	498	350	148	21
Pyrazinamide	142	72	70	12
Ethambutol	394	164	230	15
Streptomycin	498	256	242	21
Ofloxacin	494	162	332	15
Kanamycin	491	104	387	20
Capreomycin	300	66	234	15
Ethionamide	288	102	186	16

To identify novel resistance conferring mutations in this dataset, we hypothesized that resistance-conferring genotypes should follow two patterns: 1) there should be a positive correlation between presence of a resistance-conferring mutation and presence of phenotypic resistance, and 2) a resistance-conferring mutation and phenotypic resistance should evolve convergently and in concert, i.e., they should evolve repeatedly and at the same times in evolutionary history. We therefore conducted two tests: a Fisher’s exact test to measure the direct association between genotype and phenotype, and a test for correlated evolution[34] of genotype and phenotype, and then combined the p-values from both tests into an F-score (Methods).

Analysis of this combined set of genotypic features with the Fisher’s exact and correlated evolution tests for each drug revealed that, in almost all cases, mutations already known to confer drug resistance were the most significant genotypic features resulting in the highest combined F-scores (Figure 1, Supplementary Table 3). However, for most drugs, genotypic features known to confer resistance to drugs other than the one being analyzed in the association tests were also often highly significant across both tests yielding high F-scores. For example, mutations in katG, pncA, rpoB, embB, and rpsL were significantly associated with isoniazid resistance despite the fact that katG is the only gene on this list which confers resistance to isoniazid. This is due to the fact that drugs to treat M. tuberculosis are prescribed only in combination and, therefore, resistances frequently co-occur.

Figure 1

Significance of associations between genotypic variants and drug resistance phenotypes. Each circle represents a genotypic feature that is plotted at the intersection of the negative log p-values from Fisher’s exact and correlated evolution tests for each drug. P-values were corrected for multiple comparisons using the Benjamini-Hochberg method[59]. Genotypic variants known to confer resistance are colored according to the drug to which they confer resistance, and genotypes with no known effect on drug resistance are grey. Genotypes scoring well in both tests appear in the upper right quadrants.

Remaining high scoring features included genes encoding PE-PGRS proteins, polyketide synthases, secretion-associated proteins, and hypothetical proteins (Supplementary Table 3). However, when further examined, none of these other high scoring features could account for a large fraction of phenotypic resistance in strains without known genotypic resistance mechanisms. Rather, these features most often appeared in strains where phenotypic resistance could already be explained by known resistance-conferring mutations.

Loss-of-function correlated with unexplained resistance

Given that the list of highest scoring features contained a significant number of genes already associated with resistance, and that the remaining high scoring features co-occurred with known resistance-conferring mutations, we hypothesized that the strong statistical signal contributed by mutations in these genes could make the identification of novel, low frequency resistance mechanisms difficult. To amplify signal from strains with potential for novel resistance genotypes, we repeated the association analysis after removing strains with known resistance-conferring mutations (Supplementary Figure 1, Supplementary Table 4).

From this analysis, we identified 21 genes with loss-of-function mechanisms that associated with phenotypic resistance to at least one of nine drugs (excluding PE/PPE genes and transposable elements; Supplementary Tables 4, 5). Eight of these drug resistance-associated genes were represented in a library of M. tuberculosis single gene knockout mutants with gene interruption conferred by a kanamycin resistance cassette[35]. We screened all available knockout mutants representing these eight genes against a panel of five antitubercular drugs (plus the positive control kanamycin). We observed no resistance above the critical concentration to any of the tested drugs other than kanamycin, indicating that loss-of-function of each of these genes alone is insufficient to confer resistance to this particular panel of drugs (Supplementary Table 6).

One gene, ald (Rv2780), was notable in that it had an especially high number of unique loss-of-function mutations that appeared to evolve in a pattern consistent with convergent evolution. We identified 11 different ald loss-of-function mutations in 25 South African strains across four different spoligotypes, including one nonsense mutation and ten distinct indels (Figure 2A, Supplementary Table 7). When examined within a phylogenetic context, ald loss-of-function mutations appeared to evolve exclusively and repeatedly in MDR and XDR genetic backgrounds, suggesting a selection for loss-of-function of this gene among patients receiving second-line drug regimens (Figure 2B). Four other genes had loss-of-function mutations occur exclusively within an MDR and XDR genetic background (Rv1250, Rv2042c, Rv3197A: whiB7, and Rv3823c: mmpl8), although in all cases at a much lower frequency than ald (Supplementary Table 5). The role of loss-of-function mutations in these genes is unclear, as Rv1250 expression is induced by drug exposure[36], mmpl8 expression has been positively correlated with virulence[37], and overexpression of whiB7 has been shown to confer cross-resistance to aminoglycosides[38,39] while null mutants have been shown to be hypersusceptible to multiple antibiotics including macrolides[40,41] (Supplementary Note).

Figure 2

Convergent evolution of loss-of-function mutations in ald in MDR and XDR M. tuberculosis. A) Positions of insertions, deletions, and nonsense mutations in ald. Insertions are marked as blue triangles while deletions and nonsense mutations are marked as red triangles. Size of the indel or details of nonsense mutation are shown above or below each triangle. The location of the frameshift in BCG is also included. B) Phylogeny of lineage 4 strains. Inner tick marks denote MDR (gray) and XDR (black) M. tuberculosis. Outer tick marks denote insertions (blue) and deletions or nonsense mutations (red).

ald loss-of-function conferred resistance to D-cycloserine

The ald gene encodes L-alanine dehydrogenase, a component of the alanine metabolism pathway, which includes the targets of the second-line drug D-cycloserine (Figure 3)[42]. D-cycloserine resistance was not prospectively phenotypically assayed in the 498 strains analyzed here; however, D-cycloserine and terizidone, a dimeric form of D-cycloserine, are routinely prescribed as part of the South African drug-resistant M. tuberculosis treatment guidelines (see URLs). Therefore, MDR and XDR M. tuberculosis in this region could have developed D-cycloserine resistance due to selective drug pressure.

Figure 3

Alanine metabolism pathway in M. tuberculosis. L-alanine-dehydrogenase catalyzes the NAD-dependent interconversion between L-alanine and pyruvate. Under aerobic conditions this reaction allows for the utilization of L-alanine as a nitrogen source, whereas under hypoxic conditions this reaction allows for NADH recycling to NAD+. Separately, the initial step of peptidoglycan synthesis involves conversion of L-alanine to D-alanine by alanine racemase. Two D-alanine molecules are joined by D-ala-D-ala ligase to produce the dipeptide D-alanine-D-alanine, which is subsequently incorporated into the pentapeptide chain of the peptidoglycan cell wall. Both alanine racemase and D-ala-D-ala ligase are competitively inhibited by D-cycloserine, an analog of D-alanine.

To investigate whether ald loss-of-function conferred resistance to D-cycloserine, we determined the minimum inhibitory concentration (MIC) on Löwenstein-Jensen (LJ) media (Methods) for four strains: 1) CDC1551, a fully drug-susceptible reference strain of M. tuberculosis (with wild type alleles of ald, alr, ddlA, and cycA), 2) an ald knockout strain (Δald), 3) the ald knockout complemented with wild type ald (Δald-comp) from M. tuberculosis, and 4) M. bovis BCG Danish 1331 (BCG), a strain known to be naturally resistant to D-cycloserine[23,43]. Testing was repeated in multiple independent experiments to establish reproducibility of the assay; results are shown in Table 2. Wild type M. tuberculosis was the most susceptible to D-cycloserine (MIC = 15 μg/mL) while the positive control BCG was the most resistant (MIC = 40–60 μg/mL). The ald knockout had increased resistance relative to the wild type strain (MIC = 25–30 μg/mL), which overlaps the WHO-determined critical concentration for resistance on LJ of 30 μg/mL (see URLs). Complementation of ald partially restored susceptibility to D-cycloserine (MIC = 20 μg/mL). These results demonstrate that in vitro, ald loss-of-function confers a moderate level of resistance to D-cycloserine.

Table 2

Loss-of-function of L-alanine dehydrogenase (ald) confers increased resistance to D-cycloserine in vitro. The minimum inhibitory concentration (MIC) to D-cycloserine was determined for each strain on Löwenstein–Jensen media. Testing was done in three independent experiments (assays #1–#3) to establish reproducibility. Dashes indicate those strains were not tested during that experiment.

Strain	D-cycloserine MIC (μg/mL)
	Assay #1	Assay #2	Assay #3
wild type M. tuberculosis	15	15	-
Δald	25	30	25
ald-comp	20	20	20
BCG	40	40	60
BCG-ald-comp	40	-	40
M. bovis	-	25	25

These results were also confirmed in a growth advantage assay, in which all strains were cultured in a mycobacterial growth indicator tube (MGIT) system with varying concentrations of D-cycloserine (0, 7.5 and 15 μg/mL) and the resulting time to positivity in the presence of the drug was normalized to the no-drug control (Figure 4, Methods, Supplementary Figure 2). The ald knockout had a significant growth advantage relative to the wild type strain at 7.5 μg/mL D-cycloserine, (p < 0.0001, ANOVA), and no significant growth difference when compared to the resistant control strain, BCG, at the same drug concentration (Figure 4A). At a higher concentration of D-cycloserine (15 μg/mL), the ald knockout still retained a growth advantage over the wild type strain (Figure 4B; p < 0.0001, ANOVA), but grew more slowly relative to BCG (p < 0.0001, ANOVA), supporting the MIC data, which showed that ald loss-of-function confers less resistance than the mechanism present in BCG. At both drug concentrations, the ald complement partially reverted the growth advantage conferred by the ald knockout (Figure 4A, B; p < 0.0001, ANOVA).

Figure 4

A single gene knockout of ald confers a growth advantage relative to wild type M. tuberculosis when cultured in the presence of D-cycloserine. Four laboratory strains—wild type M. tuberculosis (WT), ald knockout (Δald ), ald complemented back into Δald (ald-comp) and BCG—were cultured in triplicate in 0, 7.5 and 15 μg/mL of D-cycloserine. The resulting time to positivity (ttp) in MGIT was recorded as days since inoculation, and the resulting ttp was normalized to the no drug control for each strain to calculate the growth inhibition index for each drug concentration. The mean growth inhibition index is plotted with error bars to represent the SEM. P-values were calculated using a two-way ANOVA.

Within the M. tuberculosis complex (MTBC), BCG is uniquely considered resistant to D-cycloserine[44]. This resistance was previously shown to be, at least partly, due to a unique SNP in cycA, which encodes an alanine transporter[24,45]. However, all of the RD9 lineages of MTBC from M. africanum to M. bovis, including BCG, have a frameshift in ald causing loss-of-function (Supplementary Figure 3)[46,47]. We therefore sought to evaluate the effect of the frameshift in both the BCG and more general RD9 genetic background by determining the MIC of two additional strains: 1) M. bovis ATCC 19210, a reference strain of M. bovis, and 2) M. bovis BCG complemented with M. tuberculosis ald (BCG-ald-comp). M. bovis ATCC 19210 had a D-cycloserine MIC of 25 μg/mL (Table 2), which is similar to the level of resistance observed in Δald and greater than the level of resistance observed in the M. tuberculosis reference strain, CDC1551. The ald-complemented BCG strain had an MIC of 40 μg/mL, which was not different from the wild type BCG strain (Table 2), although in a growth advantage assay at low concentrations of D-cycloserine, BCG had a significant growth advantage when compared to the ald-complemented BCG strain, (p < 0.02 and p < 0.0001; Supplementary Note, Supplementary Figures 4 and 5). Taken together, these results suggest that ald loss-of-function may have both a minor role in the resistance of BCG to D-cycloserine, and provide M. bovis with an elevated resistance to D-cycloserine relative to M. tuberculosis.

ald and alr in D-cycloserine resistant clinical isolates

To more specifically ascertain the distribution and diversity of mutations in ald, along with mutations in alr, ddlA, and cycA, we characterized all promoter, nonsynonymous, and loss-of-function mutations in these four genes across all strains in our dataset (Table 3, Supplementary Table 8). Mutations were most frequent in ald (25 strains), followed by alr (17 strains). Only seven strains were identified to have mutations in either ddlA or cycA; none of these mutations corresponded to previously characterized resistance-conferring mutations. We observed that many more strains had ald loss-of-function mutations than had ald nonsynonymous mutations (24 versus 5, respectively). By contrast, mutations in alr were exclusively nonsynonymous, with the exception of one strain with a SNP in the promoter region, supporting the hypothesis that alr is essential[48].

Table 3

Counts of clinical isolates with mutations in ald, alr, ddlA, and cycA. Mutation types include wild type (no mutations; WT) loss-of-function (LOF), nonsynonymous (NSY), and promoter SNPs (PRO); synonymous SNPs are excluded. Specific strains and mutations are shown in Supplementary Table 8.

# of Clinical Isolates	ald	alr	ddlA	cycA
449	WT	WT	WT	WT
20	LOF	WT	WT	WT
5	NSY	WT	WT	WT
4	LOF	NSY	WT	WT
12	WT	NSY	WT	WT
1	WT	PRO	WT	WT
1	WT	WT	NSY	WT
1	WT	NSY	WT	NSY
5	WT	WT	WT	NSY

We next sought to verify our prediction that clinical strains with ald loss-of-function mutations would exhibit resistance to D-cycloserine, and to determine how that level of resistance compared to strains with mutations in alr. We selected 44 clinical isolates from South Africa for testing, including 22 with wild type ald and alr, 13 with mutations in ald and wild type alr, 7 with wild type ald and mutations in alr or its promoter, and 2 with mutations in both ald and alr (Supplementary Table 9). These strains had diverse drug susceptibility phenotypes and varying genetic backgrounds (Supplementary Table 9).

Clinical isolates with ald mutations showed significantly increased resistance to D-cycloserine relative to strains with wild type alleles of both ald and alr (Figure 5A, Supplementary Table 9; p < 0.0002, Mann-Whitney U test). Strains with ald loss-of-function mutations had a median MIC of 40 μg/mL, with 8 of those 11 strains exhibiting an MIC from 30–40 μg/mL. The two strains with nonsynonymous SNPs in ald had MICs of 30 μg/mL; both of these nonsynonymous SNPs were predicted to have a deleterious effect on protein function by PROVEAN[49] (Supplementary Table 9). By contrast, strains with wild type alleles in both genes had a median MIC of 20 μg/mL, with 91% of strains exhibiting MICs <= 20 μg/mL. These data suggest that both deleterious nonsynonymous SNPs and loss-of-function mutations in ald can cause resistance to D-cycloserine in clinical isolates of M. tuberculosis at or above the current critical concentration (30 μg/mL; URLs).

Figure 5

Clinical strains with mutations in L-alanine dehydrogenase (ald) and alanine racemace (alr) exhibit increased resistance to D-cycloserine. A set of 44 clinical strains of M. tuberculosis with differing genotypes with respect to ald and alr were selected, representing wild type alleles (WT), loss-of-function mutations (LOF), nonsynonymous mutations (NSY) and promoter mutations (PRO). The minimum inhibitory concentration (MIC) to D-cycloserine was determined for each strain on Löwenstein–Jensen media. For wild type strains and strains with ald loss-of-function mutations, the median is shown as an open circle and quartiles are shown as connected bars. When compared to WT isolates, strains with ald loss-of-function had significantly increased resistance to D-cycloserine (p < 0.0002, Mann-Whitney U test).

Clinical isolates with alr mutations also showed increased resistance to D-cycloserine relative to strains with wild type alleles of both ald and alr (Figure 5A, Supplementary Table 9). Of the strains with nonsynonymous mutations in alr, 4 of 6 exhibited high level resistance to D-cycloserine, with MICs of >=60 μg/mL, while one, plus the two strains with both alr and ald mutations, exhibited MICs of 40 μg/mL. As with ald, all of the nonsynonymous SNPs seen in alr were predicted to have a deleterious effect on protein function by PROVEAN[49] (Supplementary Table 9). Additionally, the one strain with an alr promoter mutation showed high level resistance (MIC > 60 μg/mL), consistent with previous results in M. smegmatis that showed upregulation of alr through a promoter mutation conferred resistance to D-cycloserine[19]. These results suggest that both nonsynonymous SNPs in alr and SNPs in the alr promoter can confer high-level resistance to D-cycloserine in clinical isolates of M. tuberculosis, and that alr and ald mutations are non-additive for resistance.

Discussion

To identify novel drug resistance-conferring mechanisms in M. tuberculosis, we analyzed whole-genome sequencing data from two large collections of clinical isolates from South Africa and China with diverse phenotypic drug resistance profiles. Combining a direct association test with a correlated evolution test followed by functional screening of the top scoring genes with loss-of-function mutations, we identified a novel drug resistance mechanism to D-cycloserine, ald loss-of-function, which we then validated experimentally. We further showed that clinical strains with mutations in ald and alr, identified through our computational analyses, were resistant to D-cycloserine. Collectively, this work expands our understanding of D-cycloserine resistance in M. tuberculosis and provides new opportunities to detect and manage D-cycloserine resistant TB.

The strong linkage disequilibrium and stratification present in bacterial populations has long been highlighted as a major hurdle for bacterial genome-wide association studies[50]. However, the problem of interdependence is not restricted to genotypes, but extends to drug resistance phenotypes as well[26]. For most TB drugs, mutations known to confer resistance to a single drug were significantly associated with phenotypic resistance to numerous unrelated drugs. This is because TB is treated exclusively with multi-drug regimens, which makes it difficult to collect clinical isolates of M. tuberculosis with independent drug resistance phenotypes. However, combining a basic association test with a correlated evolution test appeared to mitigate this phenotypic interdependence, as it could often correctly match known resistance mechanisms with the correct drug; the method did have difficulty when the resistance-conferring genotype evolved infrequently (e.g., pyrazinamide; Supplementary Note) and had other limitations (Supplementary Note). Ultimately, to robustly overcome phenotypic interdependence in association studies of drug resistant M. tuberculosis, further methods development is needed to build upon this and other work[25,26].

After we identified loss-of-function mutations in ald as a candidate novel drug resistance mechanism, we showed experimentally that ald loss-of-function confers resistance to D-cycloserine in M. tuberculosis. Complementation of ald back into the ald knockout strain partially restored susceptibility to D-cycloserine (partial restoration by complementation has been observed before in M. tuberculosis, and could be due to altered expression levels through usage of a non-native promoter[51]). We hypothesize that strains lacking functional Ald are unable to convert L-alanine to pyruvate, thereby increasing the pool of available L-alanine (Figure 3). As L-alanine is the precursor to the pathway steps competitively inhibited by D-cycloserine, abundant L-alanine in the setting of non-functional Ald may allow for continued alanine racemization and D-alanine ligase activity, and therefore continued peptidoglycan production despite competitive inhibition by D-cycloserine. However, the observation in clinical strains that ald loss-of-function conferred a lower level of resistance than nonsynonymous or promoter mutations in alr (Supplementary Table 9), the direct target of D-cycloserine, suggests that the resistance conferred by ald could more easily be overcome through increased dosage of D-cycloserine.

M. bovis and other members of the MTBC that lack RD9[52] share a frameshift in ald. We sought to determine the effect of this frameshift in M. bovis strains traditionally classified as D-cycloserine resistant (BCG) and susceptible (all other strains)[44]. We observed that the M. bovis reference strain, ATCC 19210, had elevated levels of D-cycloserine resistance as compared to clinical isolates of M. tuberculosis lacking resistance-conferring mutations. Our results suggest that M. bovis, and potentially the entire RD9 lineage, which is generally considered D-cycloserine susceptible, may actually have elevated resistance levels (Figure 5B, URLs). However, much of the RD9 lineage (including M. bovis, but excluding BCG) also shares an inactivating mutation in the gene encoding pyruvate kinase (pykA), which provides pyruvate to Ald[53,54] and could, therefore, potentially confound analysis of D-cycloserine resistance among these strains. Evolutionary analyses suggest that in animal-infecting species of the MTBC there is more ancient functional diversity of these genes than in M. tuberculosis sensu stricto (Supplementary Note, Supplementary Figure 3). Future experiments to disentangle the roles of ald, pykA and other related genes in drug resistance across the other species of the MTBC will be needed.

The identification of so many unique ald loss-of-function mutations among drug-resistant clinical strains of M. tuberculosis suggests that, despite conferring only a moderate level of resistance, this mechanism provides sufficient resistance to be selected for under current drug exposure pressures in South Africa (Supplementary Note). This further suggests that ald loss-of-function mutations confer resistance to the D-cycloserine analog, terizidone, which is routinely prescribed for treatment of drug-resistant M. tuberculosis in South Africa (see URLs), although this hypothesis still needs to be tested. While ald loss-of-function may prove advantageous during exposure to D-cycloserine or terizidone, harboring such a mutation may prove to be deleterious in other settings. Strains with non-functional Ald may be unable to increase pyruvate availability under hypoxic conditions, and further laboratory investigation should determine whether ald loss-of-function yields decreased survival in hypoxia. However, compared to the M. bovis genetic background where pykA is also non-functional, clinical M. tuberculosis strains with loss-of-function of ald may have greater ability to tolerate hypoxic environments (such as in a granuloma) where reliance on pyruvate (rather than aerobic metabolism) is greater.

Given the severe psychiatric and central nervous system toxicities of D-cycloserine[9,10,12], it is critical to both incorporate D-cycloserine resistance mechanisms into molecular diagnostics and optimize phenotypic resistance assays to prevent the drug from being inappropriately prescribed to patients who harbor resistant strains. In order to maximize the effectiveness of phenotypic resistance assays, the range of wild type MICs for an organism must be determined to define drug resistance and to set appropriate critical concentrations[55,56]. Recently, the Stop TB Partnership recommended to the WHO that the critical concentration of D-cycloserine be lowered from 40 to 30 μg/mL (URLs). Assuming that the two isolates exhibiting increased resistance to D-cycloserine, but having wild type alleles in both ald and alr, harbor an unknown resistance mechanism, our results suggest a tentative epidemiological cutoff (ECOFF)[57] of 20 μg/mL (Figure 5B). Although phenotyping of a larger number of strains is required, our results support a lowering of the critical concentration below 30 μg/mL. Selection of a critical concentration would also need to consider whether it should be applied to M. tuberculosis sensu stricto or the entire MTBC including M. bovis.

D-cycloserine resistance is an ideal target for molecular diagnostics, given technical difficulties with reproducibility and reliability of the phenotypic resistance assay[13-15]. However, current SNP-based diagnostics are not amenable to capturing loss-of-function mutations, as the mutations can occur anywhere within the gene and thus are difficult to capture with a small number of targeted SNPs. Other whole-gene-based diagnostic strategies will instead need to be employed, some of which have already been developed for detection of drug-resistant TB, such as high-resolution melt (HRM) analysis[58] and rapid whole-genome sequencing of clinical isolates[21]. Incorporation of resistance mechanisms to D-cycloserine and other highly toxic drugs into novel molecular diagnostics could prevent those drugs from being prescribed in a treatment setting where they may be ineffective and needlessly toxic.

Online Methods

Bioinformatic Analyses

We used two sets of strains in this study: one set of 337 strains from South Africa[28] (NCBI Bioproject ID numbers: PRJNA183624, PRJNA235615) and one set of 161 strains from China[26] (NCBI Sequence Read Archive (SRA) accession number: SRA065095). For both sets, reads were mapped onto a reference assembly of H37Rv (GenBank accession number: CP003248.2) using BWA version 0.5.9[60]. In cases where read coverage of the reference was greater than 200x, reads were down-sampled using Picard (see URLs) prior to mapping. Variants were identified using Pilon version 1.5 as described[61]. All sites with unambiguous SNPs were used to generate a phylogenetic tree using FastTree version 2.1.3[62]. To assess the amount of unexplained phenotypic resistance in our dataset, we identified genomic polymorphisms known to confer drug resistance in each strain using variants from current diagnostics supplemented with a literature review (Supplementary Note).

In order to identify novel variants conferring drug resistance, we hypothesized that the presence of causal genotypes in bacteria should be both positively correlated with phenotypes and not evolutionarily independent of them. We therefore conducted two tests, the first of which was a Fisher’s exact test to measure the association between genotype and phenotype. The second test was a test for correlated evolution[34] implemented in the software package BayesTraits (see URLs). Briefly, the correlated evolution test calculates the likelihood of two traits evolving dependently and independently and then compares the two models with a likelihood ratio test. If the dependent model is significantly more likely than the independent model, the hypothesis of evolutionary independence is rejected. All test results were corrected for multiple comparisons using the Benjamini-Hochberg method[59] for plotting purposes. Once both the Fisher’s exact test and the correlated evolution test were calculated, their uncorrected p-values were combined into an F-score by taking the harmonic mean of the negative log p-values. Uncorrected values were used in the F-score calculation to provide maximal discrimination between variants with poor p-values.

In order to compare phenotypic drug resistances to genotypic features, two matrices of these features were constructed from the variant calls of all strains. In the first matrix, rare variants (those at ≤ 2% frequency) were collapsed by gene or intergenic region, while common variants (those at > 2% frequency) remained uncollapsed and were considered independently. In the second matrix, we assessed whether each gene had a loss-of-function mutation (defined as a nonsense mutation, frameshift, or large (>30 bp) in-frame insertion or deletion (indel)), regardless of the frequency of the mutation. We then combined the two matrices into a single matrix for association analyses. All analyses were conducted twice: once with all samples, and once after removing strains with known resistance-conferring mutations, in order to amplify signal from strains with potential for novel resistance genotypes (Results).

DST screen of transposon mutants

Single gene knockout transposon mutants were sourced from the Johns Hopkins TARGET mutant collection[35]. First and second-line DST was performed by the agar proportion method on 7H11 media, using the following WHO-recommended critical concentrations (in μg/mL): isoniazid 0.2, isoniazid 1.0, rifampicin 1.0, ethambutol 7.5, streptomycin 2.0, ofloxacin 2.0 and kanamycin 6.0.

ald complementation into CDC1551 and BCG Danish 1331

An ald knockout mutant (Δald) of M. tuberculosis (JHU2780-209) was obtained from the Johns Hopkins TARGET mutant collection[35]. This single gene knockout mutant was created in a background of wild type M. tuberculosis strain CDC1551 by Himar1 transposon mutagenesis with ald interrupted at nucleotide 209, which contains a kanamycin-resistance cassette. To complement Δald, a 1.116kb fragment spanning the ald gene was PCR amplified from M. tuberculosis H37Rv genomic DNA using primers listed in Supplementary Table 10. The fragment was directionally cloned into pMCZ-MOP, an L5 integrating zeocin-marked mycobacterial vector to create pMCZ-MOP-ald, which allowed for the expression of ald off the MOP (mycobacterial optimized promoter)[63]. The pMCZ-MOP-ald was then transformed into Δald to create the Δald complemented strain, resulting in constitutive expression of ald. Transformants were selected using the appropriate antibiotic markers and confirmed by amplicon sequencing of a 1.319kb fragment of the integrated plasmid spanning the inserted ald gene. pMCZ-MOP-ald was also transformed into BCG Danish 1331 (Statens Serum Institut) to create the BCG-ald-complement.

M. bovis amplicon sequencing

Amplicon sequencing of ald and pncA[64] in M. bovis ATCC 19210 confirmed the presence of the ald mutation (a266del), and a mutation in pncA (g169c) that is common to all M. bovis and BCG strains.

D-Cycloserine MIC determination

Select mycobacterial strains were subjected to D-cycloserine MIC determination on Löwenstein–Jensen (LJ) medium. LJ slants were prepared freshly as had been previously described per WHO protocols, with supplementation of D-cycloserine prior to inspissation to achieve final drug concentrations of 0μg/mL, 10μg/mL, 15μg/mL, 20μg/mL, 25μg/mL, 30μg/mL, 40μg/mL, 60μg/mL, and 120μg/mL. 103 CFU of a mid to late-log liquid culture of M. tuberculosis was inoculated onto D-cycloserine slants at each drug concentration in duplicate. Slants were incubated at 37°C for four weeks, after which CFU were counted. The MIC to D-cycloserine for each strain was determined to be the lowest drug concentration that resulted in 90% reduction in CFU when compared to the no drug control. Differences in MIC distributions were calculated using the Mann-Whitney U test; in cases where MICs were replicates, the average interval across replicates was used for the calculation.

D-cycloserine growth inhibition assay

All mycobacterial strains were grown under standard conditions, in Middlebrook 7H9 broth (Difco) supplemented with 0.5% glycerol, 10% oleic acid-albumin-dextrose-catalase (OADC), and 0.05% Tween 80 with constant shaking at 37°C. On the day of the experimental set up, D-cycloserine powder (C6880, Sigma-Aldrich) that had been stored at −20°C was freshly dissolved in a 0.1M sodium phosphate buffer (pH 8.0) and sterilized with a 0.22μm filter. BD BACTEC™ MGIT™ tubes were prepared with BD BACTEC™ MGIT™ growth supplement, and 6.5mL tween-free 7H9 media supplemented with or without D-cycloserine to achieve a final drug concentration of 0 μg/mL, 7.5 μg/mL, 15 μg/mL and 30 μg/mL of D-cycloserine. 103 CFU of a mid to late-log liquid culture of M. tuberculosis that had been prepared from a single colony isolate of each strain was inoculated into three MGIT tubes for each drug concentration to achieve three technical replicates. Dilutions of the initial inoculum were plated for CFU in parallel to quantify the starting inoculum. The bar code of each MGIT tube was scanned into the BD BACTEC™ MGIT™ 960 Mycobacterial Detection System, and the MGIT tubes were placed in the instrument as directed. Samples remained in the MGIT instrument until they were flagged positive for growth, and the time to positivity (ttp) was recorded. Tubes that were not flagged positive by 42 days were considered to have no growth, and were assigned a ttp value of > 42 days. Data were reported for each strain and each drug concentration as a mean ttp for the three technical replicates ± the standard error of the mean (SEM). To determine relative growth inhibition of each drug concentration for each strain, the ttp of each D-cycloserine-containing MGIT was normalized to the mean ttp of no drug control for that strain. The mean ttp of the no drug control was considered to have a growth inhibition value of 1. Growth inhibition was calculated for each MGIT tube and reported as a mean growth inhibition for the three technical replicates ± the SEM. Adjusted p-values were calculated using a two-way ANOVA with Tukey’s multiple comparisons test in Prism GraphPad version 6.0 (GraphPad, San Diego, CA). We also confirmed that the transposon mutagenesis did not confer a growth advantage in the presence of D-cycloserine (Supplementary Note, Supplementary Tables 11 and 12).