1. Evolutionary Biology
Download icon

Compensatory evolution drives multidrug-resistant tuberculosis in Central Asia

  1. Matthias Merker  Is a corresponding author
  2. Maxime Barbier
  3. Helen Cox
  4. Jean-Philippe Rasigade
  5. Silke Feuerriegel
  6. Thomas Andreas Kohl
  7. Roland Diel
  8. Sonia Borrell
  9. Sebastien Gagneux
  10. Vladyslav Nikolayevskyy
  11. Sönke Andres
  12. Ulrich Nübel
  13. Philip Supply
  14. Thierry Wirth
  15. Stefan Niemann  Is a corresponding author
  1. Research Center Borstel, Germany
  2. German Center for Infection Research, Germany
  3. Evolution Moléculaire, Ecole Pratique des Hautes Etudes, PSL University, France
  4. Institut de Systématique, Evolution, Biodiversité, UMR-CNRS 7205, Muséum National d’Histoire Naturelle, Université Pierre et Marie Curie, Ecole Pratique des Hautes Etudes, Sorbonne Universités, France
  5. University of Cape Town, South Africa
  6. CIRI INSERM U1111, University of Lyon, France
  7. Schleswig-Holstein University Hospital, Germany
  8. Swiss Tropical and Public Health Institute, Switzerland
  9. University of Basel, Switzerland
  10. Imperial College London, United Kingdom
  11. Public Health England, United Kingdom
  12. National Tuberculosis Reference Laboratory, Research Center Borstel, Germany
  13. Leibniz-Institut DSMZ- Deutsche Sammlung von Mikroorganismen und Zellkulturen, Germany
  14. Université de Lille, CNRS UMR 8204, Inserm U1019, CHU de Lille, Institut Pasteur de Lille, Centre d'Infection et d'Immunité de Lille, France
  15. Centre National de la Recherche Scientifique, Unité Mixte de Recherche, Center for Infection and Immunity of Lille, France
  16. Université de Lille Nord de France, France
  17. Institut Pasteur de Lille, France
Research Article
  • Cited 0
  • Views 564
  • Annotations
Cite this article as: eLife 2018;7:e38200 doi: 10.7554/eLife.38200

Abstract

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.

https://doi.org/10.7554/eLife.38200.001

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.

https://doi.org/10.7554/eLife.38200.002

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 (2008–2010) (Casali et al., 2014a) to investigate clonal relatedness, resistance and compensatory evolution in both settings.

Results

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–2006) (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 CAO-isolates was similar between 2001–2002 and 2003–2004 (49% and 52%, respectively), 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).

Drug resistance and transmission success among MDR-MTBC isolates from Karakalpakstan, Uzbekistan.

Maximum likelihood phylogeny (GTR substitution model, 1000 resamples) of 277 MDR-MTBC isolates from Karakalpakstan, Uzbekistan sampled from 2001 to 2006. 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.

https://doi.org/10.7554/eLife.38200.003

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-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 Appendix—figure 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–2006). A linear regression analysis showed correlation between sampling year and root-to-tip distance and even a moderate temporal signal (p=0.00039, R2 = 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 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–1982). The time to the MRCA was confirmed with the exponential growth demographic model (1977, 95% HPD 1977–1982, Appendix—table 2). The MRCA already exhibited a streptomycin resistance mutation (rpsL K43R) (Figure 2), and acquired isoniazid resistance (katG S315T) in 1977 (95% HPD 1973 – 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–1989), and an MDR-genotype (rpoB S450L) around 1986 (95% HPD 1985–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–1996) (Figure 2).

Evolutionary history of MTBC Central Asian outbreak (CAO) strains.

Genealogical tree of CAO strains in Karakalpakstan, Uzbekistan and effective population size over time based on a (piecewise-constant) Bayesian skyline approach using the GTR substitution model and a strict molecular clock prior of 0.94 × 10−7 substitutions per nucleotide per year. Pink shaded area represents changes in the effective population size giving the 95% highest posterior density (HPD) interval with the pink line representing the mean value. Vertical lines indicate time points of the implementation of the first standardized TB treatment program (DOTS) in Karakalpakstan and of the declaration of Uzbekistan as independent republic. Symbols on branches show steps of fixation of resistance conferring mutations.

https://doi.org/10.7554/eLife.38200.004

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 1994 (95% HPD 1992–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).

Impact of compensatory variants on transmission networks

Overall, 62.1% (172/277) of all MDR-MTBC isolates carried putative compensatory mutations (Comas et al., 2012; Casali et al., 2014a) in rpoA (n = 7), rpoC (n = 126) and rpoB (n = 43) (Supplementary file 1). These mutations were almost completely mutually exclusive, as only 4/172 isolates harbored variants in more than one RNA polymerase-encoding gene. While mutations in rpoA and rpoB were equally distributed between Beijing-CAO isolates and other non-outbreak Beijing isolates, CAO-isolates had more rpoC variants (56% vs 28%, p=0.003) (Appendix—table 3). The mutation rpoC N698N accounted for 79/124 (63.7%) of CAO isolates with putative compensatory effects. The mean number of resistance mutations was higher among isolates carrying compensatory mutations (Figure 3A), 4.77 vs 3.35 mutations (two-sample t-test p=1.2×10−10). Notably, isolates with compensatory mutations also showed larger transmission indexes than isolates presenting no compensatory mutation, 37.16 vs 9.22 (Welch two-sample t-test p<2.2×10−16) (Figure 3B). CAO-isolates with compensatory mutations also had more resistance-conferring mutations than CAO-isolates lacking such mutation (ANOVA, Tukey multiple comparisons of means P adj = 0.0000012). There was no difference observed for the means of resistance-conferring mutations among non-CAO isolates; compensatory mutation present vs. absent (P adj = 0.1978623) (Figure 3C).

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.

https://doi.org/10.7554/eLife.38200.005

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.

Contributions of resistance-conferring and compensatory mutations to the transmission success of the MTBC Beijing-CAO clade, Karakalpakstan, Uzbekistan.

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.

https://doi.org/10.7554/eLife.38200.006

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.

MDR-MTBC phylogeny and resistance mutations of isolates from Samara (Russia) and Karakalpakstan (Uzbekistan) Maximum likelihood tree (with 1,000 resamples, GTR nucleotide substitution model) based on 12,567 variable positions (SNPs) among 705 MDR-MTBC isolates from Karakalpakstan and Samara.

Any resistance associated mutations (see methods) for individual antibiotics are depicted with red bars for each isolate. 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.

https://doi.org/10.7554/eLife.38200.007

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).

Similar to Karakalpakstan, MDR-MTBC isolates from Samara with compensatory mutations also accumulated more resistance-associated mutations (4.57 vs 2.30 mutations per genome; two-sample t-test p<2.2×10−16) and had higher transmission indexes (50.32 vs 0.46; Welch two-sample t-test p<2.2×10−16) compared to isolates lacking compensatory mutations (Appendix—figure 6).

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).

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.

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.

Materials and methods

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.

Drug susceptibility testing

Drug susceptibility testing (DST) was performed for five first-line drugs (isoniazid, rifampicin, ethambutol, streptomycin, pyrazinamide), and three second-line drugs (ofloxacin, capreomycin and prothionamide) for cohort 1, and six second-line drugs for cohort 2 (capreomycin, amikacin, ofloxacin, ethionamide, para-aminosalicylic acid and cycloserine) by the reference laboratory in Borstel, Germany as described previously (Kent and Kubica, 1985).

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 (Li and Durbin, 2009). Alignments were refined with GATK (McKenna et al., 2010) and Samtools (Li et al., 2009) 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 FigTree (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–2006), we extended the dataset for the model selection process with CAO isolates from Samara (2008–2010) and ‘historical’ CAO isolates 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 Δ 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.

Appendix 1

Supplementary Methods

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 Δ 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.

Appendix Figures and Tables

Appendix 1—figure 1
Box-Plot showing pairwise SNP distances among identified Beijing clades and other Beijing isolates in comparison to non-Beijing isolates from Karakalpakstan, Uzbekistan.

Box represents inter quartile range, whiskers represent 95% of the data, outliers shown as black dots; solid black line represents the median.

https://doi.org/10.7554/eLife.38200.012
Appendix 1—figure 2
Proportions of different Beijing clades, other Beijing isolates and non-Beijing isolates in Karakalpakstan, Uzbekistan stratified to the years 2001/02, 2003/04, 2005/06.

P-values for pairwise comparisons within groups were calculated with Fisher exact test (two-sided). Beijing CAO 2001/2002 vs 2005/2006 p=0.0018, Beijing CAO 2003/2004 vs 2005/2006 p=0.0002.

https://doi.org/10.7554/eLife.38200.013
Appendix 1—figure 3
Box-Plot showing median number of (A) phenotypic and (B) genotypic drug resistances (in addition to the MDR classification,that is isoniazid and rifampicin resistance) of all isolates from Karakalpakstan.

Box represents inter quartile range, whiskers represent 95% of the data, outliers shown as black dots; solid black line represents the median. Beijing CAO 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).

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.
https://doi.org/10.7554/eLife.38200.015
Appendix 1—figure 5
Number of drug resistance mutations among different MDR-MTBC groups from Samara (n = 428) and Karakapakstan (n = 277).

Box-Plot with mean (diamond) and median (horizontal line) number of genotypic drug resistances (see methods) to additional anti-TB drugs (beyond MDR defining rifampicin and isoniazid resistance). Box represents inter quartile range, whiskers represent 95% of the data, outliers shown as black dots. P-values for three major Beijing outbreak clades (A, B and CAO), and non-Beijing isolates (mainly lineage four isolates) were calculated with unpaired t-tests with Welch correction compared to the group ‘other Beijing’. Color codes according to Figure 5. P-values for pairwise comparison with reference group calculated with unpaired t-test (two-tailed, Welch’s correction). Clade A (p≤0.0001), Clade B (p=0.0143), CAO (p≤0.0001), and non-Beijing (p=0.0009).

https://doi.org/10.7554/eLife.38200.016
Appendix 1—figure 6
Comparisons between isolates carrying compensatory mutations (in orange) and isolates with no-compensatory mutations (in blue), from the Samara 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<2.2×10−16). (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 function of the number of isolates. (C) Density plot of the number of resistance-conferring mutations for isolates carrying compensatory mutations (orange) and isolates that don’t carry compensatory mutation (blue) from Samara dataset. Proportions are adjusted by using Gaussian smoothing kernels.

https://doi.org/10.7554/eLife.38200.017
Appendix 1—figure 7
Contributions of resistance-conferring and compensatory mutations to the transmission success of M.tuberculosis of the Beijing-A clade from Samara, Russia.

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. Compensatory mutations were independently associated with transmission success, with a maximum association strength found for transmission clusters beginning around 1999.

https://doi.org/10.7554/eLife.38200.018
Appendix 1—figure 8
Exemplary minimum spanning tree to visualize the determination of a transmission index for each isolate.

Each node/isolate is labelled with its transmission index, i.e. number of other isolates with a maximum distance of 10 SNPs. Branches are indicated with number of unique SNPs.

https://doi.org/10.7554/eLife.38200.019
Appendix 1—table 1
Main characteristics of patients from cohorts 1 and 2 in Karakalpakstan, Uzbekistan.
https://doi.org/10.7554/eLife.38200.020
Cohort 1Cohort 2P value
Year of isolate collection
(patient diagnosis with MDR-TB)
2001–20022003–2006
No. MDR cases diagnosed within time period57300
No. Included in this analysis
Reasons for non-inclusion:
Multiple strain infection
No DNA available
Patient already in cohort 1
Low DNA quantity
49 (86%)
6
2
NA
0
228 (76%)
1
40
11
20
0.094
Patient residence (within Karakalpakstan)
Nukus
Chimbay
Other
Unknown
34 (69%)
6
9
0
146 (64%)
64
1
17
0.49
Male27 (55%)119 (52%)0.72
Age (median, IQR)
Missing age
32, 27–38
0
31, 24–41
49 (21%)
0.40
Previous TB treatment38 (78%)228 (100%)<0.0001
First-line resistance profile:
HR
HRE
HRS
HRES
HRSZ
HREZ
HRESZ
1
0
12 (24%)
28 (57%)
1 (2%)
1
7 (14%)
2
1
41 (18%)
49 (21%)
27 (12%)
1
107 (47%)
<0.0001
No. of first-line drugs resistant
2
3
4
5
1
12 (24%)
30 (61%)
7 (14%)
2
42 (18%)
77 (34%)
107 (47%)
<0.0001
Availability of second-line drug susceptibility testing (DST)Ofx, Cap, ProthOfx, Cap, Ami, Eth, Cyc, PAS
Ofx resistance5 (10%)6 (3%)0.033
Cm resistance1 (2%)53 (23%)0.0001
  1. Abbreviations: H = isoniziad, R = rifampicin, E = ethambutol, S = streptomycin, Z = pyrazinamide, Ofx = ofloxacin, Cm = capreomycin

Appendix 1—table 2
Path sampling results and model selection based on Δ marginal L estimates (relative to best model (ref)) considering 75 path sampling steps and chain lengths of 15 million.

Comparisons of Beast runs using a combined dataset of Central Asian outbreak (CAO) isolates originated from Germany (1995 – 2000), Karakalpakstan (2001–2006), and Samara (2008 – 2010) as well as comparisons of Beast runs from CAO-Karakalpakstan data employing the best clock model/substitution rate estimate and runs with the 95% HPD intervals for the substitution rate. Mean node age of CAO-MRCA and acquisition of putative compensatory mutation rpoC N698S of CAO-Karakalpakstan clade are given for each model.

https://doi.org/10.7554/eLife.38200.021
Clock modelDemographic modelMarginal L estimateMean
ESS
Δ marginal L estimateSubst rate x 10−7 (95%HPD)MRCA and rpoC N698S mean node age
(95%HPD)
Combined CAO dataset for clock model comparison
Strict
(no tip dating)
Coalescent constant size−10131.67430232.211.0
(fixed)
41.5 (30.6–49.1)
NA
Strict
(tip dating)
Coalescent constant size−10099.464041ref1.0
(fixed)
42.9 (34.3–50.3)
NA
Relaxed, lognormalCoalescent constant size−10117.21130317.750.96
(0.65–1.24)
57.6 (34.4–84.5)
NA
Combined CAO dataset for molecular clock estimate among CAO strains
Relaxed, lognormalCoalescent constant size−10117.21130378.280.96
(0.65–1.24)
57.6 (34.4–84.5)
NA
Relaxed, lognormalExponential−10044.4112665.480.88
(0.58–1.21)
40.5 (26.4–53.2)
NA
Relaxed, lognormalBayesian skyline−10038.93924ref0.94
(0.72–1.15)
37.1 (25.7–44.0)
NA
CAO-Karakalpakstan dataset for demographic model comparison (under best clock estimate)
Strict
(tip dating)
Bayesian skyline−7617.0928743.790.94
(fixed)
32.2 (23.9–37.3)
16.1 (11.6–16.9)
Strict
(tip dating)
Coalescent constant size−7667.92423154.620.94
(fixed)
37.5 (30.2–45.1)
15.8 (12.3–18.6)
Strict
(tip dating)
Exponential−7613.304003ref0.94
(fixed)
29.3 (23.5–33.7)
15.4 (12.4–20.1)
CAO-Karakalpakstan dataset, exponential growth model, upper and lower 95% HPD values (from best clock estimate)
Strict
(tip dating)
Exponential−7610.943926ref0.72
(fixed)
36.4 (30.4–43.3)
18.4 (15.8–21.6)
Strict
(tip dating)
Exponential−7613.3040032.360.94
(fixed)
29.3 (23.5–33.7)
15.4 (12.4–20.1)
Strict
(tip dating)
Exponential−7621.22403110.281.15
(fixed)
24.4 (20.0–25.8)
12.9 (10.8–14.4)
CAO-Karakalpakstan dataset, skyline model, upper and lower 95% HPD values (from best clock estimate)
Strict
(tip dating)
Bayesian skyline−7611.122694ref0.72
(fixed)
39.5 (30.7–47.8)
18.8 (14.6–21.5)
Strict
(tip dating)
Bayesian skyline−7617.0928745.970.94
(fixed)
32.2 (23.9–37.3)
16.1 (11.6–16.9)
Strict
(tip dating)
Bayesian skyline−7619.7127638.591.15
(fixed)
25.2 (20.2–30.7)
11.8 (10.1–13.9)
  1. Abbreviations: HPD = Highest posterior density interval

Appendix 1—table 3
Mutations in rpoB, rpoA and rpoC associated with a putative compensatory effect in rifampicin resistant MTBC strains.

Data from 277 MDR-MTBC isolates from Karakalpakstan, Uzbekistan, stratified to the particularly successful variant termed Central Asian outbreak (CAO) and other Beijing isolates. Pairwise differences between the two groups calculated with Fisher exact test; two-tailed P-values are reported.

https://doi.org/10.7554/eLife.38200.022
Beijing CAO
(n = 173)
Other Beijing
(n = 64)
P-valueAll
(n = 277)
rpoB mutations outside RRDR, excluding codon 170,400,491 variants
wild type
25 (14.5%)
147 (85.0%)
12 (18.8%)
52 (81.3%)
0.4343 (15.5%)
234 (84.5%)
rpoC variants
wild type
95 (54.9%)
78 (45.1%)
18 (28.1%)
46 (71.2%)
0.0002126 (45.5%)
151 (54.5%)
rpoA variants
wild type
5 (2.9%)
168 (97.1%)
2 (3.1%)
62 (96.9%)
1.007 (2.5%)
270 (97.5%)
  1. Abbreviations: CAO = Central Asian outbreak, RRDR = rifampicin resistance determining region

Appendix 1—table 4
Proportions of genotypic drug resistance rates for different anti-TB drugs (beyond isoniazid and rifampicin resistance) and pre-XDR/XDR-TB classification among 705 MDR-MTBC isolates from Samara (n = 428) and Karakalpakstan (n = 277), stratified to three identified major phylogenetic clades within the Beijing genotype/lineage and to other Beijing isolates, and to non-Beijing isolates (mainly lineage 4, Euro-American).
https://doi.org/10.7554/eLife.38200.023
GroupSEZKmAmCmFqThioPASPre-XDR
XDR
Beijing CAO
(n = 201)
201/201
100.0%
195/201
97.0%
152/201
75.6%
97/201
48.3%
37/201
18.4%
37/201
18.4%
6/201
3.0%
121/201
60.2%
99/201
49.3%
100/201
49.8%
Beijing clade B (W148) (n = 103)103/103
100.0%
83/103
80.6%
44/103
42.7%
61/103
59.2%
18/103
17.5%
18/103
17.5%
23/103
22.3%
75/103
72.8%
12/103
11.7%
64/103
62.1%
Beijing clade A (n = 187)184/187
98.4%
183/187
97.9%
163/187
87.2%
177/187
94.7%
0/187
0.0%
0/187
0.0%
33/187
17.6%
180/187
96.3%
7/187
3.7%
179/187
95.7%
Other Beijing (n = 100)91/100
91.0%
73/100
73.0%
52/100
52.0%
39/100
39.0%
20/100
20.0%
23/100
23.0%
14/100
14.0%
32/100
32.0%
15/100
15.0%
45/187
24.1%
Non-Beijing (n = 114)69/114
60.5%
63/114
55.3%
30/114
26.3%
39/114
34.2%
14/114
12.3%
14/114
12.3%
3/114
2.6%
34/114
29.8%
34/114
29.8%
40/114
35.1%
  1. Abbreviations: S = streptomycin, E = ethambutol, Z = pyrazinamide, Km = kanamycin, Am = amikacin, Cm = Capreomycin, Fq = fluoroquinolone, Thio = thioamide, PAS = para aminosalicylic acid

Appendix 1—table 5
Likelihood scores for different substitution models calculated with Jmodeltest 2.1 and statistical model selection based on Akaike and Bayesian Information Criteration (AIC and BIC).

Best model is assumed to have the lowest criteration value. Shown are the top 10 AIC models. Substitution model used for Bayesian inference marked in bold.

https://doi.org/10.7554/eLife.38200.024
Subst. model-lnLAICΔ AICBICΔ BIC
GTR8837.643718567.28750.021041.00257.0748 (2)
GTR + I8837.674718569.34942.0619 (2)21048.610914.6832 (5)
GTR + G8838.984218571.96844.6809 (3)21051.229917.3022 (6)
GTR + I + G8839.007718574.01536.7278 (4)21058.823324.8955 (8)
TPM1uf8845.42618576.8529.5645 (5)21033.92770.0
TPM1uf + I8845.444618578.889111.6016 (6)21041.51137.5836 (3)
TPM1uf + G8846.735418581.470914.1834 (7)21044.09310.1653 (4)
TPM1uf + I + G8846.769718583.539516.252 (8)21051.708117.7804 (7)
SYM8860.647818607.295540.008 (9)21064.371230.4435 (9)
SYM + I8860.682618609.365242.0777 (10)21071.987438.0596 (12)

References

  1. 1
  2. 2
  3. 3
  4. 4
    Infectiousness, reproductive fitness and evolution of drug-resistant Mycobacterium tuberculosis
    1. S Borrell
    2. S Gagneux
    (2009)
    The International Journal of Tuberculosis and Lung Disease : The Official Journal of the International Union Against Tuberculosis and Lung Disease 13:1456–1466.
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
    The New Profile of Drug-Resistant Tuberculosis in Russia: A Global and Local Perspective: Summary of a Joint Workshop
    1. Institute of Medicine Forum on Drug Discovery, Development, and Translation
    2. Russian Academy of Medical Science
    (2011)
    29–36, Drug-Resistant Tuberculosis in the Russian Federation, The New Profile of Drug-Resistant Tuberculosis in Russia: A Global and Local Perspective: Summary of a Joint Workshop, Washington DC, National Academies Press.
  22. 22
    Public Health Mycobacteriology: A Guide for the Level III Laboratory
    1. P Kent
    2. G Kubica
    (1985)
    US Department of Health and Human Services, Centres for Disease Control: Atlanta, GA.
  23. 23
    Treatment outcomes in Multidrug-Resistant TB pateints in Uzbekistan
    1. M Lalor
    2. S Allamuratova
    3. Z Tiegay
    4. A Khamraev
    5. J Greig
    6. K Braker
    7. P De Cros
    8. A Telnov
    (2011)
    42nd Annual International Union Against Tuberculosis and Lung Disease (IUATLD).
  24. 24
  25. 25
  26. 26
  27. 27
    International Activity Report
    1. Medecins Sans Frontiere
    (2013)
    Medecins Sans Frontieres.
  28. 28
  29. 29
  30. 30
    TB chemotherapy review
    1. VY Mishin
    (2008)
    Pulmonologiya pp. 5–14.
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
    Multidrug-resistant tuberculosis in Uzbekistan: results of a nationwide survey, 2010 to 2011
    1. DJ Ulmasova
    2. G Uzakova
    3. MN Tillyashayhov
    4. L Turaev
    5. W van Gemert
    6. H Hoffmann
    7. M Zignol
    8. K Kremer
    9. T Gombogaram
    10. J Gadoev
    11. P du Cros
    12. N Muslimova
    13. A Jalolov
    14. A Dadu
    15. P de Colombani
    16. O Telnov
    17. A Slizkiy
    18. B Kholikulov
    19. M Dara
    20. D Falzon
    (2013)
    Euro Surveillance : Bulletin Europeen Sur Les Maladies Transmissibles = European Communicable Disease Bulletin, 18, 24176581.
  38. 38
    Organization and Methodology of Pulmonary TB Chemotherapy in Outpatients
    1. USSR Ministry of Health
    (1976)
    Moscow, USSR: USSR Ministry of Health.
  39. 39
    Chemotherapy of Pulmonary TB. Methodological Recommendations
    1. USSR Ministry of Health
    (1983)
    Moscow, USSR: USSR Ministry of Health.
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45

Decision letter

  1. Gisela Storz
    Senior and Reviewing Editor; National Institute of Child Health and Human Development, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Compensatory evolution drives multidrug-resistant tuberculosis in Central Asia" for consideration by eLife. Your article has been reviewed by four peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Gisela Storz as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

The manuscript by Merker et al. presents a detailed analysis of circulating MDR/pre-XDR/XDR strains in a region of Uzbekistan. By expanding the analyses to a previously reported dataset from Samara (Russia) the authors generalize the conclusions to "central Asia". The authors found high transmission of MDR/XDR strains and that high transmission is linked to compensatory mutations. The authors also show that population sizes of the main clade changed over time in parallel to important changes on TB control policies or political/historical events. One major conclusion is that the newly endorsed WHO regimen for MDR-TB will have very limited impact on the region given that strains circulating there are already resistant to many of the relevant antibiotics.

Overall, the manuscript is very well written and the phylodynamic approach to addressing these pertinent questions is timely in terms of both its methodology and its conclusions. We have the following suggestions to improve the manuscript.

Essential revisions:

1) Bayesian model selection

The use of Path Sampling to correctly identify the model is commendable. However, while the model used was a strict clock with a Bayesian skyline, this model was never tested. Thus, the substitution rate selected for the analyses was also not from the best model; it was part of the relaxed lognormal clock with the Bayesian skyline demographic. This may lead to mutation rate differences as outlined below. Related to this, the ESS for each model comparison is quite low. Generally, BEAST analyses should aim for an ESS > 200. For such an important foundation of the paper, the authors should test the utilised model and ensure there is sufficient sampling for higher ESS values.

It is also not stated if the authors ran the finally selected model (strict + skyline) under the prior for comparison to their posterior runs. This should be undertaken to ensure the model is not driving the output, especially in the face of a moderate time signal.

2) Mutation rate

Stemming from the above point, the estimates of the mutation rate seem odd compared to previous estimates. While the overall dating analysis is robust and root-to-tip distance is significant, it looks rather modest in terms of R2. This value indicates a weak clock-like structure and should be noted in the manuscript (e.g. see Discussion and Figure 1 in Duchêne et al. 2016). This weak signal is rightly expected for MTB though, and consequently a mutation rate with broad HPDs is inferred in the Bayesian analyses. This uncertainty surrounding the mutation rate is unfortunately ignored for subsequent analyses due to the strict clock selection. Authors should comment on how this uncertainty is accounted for and why their rate is so much faster than previous estimates, which tend to be closer to 10-8.

3) Phylogenetics and SNP alignment

The higher mutation rate above may be a result of the way the alignment was input for the phylogenetics. It is unclear if authors used SNP alignments or reconstituted whole genome alignments. If the former, the SNP alignments should be corrected for invariant site counts in their ML tree (e.g. by the Stamatakis ascertainment bias correction method) and in the BEAST analyses (e.g. using the constant sites parameter). If not, the branch lengths will likely be incorrect, potentially leading to incorrect date estimates as outlined by Leaché et al. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4604835/). This would also lead to much higher mutation rates compared to whole genome estimates. Authors should redo the analyses with these corrections to ensure their time estimates are correct.

In line with this, while unlikely to be significant, the removal of 28 complete genes from the data is a bit odd and may affect the mutation rate. What of non-DR causing mutations in these genes? Certain genes, such as gyrA contain lineage defining SNPs and would skew SNP distances between isolates if removed (affecting transmission clusters as discussed below). Authors should comment on why this approach was chosen.

4) Transmission clusters

The rationale for the transmission success score and transmission index is not clear. With a SNP rate of 0.5/year we would expect that N SNPs would have evolved over N/0.5 years and therefore a 10 SNP threshold indicates a timeline of 20 years. Please clarify the rationale or where necessary adjust the calculations and figures accordingly for transmission. Additionally, in a setting like this where most of the strains are clustered, it would be of benefit to test if the different transmission clusters are monophyletic, so distance and phylogeny converge to the same delineated clusters.

Authors should also outline how a transmission cluster was defined (i.e. how did transmission indices group together to form delineated clusters). This is important as in the subsection “Impact of compensatory variants on transmission networks”, the authors identify an association between CAO + compensatory with higher number of DR mutations. This is not identified in non-CAO strains even with compensatory mechanisms. This may be because the numbers for CAO strains is "inflated" by including strains from the same transmission cluster.

5) Statistics

The use of statistics throughout the manuscript is very appreciated. However, the authors should justify the use of the t-test when it is unlikely that the underlying data is normally distributed. Authors should either test for normality or apply a Mann-Whitney test instead.

Also, while the CAO higher transmission potential and link to higher numbers of resistant mutations/phenotypes is very clear, the way it is calculated may not be correct. Do the authors take into account every strain in the dataset irrespective to whether they belong to the same cluster of transmission? Clustering may inflate the number of strains with resistance and thus the clades with higher transmission will be more likely to have more resistances. It may be better to choose from each transmission cluster one strain representative of each resistant profile and then re-run the analyses with those cluster representatives plus the unique cases.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Compensatory evolution drives multidrug-resistant tuberculosis in Central Asia" for further consideration at eLife. Your revised article has been favorably evaluated by Gisela Storz (Senior Editor), a Reviewing Editor, and three reviewers.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

While the authors response to reviewer point (4) Transmission clusters clarifies the 10-year timeline based on 0.5 SNPs/year this was not incorporated into the manuscript in the same way. It appears that the key assumption is that the study isolates are not the result of direct transmission and the minimum pairwise SNP difference is from a common ancestor of any two isolates. Including this assumption and indicating this distance is from a common ancestor and therefore the calculation is 0.5 SNPs/yr x 2 genomes would make it easier for readers to follow the calculation of the transmission index.

Similarly, the "transmission index" is also a little bit confusing as they do it per strain, not clade, and thus by using pairwise distances this means that the same case can contribute to the transmission index of several strains what seems at least weird. I may have been missing something, also the new sketch is not clear about whether they use the threshold per strain or per cluster although text suggests per strain.

The number and size of CAO clusters is not reported in the manuscript which would be helpful to understand how inflated some statistics may be due to a particular strain. These could be included in the appendix or where key statistics related to CAO strains are reported indicate how much is due to a specific number of clusters. For example, in the subsection “Impact of compensatory variants on transmission networks”, 56% of CAO-isolates had rpoC variants. Is this due to a single large transmission cluster over-representing rpoC mutations? Or a number of different strains?

A citation of Duchêne 2016 to support the proposed moderate temporal signal would be appreciated.

The first sentence of the subsection “MTBC population structure and transmission rates” differentiates between variants and polymorphisms; however, throughout the remainder of the manuscript it appears these terms are sometimes used interchangeably. Please clarify if there is indeed a difference and adjust terminology if necessary for consistency.

https://doi.org/10.7554/eLife.38200.028

Author response

Essential revisions:

1) Bayesian model selection

The use of Path Sampling to correctly identify the model is commendable. However, while the model used was a strict clock with a Bayesian skyline, this model was never tested. Thus, the substitution rate selected for the analyses was also not from the best model; it was part of the relaxed lognormal clock with the Bayesian skyline demographic. This may lead to mutation rate differences as outlined below. Related to this, the ESS for each model comparison is quite low. Generally, BEAST analyses should aim for an ESS > 200. For such an important foundation of the paper, the authors should test the utilised model and ensure there is sufficient sampling for higher ESS values.

It is also not stated if the authors ran the finally selected model (strict + skyline) under the prior for comparison to their posterior runs. This should be undertaken to ensure the model is not driving the output, especially in the face of a moderate time signal.

To better represent the uncertainties with different substitutions rates and demographic models (associated to Figure 2) we further tested the upper and lower 95% HPD interval for the two best models, i.e. exponential growth and skyline with the strict molecular clock prior. The mean node ages of the most recent common ancestor of the CAO-clade from Karakalpakstan and one major split (acquisition of rpoC N698S mutation, described in the main text) are mentioned in an amended Appendix—Table 2. The mean age of the MRCA was in the range of 1966-1982 with 1974 reported in the main text. Also the subclade rpoC N698S mean node age differed marginally between 1987-1994 with 1990 reported in the main text.

The mentioned ESS values were derived from the path sampling approach. We now mention the mean ESS of all parameters from the respective Beast run. All model comparisons were further run with the same priors. To clarify these points we amended Appendix—Table 2 and the Materials and methods section as follows:

“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.”

Further sentences were added to describe the new results in 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.”

“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 1994 (95% HPD 1992-1996), still years before implementation of the systematic DOTS-program in Karakalpakstan in 1998.”

2) Mutation rate

Stemming from the above point, the estimates of the mutation rate seem odd compared to previous estimates. While the overall dating analysis is robust and root-to-tip distance is significant, it looks rather modest in terms of R2. This value indicates a weak clock-like structure and should be noted in the manuscript (e.g. see Discussion and Figure 1 in Duchêne et al., 2016). This weak signal is rightly expected for MTB though, and consequently a mutation rate with broad HPDs is inferred in the Bayesian analyses. This uncertainty surrounding the mutation rate is unfortunately ignored for subsequent analyses due to the strict clock selection. Authors should comment on how this uncertainty is accounted for and why their rate is so much faster than previous estimates, which tend to be closer to 10-8.

We are referring to a short-term mutation rate of 10-7 which has been confirmed several times for M. tuberculosis outbreak scenarios (Eldholm et al., PNAS 2017, Cohen et al., 2015, Walker et al., 2013, Ford et al., Natgen 2013, Merker et al., 2015, Roetzer et al., 2013) and experimental models (Fortune et al., Natgen 2011). In the above mentioned paper by Duchêne et al. the values for lineage 2 and lineage 4 MTB strains also rather span 10-7, with a slightly higher rate for lineage 2. Slower rates are reported for ancient DNA samples or other dating approaches for instance from Comas et al., Natgen 2013 aiming to infer the evolutionary history of MTB strains. Our short-term rate is in line with other reports (see Figure 5 as overview in Eldholm et al., PNAS 2016) and we confirmed this rate with the extended CAO dataset covering a larger sampling interval and partially overcome the mentioned weak molecular clock signal. We state that this is indeed a ‘moderate’ signal, however the mentioned paper by Duchêne et al., 2016 also concludes that: “Importantly, however, even though they [i.e. M. tuberculosis substitution rates] were relatively low, these rates were sufficiently rapid to be accurately estimated using genomic-scale data.”

Furthermore, the main issue here is not to wonder if M. tuberculosis is a measurably evolving population or not. Simply, the rather narrow time-window of the outbreak limits the power of detection of this signal. This can be easily observed as a general process. When one does plot the regressions of the root-to-tip genetic distance against the sampling time intervals for 36 whole genome data sets from Duchêne et al., 2016 (removing the deviant ancient DNA containing data-sets), there is a highly significant positive correlation (R2 = 0.25). This highlights the fact that the chances and power to detect a molecular clock decrease with narrower sampling dates, but does not imply that there is no MEP signature.

With regard to the overall manuscript length we prefer not to elaborate more on this point further. As mentioned in the response to point 1 we included further results to take uncertainties of substitution rates more into account with respect to node ages.

3) Phylogenetics and SNP alignment

The higher mutation rate above may be a result of the way the alignment was input for the phylogenetics. It is unclear if authors used SNP alignments or reconstituted whole genome alignments. If the former, the SNP alignments should be corrected for invariant site counts in their ML tree (e.g. by the Stamatakis ascertainment bias correction method) and in the BEAST analyses (e.g. using the constant sites parameter). If not, the branch lengths will likely be incorrect, potentially leading to incorrect date estimates as outlined by Leaché et al. (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4604835/). This would also lead to much higher mutation rates compared to whole genome estimates. Authors should redo the analyses with these corrections to ensure their time estimates are correct.

In line with this, while unlikely to be significant, the removal of 28 complete genes from the data is a bit odd and may affect the mutation rate. What of non-DR causing mutations in these genes? Certain genes, such as gyrA contain lineage defining SNPs and would skew SNP distances between isolates if removed (affecting transmission clusters as discussed below). Authors should comment on why this approach was chosen.

The correction methods suggested above are indeed important, particularly when a large fraction of the SNPs are not considered, if the sampling is not representative and the taxa rather divergent. Yet this is not the case here. Furthermore, the implementation of such extra-analyses will probably only affect marginally our results and the corrected parameters would stay within our HPD intervals.

We used a SNP alignment and similar exclusion criteria (e.g. PPE/PGRS genes) as previous manuscripts utilizing M. tuberculosis dating approaches, e.g. Eldholm et al., PNAS 2017, Cohen et al., 2015, Walker et al., 2013, or Ford et al., Natgen 2013. The reported mutation rates are all covering 10-7 substitutions per site per year or translated into an outbreak scenario 0.3-0.5 mutations per genome per year. Thus we’re confident that our variant calling and the derived alignment is robust and results in well supported branch lengths and node ages.

The exclusion of entire drug resistance genes removed 85 phylogenetic informative sites, plus 537 positions (including InDels) that were used for resistance prediction purposes and might have further influenced the tree topology. In order to have the most robust phylogeny not affected by unknown resistance marker under positive selection, we therefore chose to remove all positions within these genes. To clarify this part, we modified the section as follows:

“We considered variants that were covered by a minimum of 4 reads in both forward and reverse orientation, 4 reads calling the allele with at least a phred score of 20, and 75% allele frequency. […] The remaining single nucleotide polymorphisms (SNPs) were considered as valid and used for concatenated SNP alignments.

4) Transmission clusters

The rationale for the transmission success score and transmission index is not clear. With a SNP rate of 0.5/year we would expect that N SNPs would have evolved over N/0.5 years and therefore a 10 SNP threshold indicates a timeline of 20 years. Please clarify the rationale or where necessary adjust the calculations and figures accordingly for transmission. Additionally, in a setting like this where most of the strains are clustered, it would be of benefit to test if the different transmission clusters are monophyletic, so distance and phylogeny converge to the same delineated clusters.

Authors should also outline how a transmission cluster was defined (i.e. how did transmission indices group together to form delineated clusters). This is important as in the subsection “Impact of compensatory variants on transmission networks”, the authors identify an association between CAO + compensatory with higher number of DR mutations. This is not identified in non-CAO strains even with compensatory mechanisms. This may be because the numbers for CAO strains is "inflated" by including strains from the same transmission cluster.

We rather consider the “transmission index” as a gradual measurement of the relatedness of patient isolates in a setting as a proxy or surrogate marker for their recent transmission success. Samples within a defined distance, likely are secondary cases which have been infected within the same network some time ago. It’s not strictly a monophyletic cluster/clade, but the higher the mean “transmission index” values are within a defined clade for instance, the lower the genetic diversity, and the more recent transmission events likely have occurred.

We chose 10 SNPs as maximum pairwise distance which would translate in a recent common ancestor node age of 10 years (assuming 5 individual SNPs per isolate differentiating from their progenitor strain). Thus, this would cover the extent of transmission events that occurred within the last 10 years.

To clarify this term further, we added a sketch to the Appendix 1 (extended Materials and methods section) and one sentence to the main Materials and methods section as follows:

“The transmission index was implemented as a proxy for the recent transmission success of defined clades, i.e. the higher the mean transmission index of a clade, the more transmission events have occurred in the past.”

By its nature the measurement is indeed linked with the number of strains. As mentioned by the reviewer, the resolution of phylogenetic trees, respectively the number of well supported branches declines towards the tips of such outbreaks. To set cut-offs here with regard to genetic distances would also lead to an arbitrary cluster definition. Thus we, decided in favour for this gradual measurement to count the number of phylogenetically related strains without considering the topology at the tips but highlighting hotspots of very low diversity in the outbreak itself. In addition this allowed us to implement a parameter linked to cluster size/low genetic diversity/transmission success for further calculations.

5) Statistics

The use of statistics throughout the manuscript is very appreciated. However, the authors should justify the use of the t-test when it is unlikely that the underlying data is normally distributed. Authors should either test for normality or apply a Mann-Whitney test instead.

Also, while the CAO higher transmission potential and link to higher numbers of resistant mutations/phenotypes is very clear, the way it is calculated may not be correct. Do the authors take into account every strain in the dataset irrespective to whether they belong to the same cluster of transmission? Clustering may inflate the number of strains with resistance and thus the clades with higher transmission will be more likely to have more resistances. It may be better to choose from each transmission cluster one strain representative of each resistant profile and then re-run the analyses with those cluster representatives plus the unique cases.

For the first point, in all comparisons where statistics were used, we compared larger sample sizes (n>40) so the violation of a normal distribution should not cause major problems. In Figure 3 the data almost follow normal distributions (perfect normal distributions do not exist).

For the second point, these are real data, and in an epidemic settings, indeed MDR strains spread faster and therefore generate higher transmission indices. We understand the kind of tautological concern raised by the reviewer. However, applying the strategy he mentioned would probably be misleading. To provide a proxy it is like evaluating the demographic change of a suite of let’s say 20 successful (expanding) populations by only considering one strain per population. This new Bayesian tree analysis with 20 strains would unlikely detect an expansion.

We would also like to enhance that the statistics were done by pooling multiple clusters, an approach that should minimize local deviations.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

While the authors response to reviewer point (4) Transmission clusters clarifies the 10-year timeline based on 0.5 SNPs/year this was not incorporated into the manuscript in the same way. It appears that the key assumption is that the study isolates are not the result of direct transmission and the minimum pairwise SNP difference is from a common ancestor of any two isolates. Including this assumption and indicating this distance is from a common ancestor and therefore the calculation is 0.5 SNPs/yr x 2 genomes would make it easier for readers to follow the calculation of the transmission index.

Similarly, the "transmission index" is also a little bit confusing as they do it per strain, not clade, and thus by using pairwise distances this means that the same case can contribute to the transmission index of several strains what seems at least weird. I may have been missing something, also the new sketch is not clear about whether they use the threshold per strain or per cluster although text suggests per strain.

To clarify, the chosen threshold would cover direct transmission events among the study population but also the connection via hypothetical cases or contacts that could have occurred 10 years ago and which are not sampled. To improve the understanding for this parameter with regard to the timing and pairwise distances we added the following lines in the Materials and methods section:

“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).”

In addition we extended the Appendix section with further explanations on the rationale to use the strain specific “transmission index” instead of a clade dependant cluster definition:

“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 central isolate has 5 other sampled isolates in proximity, which might indicate a super spreader patient and/or a particularly transmissible strain.”

The number and size of CAO clusters is not reported in the manuscript which would be helpful to understand how inflated some statistics may be due to a particular strain. These could be included in the appendix or where key statistics related to CAO strains are reported indicate how much is due to a specific number of clusters. For example, in the subsection “Impact of compensatory variants on transmission networks”, 56% of CAO-isolates had rpoC variants. Is this due to a single large transmission cluster over-representing rpoC mutations? Or a number of different strains?

As mentioned above in the new text passages, we did not define strict clusters within the outbreak but a gradual measurement (i.e. transmission index) per isolate that indicates in which parts of the outbreak topology the majority of transmission events might have occurred over the last 10 years. With regard to putative compensatory mutations, the two variants rpoC N698S (n=79) and rpoB I488V (n=18) deeply rooted in the CAO phylogeny (see Figure 2), accounted for a large proportion of the CAO cases with putative compensation (n=124). A possible causative link between rpoC N698S with transmission success is then discussed in the manuscript. To address the point from comment 2 we added the following sentence:

“The mutation rpoC N698N accounted for 79/124 (63.7%) of CAO isolates with putative compensatory effects.”

A citation of Duchene 2016 to support the proposed moderate temporal signal would be appreciated.

Manuscript cited:

“A linear regression analysis showed correlation between sampling year and root-to-tip distance and even a moderate temporal signal (P=0.00039, R2= 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).”

The first sentence of the subsection “MTBC population structure and transmission rates” differentiates between variants and polymorphisms; however, throughout the remainder of the manuscript it appears these terms are sometimes used interchangeably. Please clarify if there is indeed a difference and adjust terminology if necessary for consistency.

We are now referring to polymorphism only while introducing the abbreviation of SNP (single nucleotide polymorphisms).

https://doi.org/10.7554/eLife.38200.029

Article and author information

Author details

  1. Matthias Merker

    1. Molecular and Experimental Mycobacteriology, Research Center Borstel, Borstel, Germany
    2. German Center for Infection Research, Partner site Hamburg-Lübeck-Borstel-Riems, Germany
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    Contributed equally with
    Maxime Barbier
    For correspondence
    mmerker@fz-borstel.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1386-2331
  2. Maxime Barbier

    1. Laboratoire Biologie Intégrative des Populations, Evolution Moléculaire, Ecole Pratique des Hautes Etudes, PSL University, Paris, France
    2. Institut de Systématique, Evolution, Biodiversité, UMR-CNRS 7205, Muséum National d’Histoire Naturelle, Université Pierre et Marie Curie, Ecole Pratique des Hautes Etudes, Sorbonne Universités, Paris, France
    Contribution
    Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—review and editing
    Contributed equally with
    Matthias Merker
    Competing interests
    No competing interests declared
  3. Helen Cox

    Division of Medical Microbiology, Institute of Infectious Disease and Molecular Medicine, University of Cape Town, Cape Town, South Africa
    Contribution
    Conceptualization, Formal analysis, Investigation, Writing—review and editing
    Competing interests
    No competing interests declared
  4. Jean-Philippe Rasigade

    1. Laboratoire Biologie Intégrative des Populations, Evolution Moléculaire, Ecole Pratique des Hautes Etudes, PSL University, Paris, France
    2. Institut de Systématique, Evolution, Biodiversité, UMR-CNRS 7205, Muséum National d’Histoire Naturelle, Université Pierre et Marie Curie, Ecole Pratique des Hautes Etudes, Sorbonne Universités, Paris, France
    3. CIRI INSERM U1111, University of Lyon, Lyon, France
    Contribution
    Software, Formal analysis, Investigation, Visualization, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8264-0452
  5. Silke Feuerriegel

    1. Molecular and Experimental Mycobacteriology, Research Center Borstel, Borstel, Germany
    2. German Center for Infection Research, Partner site Hamburg-Lübeck-Borstel-Riems, Germany
    Contribution
    Formal analysis, Investigation, Writing—review and editing
    Competing interests
    No competing interests declared
  6. Thomas Andreas Kohl

    1. Molecular and Experimental Mycobacteriology, Research Center Borstel, Borstel, Germany
    2. German Center for Infection Research, Partner site Hamburg-Lübeck-Borstel-Riems, Germany
    Contribution
    Data curation, Software, Validation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  7. Roland Diel

    Institute for Epidemiology, Schleswig-Holstein University Hospital, Kiel, Germany
    Contribution
    Data curation, Formal analysis, Writing—review and editing
    Competing interests
    No competing interests declared
  8. Sonia Borrell

    1. Department of Medical Parasitology and Infection Biology, Swiss Tropical and Public Health Institute, Basel, Switzerland
    2. University of Basel, Basel, Switzerland
    Contribution
    Writing—review and editing
    Competing interests
    No competing interests declared
  9. Sebastien Gagneux

    1. Department of Medical Parasitology and Infection Biology, Swiss Tropical and Public Health Institute, Basel, Switzerland
    2. University of Basel, Basel, Switzerland
    Contribution
    Writing—review and editing
    Competing interests
    No competing interests declared
  10. Vladyslav Nikolayevskyy

    1. Imperial College London, London, United Kingdom
    2. Public Health England, London, United Kingdom
    Contribution
    Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9502-0332
  11. Sönke Andres

    Division of Mycobacteriology, National Tuberculosis Reference Laboratory, Research Center Borstel, Borstel, Germany
    Contribution
    Formal analysis, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  12. Ulrich Nübel

    1. Microbial Genome Research, Leibniz-Institut DSMZ- Deutsche Sammlung von Mikroorganismen und Zellkulturen, Braunschweig, Germany
    2. German Center for Infection Research, Braunschweig, Germany
    Contribution
    Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  13. Philip Supply

    1. Université de Lille, CNRS UMR 8204, Inserm U1019, CHU de Lille, Institut Pasteur de Lille, Centre d'Infection et d'Immunité de Lille, Lille, France
    2. Centre National de la Recherche Scientifique, Unité Mixte de Recherche, Center for Infection and Immunity of Lille, Lille, France
    3. Center for Infection and Immunity of Lille, Université de Lille Nord de France, Lille, France
    4. Center for Infection and Immunity of Lille, Institut Pasteur de Lille, Lille, France
    Contribution
    Conceptualization, Investigation, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
  14. Thierry Wirth

    1. Laboratoire Biologie Intégrative des Populations, Evolution Moléculaire, Ecole Pratique des Hautes Etudes, PSL University, Paris, France
    2. Institut de Systématique, Evolution, Biodiversité, UMR-CNRS 7205, Muséum National d’Histoire Naturelle, Université Pierre et Marie Curie, Ecole Pratique des Hautes Etudes, Sorbonne Universités, Paris, France
    Contribution
    Conceptualization, Formal analysis, Supervision, Validation, Investigation, Methodology, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared
  15. Stefan Niemann

    1. Molecular and Experimental Mycobacteriology, Research Center Borstel, Borstel, Germany
    2. German Center for Infection Research, Partner site Hamburg-Lübeck-Borstel-Riems, Germany
    Contribution
    Conceptualization, Resources, Software, Supervision, Funding acquisition, Investigation, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    stniemann@yahoo.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6604-0684

Funding

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.

Acknowledgements

We thank: I Razio, P Vock, T Ubben and J Zallet from Borstel, Germany for technical assistance; the national and expatriate staff of Médecins Sans Frontières, Karakalpakstan; Dr. Atadjan and Dr. K Khamraev from the Ministry of Health (Karakalpakstan) for their support. Parts of this work have been supported by the European Union TB-PAN-NET (FP7-223681) project and the German Center for Infection Research and by the Leibniz Science Campus Evolutionary Medicine of the Lung (EvoLung). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Raw sequence data (fastq files) have been deposited at the European Nucleotide Archive (ENA) under the project number ERP000192.

Senior and Reviewing Editor

  1. Gisela Storz, National Institute of Child Health and Human Development, United States

Publication history

  1. Received: May 8, 2018
  2. Accepted: October 2, 2018
  3. Version of Record published: October 30, 2018 (version 1)

Copyright

© 2018, Merker et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 564
    Page views
  • 112
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Evolutionary Biology
    Daria N Shalaeva et al.
    Research Article
    1. Evolutionary Biology
    2. Genetics and Genomics
    John Grey Monroe et al.
    Research Article