Compensatory evolution drives multidrug-resistant tuberculosis in Central Asia

Bacterial factors favoring the unprecedented multidrug-resistant tuberculosis (MDR-TB) epidemic in the former Soviet Union remain unclear. We utilized whole genome sequencing and Bayesian statistics to analyze the evolutionary history, temporal emergence of resistance and transmission networks of MDR Mycobacterium tuberculosis complex isolates from Karakalpakstan, Uzbekistan (2001–2006). One clade (termed Central Asian outbreak, CAO) dating back to 1974 (95% HPD 1969–1982) subsequently acquired resistance mediating mutations to eight anti-TB drugs. Introduction of standardized WHO-endorsed directly observed treatment, short-course in Karakalpakstan in 1998 likely selected for CAO-strains, comprising 75% of sampled MDR-TB isolates in 2005/2006. CAO-isolates were also identified in a published cohort from Russia (2008–2010). Similarly, the presence of mutations supposed to compensate bacterial fitness deficits was associated with transmission success and higher drug resistance rates. The genetic make-up of these MDR-strains threatens the success of both empirical and standardized MDR-TB therapies, including the newly WHO-endorsed short MDR-TB regimen in Uzbekistan.


Introduction
Multidrug-resistant tuberculosis (MDR-TB), caused by Mycobacterium tuberculosis complex (MTBC) strains that are resistant to the first-line drugs isoniazid and rifampicin, represent a threat to global TB control. Barely 20% of the estimated annual 480,000 new MDR-TB patients have access to adequate second-line treatment regimens. The majority of undiagnosed or ineffectively treated MDR-TB patients continue to transmit their infection and suffer high mortality (WHO, 2016).
Based on early observations that the acquisition of drug resistance could lead to reduced bacterial fitness (Middlebrook and Cohn, 1953), it was hypothesized that drug-resistant MTBC-strains had a reduced capacity to transmit, and would not widely disseminate in the general population (Borrell and Gagneux, 2009;Billington et al., 1999;Burgos et al., 2003;Dye and Espinal, 2001;Andersson and Levin, 1999). This optimistic scenario has been invalidated by the now abundant evidence for transmission of MDR and extensively drug-resistant MTBC-strains (XDR-TB; MDR-TB additionally resistant to at least one fluoroquinolone and one injectable aminoglycoside) in healthcare and community settings (Borrell and Gagneux, 2009;Gagneux et al., 2006;Müller et al., 2013;Pym et al., 2002;Comas et al., 2012). In former Soviet Union countries, which experience the highest MDR-TB rates worldwide, the expansion of drug-resistant MTBC-clones is thought to be promoted by interrupted drug supplies, inadequate implementation of regimens, lack of infection control and erratic treatment in prison settings (Balabanova et al., 2004;Casali et al., 2014a). Continued transmission is thought to be aided by the co-selection of mutations in the bacterial population that compensate for a fitness cost (e.g. growth deficit) associated particularly with the acquisition of rifampicin resistance mediating mutations (Borrell and Gagneux, 2009;Andersson and Levin, 1999;Gagneux et al., 2006;Müller et al., 2013;Pym et al., 2002;Comas et al., 2012). The compensatory mechanism for rifampicin-resistant MTBC-strains is proposed to be associated with structural changes in the RNA-polymerase subunits RpoA, RpoB, and RpoC that increase transcriptional activity and as a consequence enhance the growth rate (Comas et al., 2012). However, the impact of these bacterial genetic factors on the epidemiological success of MDR-MTBC strains and implications for current and upcoming MDR-TB treatment strategies remain unexplored.
We utilized whole-genome sequencing (WGS) to retrace the longitudinal transmission and evolution of MTBC-strains toward MDR/pre-XDR/XDR geno-and phenotypes in Karakalpakstan, Uzbekistan. In this high MDR-TB incidence setting, the proportion of MDR-TB among new TB-patients increased from 13% in 2001 to 23% in 2014 despite the local introduction of the World Health Organization (WHO) recommended DOTS strategy in 1998 and an initially limited MDR-TB treatment program in 2003 (Cox et al., 2007;Ulmasova et al., 2013). We expanded our analyses by including a WGS dataset of MDR-MTBC isolates from Samara, Russia (2008Russia ( -2010 (Casali et al., 2014a) to investigate clonal relatedness, resistance and compensatory evolution in both settings.

Study population and MTBC phenotypic resistance (Karakalpakstan, Uzbekistan)
Despite differences in sampling for cohort 1 (cross-sectional, 2001-2002) and cohort 2 (consecutive enrollment of MDR-TB patients, 2003 (see Materials and methods), patients showed similar age, sex distributions, and proportion of residence in Nukus, the main city in Karakalpakstan (Uzbekistan) (Appendix-table 1). While the majority of strains from both cohorts were phenotypically resistant to additional first-line TB drugs (i.e. beyond rifampicin and isoniazid), combined resistance to all five first-line drugs was significantly greater in cohort 2 (47% in cohort 2 compared to 14% in cohort 1, p<0.0001). The same was true for resistance to the second-line injectable drug capreomycin (23% in cohort 2 compared to 2% in cohort 1, p=0.0001) (Appendix-table 1). This finding was surprising as the isolates from cohort two patients -who were treated with individualized second-line regimens predominately comprising ofloxacin as the fluoroquinolone and capreomycin as the second-line injectable -were all obtained before the initiation of their treatment. In addition, there was no formal MDR-TB treatment program in Karakalpakstan prior to 2003. These elements imply that the higher rate of resistance to capreomycin was attributable to infection by already resistant strains (i.e. to primary resistance).

MTBC population structure and transmission rates
Utilizing WGS, we determined 6979 single-nucleotide polymorphisms (SNPs) plus 537 variants located in 28 genes and upstream regions associated with drug resistance and bacterial fitness (Supplementary file 1). The corresponding phylogeny revealed a dominant clade comprising 173/ 277 (62.5%) closely related isolates within MTBC lineage 2 (particularly Beijing-genotype) ( Figure 1). This group, termed Central Asian Outbreak (CAO), showed a highly restricted genetic diversity (median pairwise distance of 21 SNPs, IQR 13-25) and was differentiated from a set of more diverse isolates by 38 specific SNPs (Appendix-figure 1, Supplementary file 1 ). The proportion of CAOisolates was similar between 2001CAOisolates was similar between -2002CAOisolates was similar between and 2003CAOisolates was similar between -2004, but increased to 76% in 2005-06 (p<0.01). Over the same time periods, the proportions of other groups remained stable or decreased (Appendix- figure 2).
We then sized transmission networks (measured by transmission indexes, see Materials and methods) supposed to reflect human-to-human transmission over the last~10 years based on a maximum of 10 differentiating SNPs between two isolates. Transmission rates varied, even among closely related outbreak isolates ( Figure 1). Beijing-CAO-isolates formed particularly large transmission networks (>50 patients; Figure 1); 96.0% (166/173) of all Beijing-CAO isolates were associated with recent transmission (i.e. transmission index !1), versus 48.4% (31/64) of non-eLife digest Multidrug-resistant tuberculosis, often shortened to MDR-TB, is a public health crisis with close to half a million patients falling ill each year globally. Some strains of the bacterium Mycobacterium tuberculosis, which causes tuberculosis disease, are resistant to the two most effective drugs used to treat the infection. As a result, patients with MDR-TB require a longer treatment of up to two years, often with severe side effects and a low chance of cure. Resistant strains of the bacteria are usually weaker than drug-susceptible strains. So, for a long time, large MDR-TB epidemics were considered to be unlikely and outbreaks of MDR-TB were often regarded as locally contained phenomenona.
Recent research has shown that MDR-TB strains are often just as likely as drug-susceptible strains to be transmitted and therefore just as likely to cause large country-wide outbreaks. It has also become clear that the resistant bacteria acquire additional mutations over time to compensate for any weakness. However, a lack of detailed history of outbreaks has meant the role of the genetics of MDR-TB bacteria has not been fully understood. Without this knowledge, prevention of future outbreaks and containment of the most successful strains in areas with a high burden of disease is difficult.
To address this, Merker, Barbier et al. reconstructed the evolutionary history of MDR-TB strains obtained in 2001-2006 from an outbreak in Uzbekistan. Whole genome sequencing followed by statistical analysis highlighted one predomininant strain that likely emerged in the mid-1970s, when the country was part of the former Soviet Union. This strain has since acquired mutations that make it resistant to eight different drugs. The most successful bacterial strains found also had compensatory mutations that seem to aid their survival.
In 1998, the health authorities implemented a TB treatment program in the region without knowing the true extent of the MDR-TB outbreak at that time. Testing for drug resistance was not routinely available, and Merker, Barbier et al. saw that MDR-TB strains resistant to the drugs used spread in the study region and were later also found independently in Russia.
A lack of routine testing for drug resistance in TB remains common in many countries with high burdens of the disease. These findings emphasize the need for universal access to tests for TB drug resistance, therapies tailored for individual patients, and access to new and repurposed drugs to reduce the risk of future outbreaks of drug-resistant TB. Columns show drug resistance associated mutations to first-and second-line drugs (different mutations represented by different colors), genetic classification of pre-XDR (purple) and XDR (pink) isolates, and putative compensatory mutations in the RNA polymerase genes rpoA, rpoB and rpoC. Transmission index represents number of isolates within a maximum range of 10 SNPs at whole genome level. MTBC lineage two isolates (Beijing genotype) are differentiated into three clades (i.e. Central Asian Outbreak (CAO), group 2 and group 3). Isolates belonging to lineage 4 (Euro-American) are colored in grey: H = isoniazid, R = rifampicin, S = streptomycin, E = ethambutol, Z = pyrazinamide, FQ = fluoroquinolone, AG = aminoglycosides, Km = kanamycin Cm = capreomycin, TA = thioamide, PAS = para aminosalicylic acid. DOI: https://doi.org/10.7554/eLife.38200.003 CAO Beijing isolates (p<0.0001) and 57.5% (23/40) of non-Beijing isolates (p<0.0001) (Supplementary file 1). In addition, the large CAO transmission network exhibited higher levels of drug resistance relative to non-Beijing strains, as reflected by the larger number of drugs for which phenotypic (p=0.0079) and genotypic drug resistance (p=0.0048) was detected Appendixfigure 3).

Evolutionary history of CAO strains in Karakalpakstan
In order to gain more detailed insights into the emergence of resistance mutations in the evolutionary history of the CAO clade, we sought to employ a Bayesian phylogenetic analysis for a temporal calibration of the CAO phylogeny and an estimation of the mutation rate. Using an extended collection of more diverse CAO isolates (n = 220) from different settings (see Materials and methods), we initially compensated for the restricted sampling time frame of the Karakalpakstan dataset (2001)(2002)(2003)(2004)(2005)(2006). A linear regression analysis showed correlation between sampling year and root-to-tip distance and even a moderate temporal signal (p=0.00039, R 2 = 5.2%, Appendix-figure 4), allowed for a further estimation of CAO mutation rates and evaluation of molecular clock models using Bayesian statistics as discussed previously (Duchêne et al., 2016). Based on the marginal L estimates collected by path sampling, we found a strict molecular clock with tip dates to be a reasonable  . Compensatory mutations and drug resistance levels. Comparisons between isolates carrying compensatory mutations (in orange) and isolates with no-compensatory mutations (in blue), from the Karakalpakstan dataset. (A) Boxplot showing number of resistance mutations for the two categories (without or with compensatory mutations). The two categories were significantly different (two-sample t-test p=1.2Â10 À10 ). (B) Bubble plots showing the transmission index (number of isolates differing by less than 10 SNPs) as a function of antibiotic resistance related mutations. Bubble sizes are proportional to the number of isolates. (C) Density plot of the number of resistance-conferring mutations for four groups of isolates sourced from the Karakalpakstan data. Proportions are adjusted by using Gaussian smoothing kernels. The four groups are composed of non-CAO isolates with no compensatory mutations; non-CAO isolates carrying compensatory mutations; CAO isolates with no compensatory mutations and CAO isolates carrying compensatory mutations. These groups are respectively colored in light blue, dark blue, light orange and light red. DOI: https://doi.org/10.7554/eLife.38200.005 model for CAO isolates (Appendix-table 2). Mutation rate estimates (under a relaxed clock model) ranged on average from 0.88 to 0.96 Â 10 À7 substitutions per site per year (s/s/y), depending on the demographic model, in favor for the Bayesian skyline model with a mutation rate of 0.94 Â 10 À7 (s/s/y) (95% HPD 0.72-1.15 Â 10 À7 (s/s/y)) (Appendix-table 2). Comparing different demographic models for the CAO-Karakalpakstan dataset (n = 173) an exponential growth model and a Bayesian skyline model were superior over the constant size demographic prior.
Employing the Bayesian skyline model with a strict molecular clock set to 0.94 Â 10 À7 (s/s/y) specifically we determined that the most recent common ancestor (MRCA) of the CAO-clade emerged around 1974 (95% highest posterior density (HPD) 1969(HPD) -1982. The time to the MRCA was confirmed with the exponential growth demographic model (1977, 95% HPD 1977-1982. The MRCA already exhibited a streptomycin resistance mutation (rpsL K43R) (Figure 2), and acquired isoniazid resistance (katG S315T) in 1977(katG S315T) in (95% HPD 1973(katG S315T) in -1983. The CAO-population size then rose contemporaneously with multiple events of rifampicin, ethambutol, ethionamide, and para-aminosalicylic acid resistance acquisition in different branches ( Figure 2). As an illustration, the most frequent CAO-clone (upper clade in Figure 2) acquired ethambutol and ethionamide resistance mutations (embB M306V, ethA T314I) around 1984 (95% HPD 1982(95% HPD -1989, and an MDR-genotype (rpoB S450L) around 1986 (95% HPD 1985(95% HPD -1992. The effective population size reached a plateau before fixation of mutations in the ribD promoter region (leading to para-aminosalicylic acid resistance) and rpoC N698S, putatively enhancing its fitness around 1990 (95% HPD 1989-1994) ( Figure 2). Independent fixation of pyrazinamide (pncA Q10P and I133T) and kanamycin (eis À12 g/ a) resistance-associated mutations was detected in 1992 and 1991 (both with 95% HPD rounded to [1991][1992][1993][1994][1995][1996] (Figure 2). Shown are the coefficients and 95% confidence bands of multiple linear regression of the transmission success score, defined as the size of clusters diverging by at most N SNPs and divided by N or, equivalently, the size of clusters that evolved over N years divided by N. The presence of compensatory mutations was independently associated with transmission success, with a maximum association strength found for SNP distances ranging from 10 to 20 SNPs, corresponding to transmission clusters beginning around 1995. DOI: https://doi.org/10.7554/eLife.38200.006 To further account for uncertainties of substitution rates and thus fixation of drug resistance within the CAO-clade we ran the best models (Bayesian skyline and exponential growth) with the upper and lower HPD interval of the best clock estimate (see above). Similarly, the most recent fixation of the putative compensatory mutation rpoC N698S was 1994N698S was (95% HPD 1992N698S was -1996, still years before implementation of the systematic DOTS-program in Karakalpakstan in 1998. Interestingly, the DOTS implementation coincided with a second effective population size increase (Figure 2). At that time, distinct CAO-clades already exhibited pre-XDR (in this context MDR plus kanamycin resistance) resistance profiles, mediating resistance to as many as eight different anti-TB drugs. Of note, only a single isolate was identified as harboring a gyrA mutation (A90V), associated with fluoroquinolone resistance (Supplementary file 1). At the end of the study period in 2006, we observed a pre-XDR rate among CAO isolates of 52.0% (90/173), compared to 35.9% (23/64) among other Beijing isolates (p=0.03) and compared to 42.5% (17/40) among non-Beijing isolates (p=0.30) (Supplementary file 1).
Regression-based analyses of transmission success scores in the Beijing-CAO clade confirmed that the presence of compensatory mutations was strongly associated with cluster sizes independent of the accumulation of resistance mutations ( Figure 4). This pattern was mostly observed for clusters initiated in the late 1980s and the 1990s.

Combined analysis of MDR-TB cohorts from Karakalpakstan and Samara (2001-2010)
To place our analyses in a broader phylogenetic and geographic context, we combined our Karakalpakstan genome set with previously published genomes of 428 MDR-MTBC isolates from Samara (Casali et al., 2014a), a Russian region located~1700 km from Nukus, Karakalpakstan. This analysis showed that Beijing-CAO isolates accounted for the third largest clade in Samara (Casali et al., 2014a). Conversely, the second largest clade in Samara, termed Beijing clade B according to Casali et al (Casali et al., 2014a;Casali et al., 2012), or European/Russian W148 (Merker et al., 2015), was represented in Karakalpakstan by a minor clade ( Figure 5). Considering a third Beijing clade (termed clade A) restricted to Samara (Casali et al., 2014a), three major Beijing outbreak clades accounted for 69.6% (491/705) of the MDR-TB cases in both regions.
The three Beijing clades (A, B, and CAO) in Samara and Karakalpakstan had more drug resistance conferring mutations (in addition to isoniazid and rifampicin resistance) with means of 5.0 (SEM 0.07), 4.2 (SEM 0.18), and 4.7 (SEM 0.11), respectively (Appendix- figure 5), than compared to only 3.6 (SEM 0.20) additional genotypic drug resistances (p<0.0001, p=0.0143, p<0.0001) for other Beijing isolates in both settings. Isolates belonging to other MTBC genotypes (mainly lineage four clades) were found with a mean of 2.6 (SEM 0.20) additional drug resistance mediating mutations, lower than any Beijing-associated group (p 0.0009) (Appendix- figure 5).
The impact of resistance conferring and compensatory mutations on the transmission success score in Beijing-A clade from Samara (Appendix- figure 7) was strikingly similar to the one observed in CAO isolates from Karakalpakstan. The presence of compensatory mutations, but not the accumulation of resistance mutations, was significantly and independently associated with network size in clusters originating in the 1980s and 1990s, with a maximum influence found in clusters starting in the late 1990s.
Critically, the high proportions of isolates detected in both settings with pre-XDR and XDR resistance profiles among the three major Beijing clades (clade A, 96%; clade B, 62%; clade CAO, 50%; Appendix-table 4, Figure 6) reveal the low proportion of patients that are or would be eligible to receive the newly WHO endorsed short MDR-TB regimen. As per definition of the WHO exclusion criteria, for example any confirmed or suspected resistance to one drug (except isoniazid) in the short regimen, only 0.5% (1/191 in Karakalpakstan) and 2.7% (8/300 in Samara) of the patients infected with either a Beijing clade A, B or CAO strain would benefit from a shortened MDR-TB therapy (Supplementary file 1).

Figure 5 continued
The presence of any putative compensatory mutation in the RNA polymerase genes rpoA, rpoB, rpoC is depicted with green bars and country of origin and a genotypic pre-XDR and XDR isolate classification is color coded. PAS = para aminosalycylic acid. DOI: https://doi.org/10.7554/eLife.38200.007 Figure 6. Percentage of drug resistance among 705 MDR-MTBC isolates from Samara (Russia) and Karakalpakstan (Uzbekistan). MDR-MTBC isolates stratified to three Beijing clades, other Beijing isolates and non-Beijing isolates. Proportions of isolates with identified molecular drug resistance mutations (see Supplementary file 1) which mediate resistance to multiple first-and second-line anti-TB drugs. Values are rounded. Drugs used in the WHO endorsed standardized short MDR-TB regimen marked with grey boxes. *The short MDR-TB regimen further includes high-dose isoniazid treatment, and clofazimine. In that regard, we identified 622/705 (85.4%) of the MDR-MTBC isolates with the well-known high-level isoniazid resistance mediating mutation katG S315T (Supplementary file 1), for clofazimine resistance mediating mutations are not well described. DOI: https://doi.org/10.7554/eLife.38200.008

Discussion
Using WGS combined with Bayesian and phylogenetic analyses, we reveal the evolutionary history and recent clonal expansion of the dominatant MDR/pre-XDR MTBC-clade in Karakalpakstan, Uzbekistan, termed the Central Asian outbreak (CAO). Strikingly, CAO-isolates were also found also in Samara, Russia, and vice versa isolates belonging to the second largest clade in Samara (Beijing clade B, i.e. European/Russian W148 (Casali et al., 2014a;Merker et al., 2015) were identified in Karakalpakstan, suggesting that the MDR-TB epidemic in this world region is driven by few outbreak clades. During the three last decades, these strains gradually accumulated resistance to multiple anti-TB drugs that largely escaped phenotypic and molecular diagnostics, and reduced treatment options to a restricted set of drugs that often cause severe side effects. In addition, our results suggest that compensatory mutations (in RNA-polymerase subunit coding genes) that are proposed to ameliorate growth deficits in rifampicin resistant strains in vitro are also crucial in a global epidemiological context allowing MDR and pre-XDR strains to form and maintain large transmission networks. The predominance of these strain networks, seen in two distant geographic regions of the former Soviet Union clearly limit the use of standardized MDR-TB therapies, for example the newly WHO endorsed short MDR-TB regimen, in these settings.
Temporal reconstruction of the resistance mutation acquisition and of changes in bacterial population sizes over three decades demonstrates that MDR outbreak strains already became resistant to both first-and second-line drugs in the 1980s. Fully first-line resistant strains massively expanded in the 1990s, a period that shortly preceded or immediately followed the end of the Soviet Union, years before the implementation of DOTS and programmatic second-line MDR-TB treatment. This is in line with the known rise in TB incidence that accompanied the economic breakdown in Russia during the 1990s (Institute of Medicine Forum on Drug Discovery, Development, and Translation and Russian Academy of Medical Science, 2011).
From a bacterial genetic point of view, our data show that particular MDR and pre-XDR clades are highly transmissible despite accumulation of multiple resistance mutations. The acquisition of compensatory mutations after introduction of low fitness cost resistance mutations (e.g. katG S315T (Pym et al., 2002), rpoB S450L (Gagneux et al., 2006), rpsL K43R (Bö ttger et al., 1998) seems the critical stage allowing for higher transmission rates. Multiple regression analyses further strengthened this hypothesis by demonstrating that the presence of fitness compensating variants was positively associated with transmission success in different settings and outbreak clades, independently of the accumulation of resistance mutations. Compensatory evolution thus appears to play a central role in driving large MDR-TB epidemics such as that seen with the Beijing CAO-clade.
A particular concern is the high prevalence of mutations conferring resistance to second-line drugs currently included in treatment regimens, among the dominant MDR-MTBC strains. Their detected emergence in a period preceding DOTS implementation, for example in Karakalpakstan, can be explained by past, largely empirical treatment decisions or self-medication. For instance, high frequencies of mutations in the ribD promoter region, and folC among Beijing-CAO isolates, associated with para-aminosalicylic acid resistance (Zheng et al., 2013;Zhao et al., 2014), are a likely consequence of the use of para-aminosalicylic acid in failing treatment regimens in the late 1970s to the early 1980s in the Soviet Union (USSR Ministry of Health, 1976;USSR Ministry of Health, 1983;Mishin, 2008). Likewise, the frequent independent emergence of mutations in the eis promoter and of rare variants in the upstream region of whiB7, both linked to resistance to aminoglycosides (mainly streptomycin and kanamycin) (Zaunbrecher et al., 2009;Reeves et al., 2013), probably reflects self-administration of kanamycin that was available in local pharmacies. Of note, prominent mutations such as katG S315T or rpoB S450L might have occurred multiple times independently in a bacterial population and inferring the common ancestor could lead to an overestimate of the TMRCA. However, this is not the case for rare and more diverse mutations, for example conferring resistance to pyrazinamide, PAS or kanamycin, thus further strengthening the historic fixation mentioned above.
The pre-existence of fully first-line resistant strain populations (e.g. CAO-Beijing in Karakalpakstan) likely contributed to the poor treatment outcomes observed among MDR-TB patients following the implementation of first-line DOTS treatment in 1998 (Cox et al., 2006). This period coincides with a detected CAO population size increase, likely reflecting the absence of drug susceptibility testing and therefore appropriate second-line treatment during extended hospitalization at the time, resulting in prolonged infectiousness of TB-patients and further spread of these strains.
The frequencies of fluoroquinolone resistance, mediated by gyrA and gyrB mutations, remained low among the Karakalpakstan MDR-MTBC isolates, which is consistent with the notion that such drugs were rarely used for treating TB in former Soviet Union countries (see Discussion (Casali et al., 2014a;USSR Ministry of Health, 1976;USSR Ministry of Health, 1983;Mishin, 2008). This observation explains the generally favorable MDR-TB treatment outcomes observed with the use of individualized second-line regimens, including a fluoroquinolone, in the latter MDR-TB treatment program in the Karakalpakstan patient population (Cox et al., 2007;Lalor et al., 2011). However, fluoroquinolone resistance, representing the last step towards XDR-TB, is already emerging as reported for strains in Beijing clade A and B (Casali et al., 2014a).
In conclusion, the (pre-) existence and wide geographic dissemination of highly resistant and highly transmissible strain populations most likely contributes to increasing M/XDR-TB incidence rates despite scaling up of the MDR-TB programs in some Eastern European and Russian regions (Ulmasova et al., 2013; Institute of Medicine Forum on Drug Discovery, Development, and Translation and Russian Academy of Medical Science, 2011; Medecins Sans Frontiere, 2013). Importantly, from the large spectrum of resistance detected among dominating strains in this study, it can be predicted that standardized therapies, including the newly WHO endorsed short MDR-TB regimen in Uzbekistan, are/will be largely ineffective for many patients in Samara and Karakalpakstan, and likely elsewhere in Eurasia. In order to successfully control the worldwide MDR-TB epidemics, universal access to rapid and comprehensive drug susceptibility testing, best supported by more advanced technologies, will be crucial for guiding individualized treatment with existing and new/ repurposed TB drugs and to maximize chances of cure and prevention of further resistance acquisition.

Study population, Karakalpakstan (Uzbekistan)
A total of 277 MDR-MTBC isolates derived from two separate cohorts were sequenced. The first cohort comprised 86% (49/57) of MDR-MTBC isolates from a cross-sectional drug resistance survey conducted in four districts in Karakalpakstan, Uzbekistan between 2001-2002 (Cox et al., 2006). An additional 228 isolates were obtained from TB-patients enrolled for second-line treatment in the MDR-TB treatment program from 2003 to 2006. These isolates represented 76% (228/300) of all MDR-TB cases diagnosed over the period. While the MDR-TB treatment program covered two of the four districts included in the initial drug resistance survey, the majority of isolates from both cohorts, 69% and 64% respectively, were obtained from patients residing in the same main city of Nukus (Appendix-table 1).

Study population, Samara (Russia)
To set the MDR-MTBC isolates from Karakalpakstan into a broader geographical perspective, raw WGS data of 428 MDR-MTBC isolates from a published cross-sectional prospective study in Samara, Russia from 2008 to 2010 (Casali et al., 2014a) were processed as described below and included into a composite MDR-MTBC dataset.

Whole genome sequencing
WGS was performed with Illumina Technology (MiSeq and HiSeq 2500) using Nextera XT library preparation kits as instructed by the manufacturer (Illumina, San Diego, CA). Fastq files (raw sequencing data) were submitted to the European nucleotide archive (see Supplementary file 1 for accession numbers). Obtained reads were mapped to the M. tuberculosis H37Rv reference genome (GenBank ID: NC_000962.3) with BWA . Alignments were refined with GATK (McKenna et al., 2010) and Samtools ) toolkits with regard to base quality re-calibration and alignment corrections for possible PCR artefact. We considered variants that were covered by a minimum of four reads in both forward and reverse orientation, four reads calling the allele with at least a phred score of 20, and 75% allele frequency. In the combined datasets, we allowed a maximum of 5% of all samples to fail the above-mentioned threshold criteria in individual genome positions to compensate for coverage fluctuations in certain genome regions; in these cases, the majority allele was considered. Regions annotated as 'repetitive' elements (e.g. PPE and PE-PGRS gene families), insertions and deletions (InDels), and consecutive variants in a 12 bp window (putative artefacts flanking InDels) were excluded. Additionally, 28 genes associated with drug resistance and bacterial fitness (see Supplementary file 1) were excluded for a conservative and robust phylogenetic reconstructions. The remaining single-nucleotide polymorphisms (SNPs) were considered as valid and used for concatenated SNP alignments. Further detailed methods of the phylogenetic reconstruction, molecular resistance prediction, strain-to-strain genetic distance, and Bayesian models are given as Appendix 1.

Transmission index
Based on the distance matrix (SNP distances), we further determined for every isolate the number of isolates that were in a range of 10 SNPs or less (in the following referred to as 'transmission index'). This 10 SNP-threshold was used to infer the number of recently linked cases, as considered within a 10-year time period, based on previous convergent estimates of MTBC genome evolution rate of » 0.5 SNPs/genome/year in inter-human transmission chains and in macaque infection models (Ford et al., 2011;Walker et al., 2013;Roetzer et al., 2013;Walker et al., 2014). This can include direct transmission events among the study population but also cases which are connected by a more distant contact which was not sampled. In the latter case, we assumed that two isolates with a maximum distance of 10 SNPs share a hypothetical common ancestor that is 5 SNPs apart from the two sampled isolates (considering a bifurcating phylogeny) and thus covers a timeframe of 5 SNPs over 0.5 SNPs/year equals 10 years between the two actual samples and a shared recent ancestor node/case (see also Appendix 1).

Genotypic drug resistance prediction
Mutations (small deletions and SNPs) in 34 resistance-associated target regions (comprising 28 genes) were considered for a molecular resistance prediction to 13 first-and second-line drugs (Supplementary file 1). Mutations in genes coding for the RNA-Polymerase subunits rpoA, rpoB (excluding resistance mediating mutations in the rifampicin resistance determining region (RRDR), and in codons 170, 400, 491), and rpoC were reported as putative fitness compensating (e.g. in vitro growth enhancing) variants for rifampicin-resistant strains as suggested previously (Comas et al., 2011;de Vos et al., 2013;Casali et al., 2014b;Cohen et al., 2015). A detailed overview of all mutations considered as genotypic resistance markers is given in Supplementary file 1. Mutations that were not clearly linked to phenotypic drug resistance were reported as genotypic non wild type and were not considered as genotypic resistance markers. When no mutation (or synonymous, silent mutations) was detected in any of the defined drug relevant target regions the isolate was considered to be phenotypically susceptible.

Phylogenetic inference (maximum likelihood)
We used jModelTest v2.1 and Akaike and Bayesian Information Criterion (AIC and BIC) to find an appropriate substitution model for phylogenetic reconstructions based on the concatenated sequence alignments (Appendix-table 5). Maximum likelihood trees were calculated with FastTree 2.1.9 (double precision for short branch lengths) (Price et al., 2010) using a general time reversible (GTR) nucleotide substitution model (best model according to AIC and second best model according to BIC), 1000 resamplings and Gamma20 likelihood optimization to account for evolutionary rate heterogeneity among sites. The consensus tree was rooted with the 'midpoint root' option in Fig-Tree (resulting in the expected tree topology of lineage 2-4 strains) and nodes were arranged in increasing order. Variants considered as drug resistance markers (see above) and putative compensatory variants were analyzed individually and mapped on the phylogenetic tree to define resistance patterns of identified phylogenetic clades.

Molecular clock model
In order to compute a time scaled phylogeny and employ the Bayesian skyline model (see below) for the identified Central Asian outbreak (CAO) clade, we sought to define an appropriate molecular clock model (strict versus relaxed clock) and a mutation rate estimate. Due to the restricted sampling timeframe of the Karakalpakstan dataset (2001)(2002)(2003)(2004)(2005)(2006), we extended the dataset for the model selection process with CAO isolates from Samara (2008)(2009)(2010) and 'historical' CAO isolates from MDR-TB patients in Germany (1995)(1996)(1997)(1998)(1999)(2000) thus allowing for a more confident mutation rate estimate. The strength of the temporal signal in the combined dataset, assessed by the correlation of sampling year and root-to-tip distance, was investigated with TempEst v1.5 (44). Regression analysis was based on residual mean squares, using a rooted ML tree (PhyML, GTR substitution model, 100 bootstraps), R-square and adjusted p-value are reported. For the comparison of different Bayesian phylogenetic models, we used path sampling with an alpha of 0.3, 50% burn-in and 15 million iterations (resulting in mean ESS values > 100), marginal likelihood estimates were calculated with BEAST v2.4.2 (45), and D marginal L estimates are reported relative to the best model.
First, we employed a strict molecular clock fixed to 1 Â 10 À7 substitutions per site per year as reported previously (Ford et al., 2011;Walker et al., 2013;Roetzer et al., 2013) without tip dating, a strict molecular clock with tip dating and a relaxed molecular clock with tip dating. BEAST templates were created with BEAUti v2 applying a coalescent constant size demographic model, a GTR nucleotide substitution model, a chain length of 300 million (10% burn-in) and sampling of 5000 traces/trees. Second, we ran different demographic models (i.e. coalescent constant size, exponential, and Bayesian skyline) under a relaxed molecular clock using tip dates and the same parameters for the site model and Markov-Chain-Monte-Carlo (MCMC) as described above.
Third, we tested and compared the best models for the Karakalpakstan CAO-clade under a strict molecular clock prior including the upper and lower 95% HPD interval (Appendix-table 2).
Inspection of BEAST log files with Tracer v1.6 showed an adequate mixing of the Markov chains and all parameters were observed with an effective sample size (ESS) >200 for the combined dataset (n = 220) and in the thousands for the Karakalpakstan CAO clade (n = 173), suggesting an adequate number of effectively independent draws from the posterior sample and thus sufficient statistical support. Other priors between the model comparisons were not changed.

Bayesian skyline plot
Changes of the effective population size of the CAO clade in Karakalpakstan over the last four decades were calculated with a Bayesian skyline plot using BEAST v2.4.2 (45) using a tip date approach with a strict molecular clock model of 0.94 Â 10 À7 substitutions per site per year (best model according to path sampling results, see above), and a GTR nucleotide substitution model. We further used a random starting tree, a chain length of 300 million (10% burn-in) and collected 5000 traces/trees. Again adequate mixing of the Markov chains and ESS values in the hundreds were observed. A maximum clade credibility genealogy was calculated with TreeAnnotator v2.

Impact of resistance-conferring and compensatory mutations on transmission success
We used multiple linear regression to examine the respective contributions of antimicrobial resistance and putative fitness cost-compensating mutations to the transmission success of tuberculosis. To take transmission duration into account, we computed, for each isolate and each period length T in years (from 1 to 40y before sampling), a transmission success score defined as the number of isolates distant of less than T SNPs, divided by T. This approach relied on the following rationale: based on MTBC evolution rate of 0.5 mutation per genome per year, the relation between evolution time and SNP divergence is such that a cluster with at most N SNPs of difference is expected to have evolved for approximately N years. Thus, transmission success score over T years could be interpreted as the size of the transmission network divided by its evolution time, hence as the average yearly increase of the network size. For each period T, the transmission success score was regressed on the number of resistance mutations and on the presence of putative compensatory mutations. The regression coefficients with 95% confidence intervals were computed and plotted against T to identify maxima, that is, time periods when the transmission success was maximally influenced by either resistance-conferring or -compensating mutations. These analyses were conducted independently on outbreak isolates of the Beijing-CAO clade in the Karakalpakstan cohort and of the Beijing-A clade in the Samara cohort.

Statistical analyses
Differences between cohorts and numbers of sampled isolates per year category were performed using Chi-squared analysis (mid-P exact) or Fisher's exact test, while comparison of median age was performed using the Mann-Whitney test. p-Values for pairwise comparisons of groups regarding pairwise genetic distances, number of resistant DST results and number of resistance related mutations were calculated with an unpaired t-test (Welch correction) or a t-test according to the result of the variances comparison using a F-test. Boxplot, bubble plots and density plots have been performed in R.

Funder Author
Leibniz Science Campus Evolutionary Medicine of the Lung

Matthias Merker Stefan Niemann
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Transmission index
In the context of this manuscript, we determined for every isolate the number of isolates that were in a range of 10 SNPs or less (in the following referred to as 'transmission index', see figure below).
The rationale to implement a 'transmission index' was the need to link each isolate with a continuous parameter (for further analysis) which reflects the number of recently linked cases, that is the extend of a putative transmission network. These networks are better reflected with a minimum spanning tree (rather than bifurcating phylogenies) which allows the visualization of super spreaders for instance as central nodes (see figure below). Thus an isolate with a high transmission index might well be linked to a patient that infected multiple secondary cases. The benefit compared to a categorical parameter like 'clustered' and'not clustered' is that the transmission index has the potential to indicate transmission hotspots within an outbreak scenario and is independent from a phylogenetic clade definition, which in turn would be difficult to assign due to the close genetic relationship in MTBC outbreaks and low bootstrap values for small sub-groups at the tips of a tree. The figure below illustrates the transmission index calculation per isolate. The upper left isolate for instance just has one sampled isolate within 10 SNP distance. A direct transmission event is relatively unlikely with that distance and both isolates rather share a common ancestor/linked case that was not sampled in the study period or study area. This common ancestor would be then in a four and five SNP distance, respectively that is translated to an infection event eight and ten years ago assuming an evolutionary rate of » 0.5 SNPs/genome/year (Ford et al., 2011;Walker et al., 2013;Roetzer et al., 2013;Walker et al., 2014). The central isolate has five other sampled isolates in proximity, which might indicate a super spreader patient and/or a particularly transmissible strain.

Genotypic drug resistance prediction
Mutations (small deletions and SNPs) in 34 resistance associated target regions (comprising 28 genes) were considered for a molecular resistance prediction to 13 first-and second-line drugs (Supplementary file 1). Mutations in genes coding for the RNA-Polymerase subunits rpoA, rpoB (excluding resistance mediating mutations), and rpoC were reported as putative fitness compensating (e.g. in vitro growth enhancing) variants for rifampicin resistant strains. A detailed overview of all mutations considered as genotypic resistance marker is given as Supplementary file 1. Mutations that were not clearly linked to phenotypic drug resistance were reported as genotypic non wild type and were not considered as genotypic resistance markers. When no mutation (or synonymous, silent mutations) was detected in any of the defined drug relevant target regions the isolate was considered to be phenotypically susceptible.

Phylogenetic inference (maximum likelihood)
We used jModelTest v2.1 and Akaike and Bayesian Information Criterion (AIC and BIC) to find an appropriate substitution model for phylogenetic reconstructions based on the concatenated sequence alignments (Appendix-table 5). Maximum likelihood trees were calculated with FastTree 2.1.9 (double precision for short branch lengths) (Price et al., 2010) using a general time reversible (GTR) nucleotide substitution model (best model according to AIC and second best model according to BIC), 1000 resamplings and Gamma20 likelihood optimization to account for evolutionary rate heterogeneity among sites. The consensus tree was rooted with the 'midpoint root' option in FigTree and nodes were arranged in increasing order. Variants considered as drug resistance marker (see above) and putative compensatory variants were analyzed individually and mapped on the phylogenetic tree to define resistance patterns of identified phylogenetic subgroups.

Molecular clock model
In order to compute a time scaled phylogeny and employ the Bayesian skyline model (see below) for the identified Central Asian outbreak (CAO) clade we sought to define an appropriate molecular clock model (strict versus relaxed clock) and a mutation rate estimate. Due to the restricted sampling timeframe of the Karakalpakstan dataset (2001 -2006) we extended the dataset for the model selection process with CAO strains from Samara (2008 -2010) and 'historical' CAO strains isolated from MDR-TB patients in Germany (1995 -2000) thus allowing for a more confident mutation rate estimate. The strength of the temporal signal in the combined dataset, assessed by the correlation of sampling year and root-to-tip distance, was investigated with TempEst v1.5 (44). Regression analysis was based on residual mean squares, using a rooted ML tree (PhyML, GTR substitution model, 100 bootstraps), R-square and adjusted P-value are reported. For the comparison of different Bayesian phylogenetic models we used path sampling with an alpha of 0.3, 50% burn-in and 15 million iterations (resulting in mean ESS values > 100), marginal likelihood estimates were calculated with BEAST v2.4.2 (45), and D marginal L estimates are reported relative to the best model.
First, we employed a strict molecular clock fixed to 1 Â 10 À7 substitutions per site per year as reported previously (Ford et al., 2011;Walker et al., 2013;Roetzer et al., 2013) without tip dating, a strict molecular clock with tip dating and a relaxed molecular clock with tip dating. BEAST templates were created with BEAUti v2 applying a coalescent constant size demographic model, a GTR nucleotide substitution model, a chain length of 300 million (10% burn-in) and sampling of 5000 traces/trees. Second, we ran different demographic models (i.e. coalescent constant size, exponential, and Bayesian skyline) under a relaxed molecular clock using tip dates and the same parameters for the site model and Markov-Chain-Monte-Carlo (MCMC) as described above. Inspection of BEAST log files with Tracer v1.6 showed an adequate mixing of the Markov chains and all parameters were observed with an effective sample size (ESS) in the hundreds, suggesting an adequate number of effectively independent draws from the posterior sample and thus sufficient statistical support.

Bayesian Skyline Plot
Changes of the effective population size of the CAO clade in Karakalpakstan over the last four decades were calculated with a Bayesian skyline plot using BEAST v2.4.2 (45) using a tip date approach with a strict molecular clock model of 0.94 Â 10 À7 substitutions per site per year (best model according to path sampling results, see above), and a GTR nucleotide substitution model. We further used a random starting tree, a chain length of 300 million (10% burn-in) and collected 5000 traces/trees. Again adequate mixing of the Markov chains and ESS values in the hundreds were observed. A maximum clade credibility genealogy was calculated with TreeAnnotator v2.

Impact of resistance-conferring and compensatory mutations on transmission success
We used multiple linear regression to examine the respective contributions of antimicrobial resistance and putative fitness cost-compensating mutations to the transmission success of tuberculosis. To take transmission duration into account, we computed, for each isolate and each period length T in years (from 1 to 40y before sampling), a transmission success score defined as the number of isolates distant of less than T SNPs, divided by T. This approach relied on the following rationale: based on MTBC evolution rate of 0.5 mutation per genome per year, the relation between evolution time and SNP divergence is such that a cluster with at most N SNPs of difference is expected to have evolved for approximately N years. Thus, transmission success score over T years could be interpreted as the size of the transmission network divided by its evolution time, hence as the average yearly increase of the network size. For each period T, the transmission success score was regressed on the number of resistance mutations and on the presence of putative compensatory mutations. The regression coefficients with 95% confidence intervals were computed and plotted against T to identify maxima, that is, time periods when the transmission success was maximally influenced by either resistance-conferring or -compensating mutations. These analyses were conducted independently on outbreak strains of the Beijing-CAO clade in the Karakalpakstan cohort and of the Beijing-A clade in the Samara cohort.
isolates exhibit more phenotypic drug resistances compared to non-Beijing isolates (p=0.0079) and more genotypic drug resistances compared to other Beijing isolates (p<0.0001), and non-Beijing isolates (p<0.0001). P-values for pairwise comparison with reference group calculated with unpaired t-test (two-tailed, Welch's correction). DOI: https://doi.org/10.7554/eLife.38200.014 Appendix 1-figure 4. Linear regression analysis showing correlation between root-to-tip distance and sampling years of an extended collection of 220 Beijing CAO datasets covering the period 1995 to 2009.