Genomic analysis of globally diverse Mycobacterium tuberculosis strains provides insights into the emergence and spread of multidrug resistance.

Multidrug-resistant tuberculosis (MDR-TB), caused by drug-resistant strains of Mycobacterium tuberculosis, is an increasingly serious problem worldwide. Here we examined a data set of whole-genome sequences from 5,310 M. tuberculosis isolates from five continents. Despite the great diversity of these isolates with respect to geographical point of isolation, genetic background and drug resistance, the patterns for the emergence of drug resistance were conserved globally. We have identified harbinger mutations that often precede multidrug resistance. In particular, the katG mutation encoding p.Ser315Thr, which confers resistance to isoniazid, overwhelmingly arose before mutations that conferred rifampicin resistance across all of the lineages, geographical regions and time periods. Therefore, molecular diagnostics that include markers for rifampicin resistance alone will be insufficient to identify pre-MDR strains. Incorporating knowledge of polymorphisms that occur before the emergence of multidrug resistance, particularly katG p.Ser315Thr, into molecular diagnostics should enable targeted treatment of patients with pre-MDR-TB to prevent further development of MDR-TB.

Introduction

Drug-resistant Mycobacterium tuberculosis is a threat to global tuberculosis (TB) control efforts. Failure to identify and appropriately treat patients with drug-resistant TB can lead to increased mortality, nosocomial outbreaks, and the expansion of drug resistance[1]. Five percent of M. tuberculosis cases worldwide are now multidrug-resistant (MDR), defined as having resistance to both isoniazid and rifampicin[2]. Therapeutic regimens for MDR-TB can exceed 18 months, and include agents that often confer significant adverse effects[3]. 0.5% of global TB cases are now considered extensively drug-resistant (XDR), defined as MDR with additional resistance both to fluoroquinolones and at least one second-line injectable drug[2]. XDR-TB has incredibly poor treatment outcomes; in one long-term cohort of XDR patients in South Africa, only 19% of patients had a favorable outcome[4].

The global front-line molecular diagnostic for drug-resistant M. tuberculosis, GeneXpert MTB/RIF (GeneXpert, Sunnyvale, California, USA, see urls), simultaneously detects the presence of M. tuberculosis and identifies rifampicin resistance[5]. While GeneXpert identifies patients harboring rifampicin-resistant strains for initiation of MDR-TB treatment, this test may not identify resistance at the earliest available opportunity. In a recent analysis of a large collection of M. tuberculosis clinical isolate genomes from South Africa[6], Cohen et al. showed that the overwhelming majority of MDR- and XDR-TB evolved resistance to isoniazid prior to resistance to rifampicin. This result was consistent with another recent genomic analysis of strains from Russia[7] and an MDR-TB outbreak in Argentina[7]. In addition, analysis of phenotypic drug susceptibility tests from a large, global collection of strains collected during TB drug resistance surveys indicated that isoniazid resistance is acquired before rifampicin resistance[8].

Collectively, these results suggest that, in order to detect resistance as soon as possible and to prevent MDR- and XDR-TB from evolving, molecular diagnostic tests for M. tuberculosis should include the earliest resistance mutations to emerge; however, the identities of these MDR “harbinger” mutations remain undefined. To close this gap in our understanding, we have undertaken a large-scale analysis of a global dataset of 5,310 M. tuberculosis whole-genome sequences, including 868 newly sequenced strains, and 4,442 previously published strains, to determine the ordering of acquisition of drug resistance mutations, and to identify which mutations occur early along the pathway toward MDR and might serve as early sentinels in the development of MDR.

Results

Drug resistance arises by universal mechanisms across the globe

In order to examine global phylogeographic patterns, including the order of drug resistance mutation evolution in M. tuberculosis, we compiled a set of 8,316 whole genome sequences from clinical M. tuberculosis strains that were either newly sequenced as part of this study or were part of fourteen published studies that used Illumina technology (Supplementary Table 1)[6,9-22]. After quality filtering (see Methods), our dataset included 5,310 genome sequences representing M. tuberculosis strains from 48 countries and 17 United Nations (UN) defined geographic regions[23] (Supplementary Note A; Supplementary Figures 1–4; Supplementary Tables 1–3). Though our dataset represented a broad diversity of TB strains from many global regions, their phylogeographic distribution did not perfectly match the actual distribution of TB burden worldwide (Supplementary Figure 1b); however, all seven known global lineages of M. tuberculosis[24] were represented (lineage 1, EAI or Indo-Oceanic lineage; lineage 2, Beijing lineage; lineage 3, CAS or Central Asian lineage; lineage 4, Euro-American lineage; lineage 5, M. africanum West African Type I; lineage 6, M. africanum West African Type II; and the deep-branching lineage 7), as well as M. bovis. Unsurprisingly, lineages 1–4 predominated (99.2%), consistent with the previously described limited geographic and host distributions of lineages 5–7 (Supplementary Table 3)[21,25,26].

To examine the distribution of drug resistance in this sample, for each isolate genome we computationally predicted resistance to eight drugs[27] using a curated list of polymorphisms associated with resistance (see Methods; Supplementary Table 4). As phenotypic drug resistance information was unavailable for most datasets, we did not incorporate phenotypic information into our analysis. We identified a total of 392 unique drug-resistance polymorphisms in at least one strain (Supplementary Tables 2,5). In comparison to expected global resistance rates, we observed higher rates of resistance, with 962 strains (18%) carrying mutations for both rifampicin and isoniazid resistance and lacking mutations for ofloxacin and kanamycin resistance (MDR sensu stricto), and 257 (5%) of strains carrying mutations for resistance to all four drugs that define XDR (rifampicin, isoniazid, ofloxacin, and kanamycin) (Supplementary Table 6). Another 409 (8%) strains carried mutations causing pre-XDR level resistance (MDR genotype plus mutations to either ofloxacin or kanamycin). Over half of the sequenced strains did not carry any resistance-conferring mutation so were predicted to be drug susceptible (Figure 1; Supplementary Figures 4–6; Supplementary Table 2).

Figure 1

A) Geographic distribution of M. tuberculosis isolates in our dataset by drug resistance pattern. This plot shows the distribution of the 5,310 M. tuberculosis isolates included in our dataset by drug resistance genotype (pie charts) and by 11 UN geographic subregions[23] (coloring), and is not meant to indicate the overall global incidence of TB or drug resistance. There were no strains in our dataset from geographic regions colored grey. UN geographic subregions with fewer than 30 strains were excluded from this figure. Map modified from a blank map of UN geographical subregions.[64] B) The overall proportion of drug-resistant strains identified among all 5,310 M. tuberculosis isolates in our dataset.

Drug-resistance was identified in nearly all (15 of 17) UN regions for which we had data, though its regional distribution varied considerably (Supplementary Note A; Supplementary Figures 7–10). In certain regions of the globe, we observed high numbers of closely related strains with nearly identical sets of resistance mutations, which can be attributed to clonal transmission. Because our dataset contained isolates from several known outbreaks[6,12], rather than focusing on the total number of strains with each mutation, we instead examined the number of times each mutation evolved in different global regions by counting independent arisals, or the number of separate evolutions of a specific mutation occurring at defined positions in the phylogeny for a specific geography. Using parsimony-based analysis to reconstruct mutation gains and losses at all nodes across the phylogeny (see Methods), we observed that the distribution of arisals of specific mutations was fairly constant across the globe, in contrast to the uneven distribution of strains with these mutations (Supplementary Note B; Supplementary Tables 7–8), suggesting that drug resistance has arisen via similar mechanisms irrespective of geography. This was also true for the evolution of MDR and XDR, which we calculated evolved independently 573 and 138 times within our dataset, respectively. Along with frequent, repeated, de novo arisals, person-to-person transmission – as predicted when all strains descending from a common ancestor in the phylogeny shared the same MDR genotype – was also an important contributor to the observed MDR cases. 360 of the 573 arisals (63%) led to a single MDR strain in our dataset (de novo evolution), whereas 213 arisals (37%) resulted in two or more descendant MDR strains, likely indicating person-to-person transmission of MDR-TB.

Isoniazid resistance overwhelmingly arises before rifampicin resistance, across all lineages, geographic regions, and time

In an earlier analysis involving a smaller, South African dataset[6], we showed that isoniazid resistance evolved prior to rifampicin resistance in almost all cases. To determine whether this ordering of mutation acquisition held in a globally diverse set of strains, we used a parsimony-based analysis to examine the order of pairwise arisals of drug resistance mutations. We filtered portions of the phylogeny with ambiguous topologies (see Methods) and only included nodes at which explicit ordering could be established[6]. In agreement with our previous results[6], we found that resistance to first-line drugs generally evolved prior to resistance to second-line drugs (Supplementary Note C), as would be expected from the order in which antituberculous drugs are utilized in clinical practice. We also observed that mutations conferring isoniazid resistance overwhelmingly arose prior to any other mutation implicated in resistance (Figure 2; Supplementary Table 9; and Supplementary Note C), despite substantial complexity in the types and ordering of evolved mutations (Figure 3). Strikingly, isoniazid resistance predated rifampicin resistance in 96% of pairwise comparisons (155 out of 162), a pattern that remained true regardless of lineage or geographic source (Figure 4). While the majority of this effect was due to katG mutations (98% in 114 out of 116 pairings), non-katG isoniazid mutations followed this same pattern (89%, 41 out of 46 pairings). Thus, the provenance of global MDR was overwhelmingly isoniazid-resistant strains. In particular, strains carrying the katG S315 mutation frequently gained rifampicin resistance, whereas only a very small minority (4%) of global MDR arisals were due to gain of isoniazid resistance within a rifampicin-resistant background, despite the presence of 48 rifampicin mono-resistant strains, and 152 (3% of total) isolates containing rifampicin resistance, but no isoniazid resistance.

Figure 2

Across the globe, isoniazid resistance was overwhelmingly the first step towards drug resistance. Acquisition of a katG S315 mutation precedes all other resistance mutations for the majority of instances in which the order of acquisition can be disambiguated. We quantified the pairwise number of evolutions in which resistance to one drug preceded resistance to a second drug. Reported numbers represent the number of independent evolution events (not the number of strains) in which the drug resistance indicated by the row labeled “first resistance” was acquired before the drug resistance indicated by the column labeled “second resistance”. The shading color indicates the percentage of evolutionary events in which the “first resistance” clearly predates the “second resistance” for that drug pair. While inhA mutations can confer resistance to both isoniazid and ethionamide[65], we defined genotypic ethionamide resistance as mutations in only ethA to simplify the analysis and avoid double counting.

Figure 3

Sequential acquisition of drug resistance mutations reveals that isoniazid resistance-conferring mutations, specifically katG S315T, most often come first in sequential pairs. This figure includes data from 71 mutations conferring drug resistance with at least 10 occurrences in our dataset, which represent 93% of all drug resistance mutations in our dataset. Using PAUP analysis to assign specific mutation gains to individual nodes on the phylogeny, we tabulated all routes of drug resistance acquisition across the full strain phylogeny, examining only those nodes on the tree where drug resistance mutations arose (i.e. node 1 [mutation A] -> node 2 [mutations B and C] -> node 3 [mutation D]). We tabulated the number of times each pair of mutations arose sequentially at adjacent nodes (i.e. mutations A->B, A->C, B->D, and C->D). We removed node pairs that did not meet specific bootstrap and branch length criteria (see Methods). The ribbons in this figure depict the number of times that each pair of mutations arose sequentially at adjacent nodes across the entire dataset. The width of the ribbon at each end is proportional to the number of times mutation A arose before mutation B, or vice versa (i.e., a ribbon with a thick end at katG S315T and a thin end at rpoB S450L indicates that katG S315T mutation arose prior to rpoB S450L much more frequently than the opposite). Each ribbon is colored according to the mutation more often occurring first in each sequential pairing.

Figure 4

In all lineages and global regions, the katG mutation S315T occurs first, and few examples of the reverse ordering are observed. We separately recalculated phylogenies for isolates from patients in each of the 11 UN subregions and 5 lineages with greater than 30 representatives (see Methods). This figure depicts the pairwise ordering of the katG S315T mutation in relation to mutations conferring resistance to the other three XDR-definings drugs (rifampicin (R), kanamycin (K), and ofloxacin (O)), within each individual M. tuberculosis lineage and geographic region. The numbers in panel B do not necessarily add up to the numbers in panel A, as the analyses of regions and lineages were performed individually, which can affect the number of arisals. Grey shading indicates that there were not sufficient pairings. Data are not shown for the following regions and lineages, as there were insufficient pairings: West Africa, Southern Europe, Central Asia, Northern America, lineage 1, and M. bovis.

One possible explanation for this striking result is that isoniazid entered into clinical use approximately 20 years before rifampicin (1971–1993, depending on geography)[28] resulting in ancestral M. tuberculosis populations having differing exposures to these drugs, and potentially impacting the ordering of acquisition of resistance in favor of isoniazid resistance before rifampicin resistance. To test this hypothesis, we predicted the date for the arisal of each isoniazid and rifampicin resistance mutation using BEAST[29], and then tallied the number of co-arisals of resistance to both drugs occurring at various posterior time cutoffs to the present (see Methods) starting at 1971 (the date of introduction of rifampicin) and ending at 2000 (a later date to account for lag in the timing of rifampicin’s widespread use) (Supplementary Tables 9–11). Our results revealed that, regardless of the date increment or evolutionary rate chosen (Methods; Supplementary Note D), katG S315-mediated resistance arose prior to rifampicin resistance 92%–98% of the time (Supplementary Table 10), indicating that, even during the era when isoniazid and rifampicin were given in combination, isoniazid resistance predated rifampicin resistance.

Diagnostics for early detection of pre-MDR M. tuberculosis

Contributors to the current global burden of MDR-TB include not only historic emergences of MDR, which led to person-to-person transmission of MDR-TB, but also ongoing de novo evolution. Of the 573 MDR arisals in our dataset, we estimated that 67% occurred since 2004 (see Methods). Thus, new strategies for curbing the emergence of MDR, such as identifying MDR precursors, will be critical to the control of MDR worldwide. GeneXpert, currently one of the front-line diagnostic tests used to exclusively identify rifampicin resistance-conferring mutations in the RRDR of rpoB[5], is commonly used globally as a proxy for detecting MDR-TB. The habitual ordering of isoniazid resistance prior to rifampicin resistance indicates that GeneXpert MTB/RIF serves as an appropriate proxy for MDR-TB, and is well suited to detect MDR in all geographic regions and all lineages of M. tuberculosis (Figure 4). However, as rifampicin mutations (detectable by GeneXpert) are rarely the first resistance mutations to emerge, by the time a mutation that is detectable by GeneXpert develops, oftentimes there is pre-existing resistance to multiple additional drugs, including second-line drugs (Figure 5). As we excluded nodes at which we were unable to disambiguate the ordering of resistance, our estimates represent a lower bound on the number of nodes where resistance to other drugs was gained prior to a GeneXpert-detectable node (Supplementary Figure 11).

Figure 5

Non-rifampicin drug resistance often precedes the arisal of GeneXpert mutations. Data are shown here for the nodes at which a GeneXpert mutation arose. A) This plot shows the percentage of GeneXpert nodes where resistance to each of eight drugs unambiguously preceded the arisal of GeneXpert mutations. Drug resistances that appeared to arise coincident to the GeneXpert node were excluded from this representation. More than one additional resistance could precede a single GeneXpert node. No strains contained additional rifampicin mutations arising before GeneXpert mutations. The bars represent a lower bound on the percentage of GeneXpert mutations preceded by additional resistance mutations, as we were unable to disambiguate ordering for a substantial number of nodes where additional mutations arose at the same node (see Supplementary Figure 11). B) Percentage of nodes where resistance to one or more other drugs unambiguously preceded the arisal of GeneXpert mutations. 13% of GeneXpert arisal nodes unambiguously had no additional drug resistance mutations arising prior.

Diagnostics that identify mutations present before MDR emerges would provide an opportunity to identify drug resistance during a period where there are both greater therapeutic options and improved treatment outcomes[30]. While the results from our pairwise ordering of resistance acquisition clearly demonstrated katG S315T-mediated isoniazid resistance as among the earliest to evolve, our pairwise approach necessarily oversimplified the complex process of MDR-TB evolution. Thus, in a complementary approach to identify other possible sentinels of complex resistance, we cataloged all of the resistance mutations that commonly evolved prior to development of MDR-TB (“pre-MDR-TB mutations”). For this set of mutations, we quantified the fraction of MDR-defining nodes where one of these pre-MDR mutations had evolved prior to the development of MDR-TB in order to determine how much resistance to other drugs had unambiguously arose prior to MDR (Figure 6). Our analysis (see Methods) revealed a set of 16 resistance mutations (of 340 total found among MDR and XDR strains) that arose prior to MDR-TB a minimum of two independent times (Supplementary Note E; Supplementary Table 12).

Figure 6

katG S315T is a commonly occurring mutation with very little resistance to other drugs arising prior to its occurrence. A) For each of the 16 “pre-MDR” mutations, the percentage of nodes where resistance to another drug unambiguously preceded the arisal of that mutation. These bars represent a lower bound on the percentage of nodes preceded by another resistance mutation, as we were unable to disambiguate ordering for a substantial number of nodes where additional mutations arose at the same node. B) The number of independent arisals for each of 16 “pre-MDR” (or harbinger) mutations. Since there are two embB M306I mutations, the nucleotide change at position 4,247,609 is also indicated for these two variants.

Surprisingly, we observed resistance mutations for all eight drugs among this set of pre-MDR mutations. However, many of these mutations evolved infrequently (Figure 6a), and thus would likely have low negative predictive value if included on a diagnostic panel aimed at identifying pre-MDR-TB. Of the more frequently occurring mutations, we observed that many arose repeatedly after resistance to multiple other drugs (Figure 6b), indicating that their typical emergence pattern is not limited to early resistance. Conversely, katG S315T stood out as a frequently occurring mutation (Figure 6a) with very few instances of other resistance arising prior to its gain (Figure 6b). Despite the high level of complexity in the stepwise acquisition of drug resistance mutations in M. tuberculosis (Figure 3), katG S315T was by far the most common mutation to evolve prior to MDR (Figure 2, Supplementary Notes E–G; Supplementary Tables 12–13). Of the 321 independent arisals of katG S315T in our dataset, 302 (94%) occurred at the earliest node where drug resistance was present.

Discussion

We constructed the largest global dataset of M. tuberculosis whole genomes analyzed to date, consisting of 5,310 diverse strains. Though the global distribution of strains in this dataset does not reflect the global incidence of TB for some regions (Supplementary Figure 1b), our unique dataset had a broad geographic distribution and deep sampling of drug-resistant strains, including MDR and XDR strains from multiple lineages and regions. We were, therefore, able to dissect the step-by-step evolution of drug resistance and identify harbinger resistance mutations that emerged prior to development of MDR-TB.

We observed that MDR- and XDR-TB evolved many independent times, in different lineages and regions of the world, suggesting that there are many “permissive” environments that have allowed MDR- and XDR-TB to emerge repeatedly. Molecular diagnostic tests for drug-resistant TB could be improved by incorporating knowledge of the global patterns of resistance emergence. We observed that the distribution of arisals of specific resistance mutations was fairly constant across the globe, indicating that drug resistance has arisen via common mechanisms worldwide. Thus, a universal diagnostic for detecting resistance to the eight drugs examined here may be achievable without need for regional specialization. Without phenotypic drug susceptibility data for all included strains, we were not able to identify novel drug resistance mutations or quantify the amount of drug resistance that remains unexplained by our curated list of polymorphisms; however, we expect it to be small[27].

By dissecting the step-by-step evolution of drug-resistant M. tuberculosis across the phylogeny, we observed that patterns in the order of emergence of drug resistance also appeared to be conserved globally. In particular, across all lineages and geographic regions, isoniazid resistance overwhelmingly arose before rifampicin resistance. Some regions of the world, such as Iran,[31] are reported to have a high incidence of rifampicin-resistance; however, our results suggest that rifampicin mono-resistance rarely leads to MDR-TB. Though the effects of convergent evolution among frequently evolving mutations could cause isoniazid resistance evolutions to be dated further back in time than when they actually arose, we have taken care to minimize such effects (Supplementary Note H). Supporting our results, this relative ordering of isoniazid and rifampicin resistance was consistent with prior findings based on genomic data in Russia[7], South Africa[6] and South America[7], as well as analysis of a large global collection of phenotypic data[8].

Why would isoniazid resistance arise first? We show that the earlier clinical introduction of isoniazid is not a major contributor to the earlier arisal of isoniazid resistance; our dating analysis indicates that isoniazid resistance arises before rifampicin resistance across all time periods, including recently (Supplementary Note D). However, there are many alternative, though not definitive, explanations for this preferential ordering. Isoniazid is a prodrug, which must first be activated by KatG (encoded by katG), the catalase-peroxidase[32] to form an adduct with NAD[33], which then inhibits InhA (encoded by inhA), an NADH-dependent enoyl-acyl carrier protein reductase[34] and, ultimately, inhibits mycolic acid biosynthesis[35]. The major mechanisms of isoniazid resistance include 1) mutations in katG, a non-essential gene, which result in failure to activate isoniazid, and 2) either upregulation or target-modification of InhA. Rifampicin inhibits the β-subunit of the mycobacterial RNA polymerase, encoded by a single, essential gene, rpoB[36]. M. tuberculosis grown in vitro have higher spontaneous mutation rates toward isoniazid resistance than rifampicin resistance[37,38], which could be due to the greater number of mutations that can lead to isoniazid resistance as compared to rifampicin resistance i.e., any inactivating mutation within katG can result in isoniazid resistance whereas only specific non-inactivating mutations in rpoB can result in rifampicin resistance. However, we observed that a single mutation in katG, S315T, accounted for the majority of isoniazid resistance arisals, and that, overall, there were ~20% more independent arisals of resistance to rifampicin than isoniazid indicating that the relative rates of resistance in vivo may differ from those calculated in vitro.

Another possible explanation for the ordering is that isoniazid resistant strains, including those carrying the prevalent katG S315T mutation, are more likely to develop resistance to other drugs. While previous in vitro studies have shown a difference in the types of rifampicin resistance mutations arising in isoniazid-resistant backgrounds[39], there is no evidence that isoniazid resistant strains are transformed into “hypermutators”[40]. Further, the sequence surrounding the katG S315 site does not appear susceptible to mutation nor does it appear to be a mutational hotspot in vitro[37]. However, as others and we have shown, this specific mutation is common among isoniazid resistant clinical isolates, indicating that it is well tolerated in vivo. This is likely due to the fact that this mutation preserves mycobacterial catalase activity while still preventing activation of isoniazid[41]. This preserved fitness may impact the evolutionary adaptive landscape[42,43] through which M. tuberculosis may acquire future resistance. Such a fitness landscape, which takes into account the relative fitness of different combinations of resistance and compensatory mutations, may produce a restricted set of evolutionary paths leading to MDR-TB.

A third possibility is that there is differential drug availability within the body, either due to pharmacokinetic effects[44,45] or differential clinical penetrations of the drugs into lesions[46], that may influence the ordering of emergence of mutations. Current treatment regimens that result in suboptimal dosing of rifampicin[47-49] may result in effective mono-exposure to isoniazid, increasing the likelihood of developing isoniazid resistance first. Isoniazid preventative therapy (IPT)[50], treatment of suspected cases of latent TB with only isoniazid, could provide an opportunity for isoniazid resistance to develop prior to exposure to other drugs. However, IPT is not commonly used in most of the countries for which we have assembled data and, therefore, is unlikely to have a major role in the early arisal of isoniazid resistance in our dataset.

Early identification and appropriate treatment of individuals with isoniazid mono-resistant strains, such as with non-isoniazid based regimens[51], may prevent the selection and eventual transmission of additional MDR strains. The worldwide case rate of isoniazid mono-resistance is estimated to be as high as 2%–15%[52-56], or 200,000–1.4 million cases per year. Several studies have shown that patients harboring isoniazid mono-resistant strains have worse clinical outcomes than those harboring susceptible strains[57-59], and enhanced treatment regimens for such strains resulted in lower rates of treatment failure and acquired drug resistance[59]. One large retrospective cohort study also points to early detection of isoniazid mono-resistance for improved outcomes[52]. In particular, the katG S315T mutation has been associated with unfavorable treatment outcome and increased relapse in one population[60]. These and other results challenge the predictions of an earlier mathematical modeling study, which prognosticated that incorporation of isoniazid resistance onto a molecular test in India would provide only a negligible benefit to the control of MDR-TB[61].

Our large dataset confirms that GeneXpert, currently the most widely used rapid molecular diagnostic for the diagnosis of M. tuberculosis and MDR-TB, performs excellently as a surrogate marker for MDR, irrespective of lineage and/or region of the globe. However, as it includes only rifampicin resistance mutations, GeneXpert does not identify drug resistance at the earliest available opportunity. Thus, diagnostic algorithms that rely upon application of GeneXpert alone will allow rifampicin-susceptible but otherwise drug-resistant strains to propagate unchecked. In fact, our evolutionary analysis revealed that, by the time a GeneXpert identifiable rifampicin resistance mutation is acquired, oftentimes multiple additional resistance-conferring mutations are already present. Additional commercially available diagnostics, such as the Hain MTBDRplus and Hain MTBDRsl ver 2.0 line probe assays, are available and detect a broader set of resistance mutations. Despite logistical considerations that make practical application of this technology more difficult—such as moderate turn-around times and need for specialized laboratory facilities—line probe assays are able to detect isoniazid resistance (HainMTBDRplus) with excellent specificity[62]. However, current diagnostic algorithms in certain TB endemic countries[63] call for the application of these tests only after rifampicin resistance has been identified by GeneXpert. Therefore, more comprehensive diagnostic tests are not being appropriately used to their full potential to identify rifampicin-susceptible but otherwise drug-resistant strains.

Through an evolutionary analysis of a diverse, global dataset encompassing 5,310 strains of M. tuberculosis, we observed that recent de novo emergence of MDR-TB in the last 10 years is a significant contributor to global MDR today. Thus, to stem the development of additional MDR strains, one should seek to identify resistant strains in the “pre-MDR” stage, during which there are additional therapeutic options and improved treatment outcomes. Identification of harbinger mutations, such as katG S315T, may serve as an early warning signal that MDR may soon develop. Focusing on common, early-occurring mutations could improve the design of diagnostic tests, aimed to target the earliest-occurring signatures of drug resistant bacteria. Future prospective research will be needed to determine whether these harbinger mutations increase the risk of MDR emergence in a given population. If substantiated, surveillance efforts for harbinger mutations may assist organizations to better allocate intensified TB control resources to at-risk areas, and to target drug resistance in the pre-MDR stage.

Online Methods

Whole Genome Sequencing and datasets

References detailing the sequencing methods for the published datasets can be found in Supplementary Table 1. Sequencing data for the TB-ARC projects (Supplementary Table 1) was generated at the Broad Institute as in Cohen et al.[6] Additional information about samples for each of these unpublished projects can be found at the Broad Institute website (http://olive.broadinstitute.org/projects/tb_arc). The goal of the TB-ARC project was to create a catalog of mutations for antibiotic resistance in M. tuberculosis to inform diagnostics. As such, strains from each of the countries represented both drug sensitive and drug resistant isolates that would enable curation of such a catalog. For the Indian dataset, 223 drug resistant and drug sensitive strains, representative of lineages found in India (particularly lineages 1 and 3) were selected and sequenced from studies in Tiruvallur and Madurai in Southern India. For the MRC dataset, 189 primarily drug-resistant strains from South African and Swaziland were selected for sequencing. For the CDRC dataset, 179 genomes from South Korea and Uganda, with a wide variety of drug resistance patterns, as well as extensive characterization of drug susceptibility profiles, were sequenced. For the Swedish dataset, 150 genomes were collected, primarily from Sweden’s immigrant population. This set includes a complete collection of all 141 MDR and XDR strains identified nationwide in Sweden in the period 2003–2013. For the dataset from Moldova, 95 genomes were selected from the countrywide specimen and data repositories. For the Romanian dataset, 34 genomes were sequenced, with the goal of describing drug resistant strains circulating in Romania and their diversity. For the Iranian dataset, 33 primarily highly drug resistant samples, including totally drug resistant (TDR) samples were sequenced.

The study protocol for these TB-ARC projects was approved by the Massachusetts Institute of Technology Committee on the Use of Humans as Experimental Subjects. Informed consent was obtained from all subjects, or else an appropriate waiver of consent was obtained.

For all of these TB-ARC projects, genomic DNA was extracted using published methods[66]. Library preparation and whole genome sequencing (WGS) were performed as previously described on the Illumina HiSeq 2000 at the Broad Institute [67]. Sequencing data were submitted to the Sequence Read Archive at NCBI under the umbrella BioProject identifiers listed in Supplementary Table 1.

Alignments

Raw read data for 8,136 strains were downloaded from the SRA (see SRA accessions in Supplementary Table 1). Reads were mapped onto a reference strain of H37Rv (GenBank accession number CP003248.2) using BWA version 0.7.10[68]. Variants were identified using Pilon version 1.11 as described [67]. The global M. tuberculosis lineage designations used in our analysis, as well as each strain’s spoligotype, were predicted using digital spoligotyping as in Cohen et al., 2015[6].

We eliminated 824 strains that did not pass our quality control filters: average sequencing depth of coverage >20X; fraction of long insertions <0.2; ambiguity rate <0.5 (to remove samples that were suspected to represent mixes of different genotypes); number of low coverage bases <250,000; and having a single match to one lineage in our lineage-prediction algorithm. We also eliminated strains for which Pilon failed three times. Of the remaining 7,312 samples, we removed 1,970 strains with no “country” metadata or description in a publication; 19 strains with substantial non-tuberculous mycobacteria contamination; as well as 13 additional duplicate patient samples. These filters resulted in a final set of 5,310 strains for analysis.

Emu[69] was run to canonicalize variants. We conducted phylogenetic analyses for the entire set of 5,310 strains, as well as for a subset corresponding to each lineage and each United Nations geographical subregion[23] with >30 strains (Supplementary Table 3). For each set, all sites with unambiguous single nucleotide polymorphisms (SNPs) in at least one strain were combined into a concatenated alignment. Ambiguous positions were treated as missing data. The concatenated alignment was then were used to generate a midpoint rooted phylogenetic tree using FastTree[70] version 2.1.8.

Drug resistance mutations

A curated list of genomic polymorphisms that confer drug resistance was defined for eight drugs: rifampicin, ethambutol, isoniazid, ethionamide, ofloxacin, pyrazinamide, streptomycin, and kanamycin. This was based on a literature review and consideration of current molecular drug resistance diagnostics. All mutations incorporated in current molecular diagnostics in standard practice were included. This included the Xpert MTB/RIF[5], the Hain Genotype MTBDRplus, the Hain Genotype MTBDRsl Line Probe Assay, as well as the Hain MTBDRsl ver 2.0 (see urls). Additional resistance mutations were selected for inclusion based on laboratory evidence and recent compelling genomic evidence that these mutations encode for resistance (Supplementary Tables 4–5). Because of greater uncertainty calling longer variants in our data, we excluded insertions and deletions longer than 10bp. Using this curated list, we identified 392 drug-resistance mutations among the 231,898 total variants observed in our full dataset across 5,310 strains.

Recent reports suggest that currently tabulated mutation sets account for the majority of phenotypic resistance[27,71,72]. While pyrazinamide was predicted at a lower sensitivity than other drugs in one recent analysis[71], including all loss-of-function mutations in pncA, as we have done in this study, would likely greatly improve sensitivity to detect pyrazinamide resistance.

Evolution of drug resistance mutations

We used PAUP[73] version 4.0b10 to reconstruct gains and losses of drug resistance mutations across the phylogenetic tree. We performed this analysis both for the full phylogeny of all 5,310 strains, as well as for individual phylogenies for each of the 11 geographic subregions and five lineages with >30 strains. PAUP was run using a cost matrix that assigned a 20x greater cost for a loss event relative to a gain event.

When examining the relative ordering of resistance mutations at two different nodes, we removed portions of the tree with potentially ambiguous topology. We removed node pairs from our analysis when the ancestral node had a bootstrap value <90%, as well as node pairs where the longest of the individual branch lengths separating them was >1e-4. 25% of the branches in our phylogeny had branch lengths >1e-4. Our combined branch length and bootstrap filtering removed a total of 48% of the node pairs.

Dating the phylogeny

BEAST[29] version 1.8.2 was used to estimate dates of acquisition of drug resistance mutations in the phylogenies of lineages 1–4. Lineages 1 and 3 contained a small enough number of strains to run BEAST (494 and 431 strains, respectively). However, because the size of lineages 2 and 4 were beyond the current capabilities of the BEAST algorithm, we subdivided these lineages into subsets of size 400–700 strains, and ran BEAST separately on each of these subsets. First, we removed very closely related strains from lineages 2 and 4. To do this, we clustered strains using simple agglomerative hierarchical clustering. For each cluster containing multiple sequences with <10 SNP differences in the core region aligning to H37RV, we kept only one strain. This reduced the number of unique strains in lineage 2 to 978, and the number of unique strains in lineage 4 to 1,556. We then manually examined the phylogenies to split the remaining lineage 2 strains into two clades (lineage 2a with 462 strains, and lineage 2b with 516 strains), and the remaining lineage 4 strains into three clades (lineage 4a with 423 strains, lineage 4b with 413 strains, and lineage 4c with 720 strains). We constructed a phylogenetic tree for each of these seven subsets using FastTree[70] version 2.1.8 (one tree from lineage 1, two trees from lineage 2, one tree from lineage 3, and three trees from lineage 4). We then ran BEAST to estimate dates of acquisition of drug resistance mutations in these seven clades representing lineages 1–4.

Because of the large size of our dataset and small spread in sample isolation dates, we used simplified parameters and a fixed evolutionary rate. We ran BEAST twice, using a fixed upper and lower bound for the evolutionary rate. Because of the wide range of rates previously observed in M. tuberculosis, including varying rates for different strains[6,11,19,38,74-76], we used a lower bound of 0.3 mutations/ genome/ year [74] and an upper bound of 0.6 mutations/ genome/ year[6], to cover the entire range of published rates across all lineages. Isolation dates for each sample were used as input to BEAST. If our metadata included only a range of isolation dates, we selected the midpoint of this date range. We enforced the topology of our input tree generated using FastTree[70]. We used the following parameters when running BEAST: GTR, empirical base frequencies, and no site heterogeneity model. BEAST was run for a minimum of 10 million iterations, sampling every 1,000 iterations. The program Tracer was used to examine mixing and effective sample size in order to assess chain length and model convergence. If the effective sample size (ESS) with 10 million iterations for all variables was not >100, then BEAST was run for additional iterations, until ESS values were all >100. The first 1 million iterations were excluded as “burn-in”. Estimated dates are given with 95% highest posterior density (HPD) intervals.

We also used the BEAST data to calculate a date for the arisal of each MDR node. In order to calculate the number of strains descending from each MDR node, we included the closely related strains that were removed after performing the hierarchical clustering.

Data Availability Statement

Newly sequenced TB-ARC data has been deposited in the SRA database. Accessions for all newly sequence data, as well as previously published data, are listed in Supplementary Table 1. Data are also available on the Broad Institute’s website (http://olive.broadinstitute.org/projects/tb_arc)