The recent rapid expansion of multidrug resistant strains of Mycobacterium tuberculosis Ural lineage 4.2 in the Republic of Moldova =================================================================================================================================== * Melanie H. Chitwood * Caroline Colijn * Chongguang Yang * Valeriu Crudu * Nelly Ciobanu * Alexandru Codreanu * Jaehee Kim * Isabel Rancu * Kyu Rhee * Ted Cohen * Benjamin Sobkowiak ## Abstract The projected trajectory of multidrug resistant tuberculosis (MDR-TB) epidemics depends on the reproductive fitness of circulating strains of MDR *M. tuberculosis (Mtb)*. Previous efforts to characterize the fitness of MDR *Mtb* have found that *Mtb* strains of the Beijing sublineage (Lineage 2.2.1) may be more prone to develop resistance and retain fitness in the presence of resistance-conferring mutations than other lineages. Using *Mtb* genome sequences from all culture-positive cases collected over two years in Moldova, we estimate the fitness of Ural (Lineage 4.2) and Beijing strains, the two lineages in which MDR is concentrated in the country. We estimate that the fitness of MDR Ural strains substantially exceeds that of other susceptible and MDR strains, and we identify several mutations specific to these MDR Ural strains. Our findings suggest that MDR Ural *Mtb* has been transmitting efficiently in Moldova and poses a substantial risk of spreading further in the region. ## Introduction Multidrug resistant tuberculosis (MDR-TB) is an important driver of TB-related morbidity and mortality and a major threat to TB control in a number of settings. In several countries of the former Soviet Union upwards of 20% of new TB cases have MDR-TB1, with the vast majority of new cases of MDR-TB resulting from direct transmission rather than through resistance acquired during treatment2. Determining the epidemiological and genomic factors that influence the spread of these strains is essential to interrupt the transmission of MDR-TB. The ability for *Mycobacterium tuberculosis* (*Mtb*) to survive, reproduce, and transmit can vary greatly among drug-resistant strains3. Drug resistance in *Mtb* is conferred through chromosomal mutation and is often associated with reductions in bacterial fitness in the absence of the selective pressure of antibiotic treatment. However, the degree and the durability of these fitness costs, as well as the rate of acquisition of resistance-conferring mutations4, appear to differ by *Mtb* strain and lineage5–7 and may be mitigated by compensatory mutations8,9. For example, strains belonging to the *Mtb* lineage 2 (L2) Beijing sub-lineage have been associated with a higher rate of drug resistance acquisition without a reduction in fitness10,11, and appear to be a major cause of MDR-TB transmission in many settings. In contrast, the globally dispersed *Mtb* lineage 4 (L4) appears to exhibit more variability in transmissibility between the distinct sub-lineages12. We previously reported extensive local transmission of MDR-TB in the Republic of Moldova driven by three large clades of of highly drug resistant strains13. Two clades of Beijing lineage 2.2.1 strains were found mainly in the semiautonomous region of Transnistria and the third Ural lineage 4.2 clade was present throughout the country. While the Ural lineage has been reported in other former countries of the Soviet Union14,15, it has not previously been associated with significant local transmission10. This apparent widespread transmission of MDR Ural strains in Moldova is a concerning finding. Here we compare the predicted fitness of MDR and pan-susceptible Ural and Beijing sub-lineage strains from the Republic of Moldova. Using phylodynamic approaches, we estimate the fitness of clades of MDR and non-MDR strains from both sub-lineages. We identify key genomic differences between these groups to explore putative mechanisms for transmission success of MDR-TB strains in the region. The genomic characterization of MDR-TB strains is critical for the identification of strains with epidemic potential and can aid in the surveillance of such strains throughout the region. ## Results ### Sample Data We use data from a previously reported prospective countrywide study of *Mtb* in the Republic of Moldova13. These data include demographic data on all culture-positive TB cases among non-incarcerated adults in the country over the period 1 January 2018 to 31 December 2019 and whole genome sequencing of clinical isolates collected at the time of diagnosis. Of 2770 individuals diagnosed over the study period, 2236 consented to participate in the study and had a *Mtb* isolate available for sequencing. Sequence data was obtained for 2220 patients; 386 individuals were excluded due to evidence of polyclonal infections and 1834 patients were included in the final dataset. An *in silico* approach was used to predict lineages and antimicrobial resistance profiles of all strains in the collection16. These predictions aligned with phenotypic testing where available (Supplementary Table 1). There were 804 *Mtb* sequences of the Beijing lineage 2.2.1 (44%), 420 sequences of the Ural lineage 4.2 (23%), and 594 sequences (32%) of other L4 strains. We focused our analyses on Beijing lineage 2.2.1 and Ural lineage 4.2 as these strains were responsible for 97% of the MDR-TB observed in Moldova over the study period. ### Clade identification Timed phylogenetic trees were constructed separately for the Beijing lineage 2.2.1 and Ural lineage 4.2 sequences, calibrated by collection dates at the tips (Figure 1). Next, we looked to identify groups in each tree with distinct demographic or epidemiological histories to conduct phylodynamic analyses and predict clade-level expansions. We used a time threshold to characterize clades that emerged within the last 120 years on the phylogeny, censoring any clades smaller than 20 taxa. The Ural lineage 4.2 taxa were divided into three clades using this time-based approach (Figure 1A), with two of these clades containing mostly non-MDR strains (Ural\_B 1% MDR and Ural_C 0% MDR), and one clade almost completely comprising MDR strains (Ural_A 95.8% MDR). This Ural_A clade also includes almost all MDR Ural strains included in the study population (252/256, 98.4%). With the time-based approach, the Beijing lineage 2.2.1 taxa were divided into multiple clades that are more heterogenous than the Ural clades in terms of the proportion of MDR strains (Figure 1B). We also used an alternative approach to define clades that detects cryptic population structure within phylogenetic trees called *TreeStructure*17. This model compares the observed ordering of ancestors to the expected ordering in a homogenous population. With this method, we identified population structure in the Ural lineage 4.2 tree (Supplementary figure 1A), but not in the Beijing lineage 2.2.1 tree (Supplementary figure 1B). ![Supplementary figure 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/11/10/2023.11.10.23298377/F5.medium.gif) [Supplementary figure 1.](http://medrxiv.org/content/early/2023/11/10/2023.11.10.23298377/F5) Supplementary figure 1. Timed phylogenetic trees of Ural (A) and Beijing (B) strains used in the study. Branches are coloured by structured population clade designations estimated with *TreeStructure*, and the position of MDR (dark blue) and non-MDR (light blue) strains shown in the coloured bar. ![Figure 1.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/11/10/2023.11.10.23298377/F1.medium.gif) [Figure 1.](http://medrxiv.org/content/early/2023/11/10/2023.11.10.23298377/F1) Figure 1. Timed phylogenetic trees of Ural (A) and Beijing (B) strains used in the study. Branches are coloured by the time-based clade designations with corresponding annotated clades names (right). MDR (dark blue) and non-MDR (light blue) strains are indicated with color. ### Markers of clade expansion To investigate comparative rates of expansion in the Ural and Beijing time-based clades, we estimated the Local Branching Index (LBI)18 from a maximum likelihood phylogeny produced using the full *M. tuberculosis* dataset, including taxa from both lineages. LBI uses the length of the branches around each internal and terminal node in a phylogenetic tree to estimate local clade expansions, with large LBI values consistent with more rapid branching and clade growth. Comparing LBIs of terminal nodes, we found that the LBIs of Ural lineage 4.2 MDR taxa were higher on average than Ural lineage 4.2 non-MDR taxa and Beijing lineage 2.2.1 taxa (ANOVA p < 0.001) (Figure 2A). LBI values for Ural lineage 4.2 MDR taxa were also higher on average than LBIs for taxa belonging to other L4 strains in the collection (Supplementary figure 2). Concordantly, the mean LBI in the majority MDR clade Ural_A was 0.026, higher than clades Ural_B (mean LBI = 0.011, p < 0.001) and Ural_C (mean LBI = 0.008, p < 0.001) (Figure 2B). In contrast, the mean LBI of the MDR clade Beijing_A was 0.022, lower than the mean LBI estimate across the remaining clades (0.025, t-test p < 0.001). ![Supplementary figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/11/10/2023.11.10.23298377/F6.medium.gif) [Supplementary figure 2.](http://medrxiv.org/content/early/2023/11/10/2023.11.10.23298377/F6) Supplementary figure 2. (A) Local branching index (LBI) estimates for all lineage 2 and lineage 4 sub-lineages present in Moldova during the study period, separated by MDR and non-MDR strain status. (B) Case counts of sequences collected during the study period of each lineage2 and lineage 4 sub-lineage, separated by MDR and non-MDR strain status. ![Figure 2.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/11/10/2023.11.10.23298377/F2.medium.gif) [Figure 2.](http://medrxiv.org/content/early/2023/11/10/2023.11.10.23298377/F2) Figure 2. (A) Local branching index (LBI) values for all MDR (orange) and non-MDR (blue) strains used in the study for the Beijing and Ural lineages. (B) LBI values for strains in each time-based clade, colored by whether the clade is majority MDR (≥ 90% - orange), non-MDR (<10% - blue), or mixed MDR and non-MDR (yellow). Next, we used SkyGrowth19 to estimate effective population sizes (Ne) of each time-based clade separately in the timed phylogenies. This nonparametric method models the growth rate of the effective population size over discrete time intervals, which has a direct relation to the effective reproduction number (Re). We predicted the time dependent Ne between 1990 and the beginning of sampling in 2018 for each clade (Figure 3). We found that the Ne of the majority MDR clade, Ural\_A, increased 3.0 fold (95% CrI: 1.30, 5.76) over the period 2010-2018, indicating the recent rapid expansion of this clade. The Ural_A increase in Ne was markedly higher than that of the majority non-MDR Ural clades (Ural_B 0.26 [-0.45, 2.1] and Ural_C -0.29 [-0.84, 2.09]), as well as the majority MDR clade, Beijing_A, (1.18 [0.16, 3.12]), over the same approximate time period. ![Figure 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/11/10/2023.11.10.23298377/F3.medium.gif) [Figure 3.](http://medrxiv.org/content/early/2023/11/10/2023.11.10.23298377/F3) Figure 3. Effective population size estimates for each time-based clade, inferred using *SkyGrowth*. Ural lineage clades are in red and Beijing lineage clades in green, with the percentage of MDR strains in each clade shown in the panel header. ### The effective reproduction number of Ural strains We estimated the effective reproduction number (Re) of MDR and Non-MDR Ural strains using a multi-type birth death (MTBD) model implemented in the Bayesian phylodynamic inference software BEAST220. Re is the expected number of secondary infections caused by a single infectious individual; a higher Re suggests a pathogen is spreading more rapidly. Under a birth-death model, new infections are ‘born’ when a new host is infected, ‘die’ when the host recovers, and are sampled with some probability (zero before the beginning of the study period). Infections are stratified into several ‘types’, and the effective reproduction number of each type is the birth rate divided by the death rate. We found that the predicted Re of the MDR Ural lineage 4.2 strains was Re = 2.5 (95% HPD: 1.59, 3.80). This was substantially higher than the non-MDR Ural 4.2 strains, with an estimated Re = 1.06 (1.02, 1.16); across model iterations, the Re of MDR Ural strains was 2.34 (1.55, 3.40) times higher than non-MDR Ural strains. ![Fig 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/11/10/2023.11.10.23298377/F4.medium.gif) [Fig 4:](http://medrxiv.org/content/early/2023/11/10/2023.11.10.23298377/F4) Fig 4: The effective reproduction number (Re) of MDR (light red) and non-MDR (dark red) Ural strains inferred using a multi-type birth death model in *BEAST2*. ### Genomic associations with Ural MDR-TB in Moldova We characterized mutational differences between the MDR Ural isolates and the non-MDR Ural lineage 4.2 and MDR Beijing strains. We identified 70 SNPs and 5 short insertions or deletions (indels) with an allele frequency of ≥ 90% in the Ural MDR strains that occurred at a frequency of < 10% in the non-MDR Ural or MDR Beijing strains (Supplementary table 2). Of the 70 SNPs, 37 were non-synonymous mutations in coding regions, 29 were synonymous coding SNPs and four were in intergenic regions. There was evidence of multiple non-synonymous SNPs in genes that have been previously linked to antimicrobial resistance. The *rpoB* S450L mutation, which was found in high proportions in both MDR Ural and Beijing isolates and in low frequencies in non-MDR strains, has been previously well-characterized as conferring rifampin resistance21. We also found a high proportion of MDR Ural strains harbored a mutation in *rpsL*, K88R, that has been implicated in streptomycin resistance22, though with an associated fitness cost23. Additionally, a mutation in *ethA*, H281P, was identified in a high proportion of the MDR Ural strains and is almost absent in the rest of the population. This SNP has been observed previously in MDR strains from Moldova24, and mutations in this gene have been linked to resistance to the second-line TB antimicrobial, ethionamide25. Interestingly, we indentified a non-synonymous SNP that was present in 98% of our Ural MDR strains in *murA*, a key gene in the synthesis of an essential componenet of the *Mtb* cell wall, peptidoglycan (PG)26. A recent study described a previously unknown mechanistic link between PG synthesis and mutations in *rpoB* 27. The location of the *murA* R110T mutation, near a known active site residuce, raises the potential that this may be a compensatory mutation associated with rifampin resistance. A synonymous SNP in *esxO* at position 2625924 (codon 83), which was fixed in the Beijing lineage, was found to be present in 94% of the MDR Ural strains and only 4% of the Ural susceptible strains. We also identified a synonymous SNP in *esxP* at codon 10 almost exclusively in the MDR Ural strains. The ESX gene family have been shown to play a role in host-pathogen interaction28 and we found previous evidence of a mutation in *esxW* in the Beijing strains circulating in Moldova13. Synonymous SNPs in the *esxO* and *esxP* have been previously identified as potentially under selection in MDR-TB strains elsewhere29. Additional non-synonymous SNPs were identified in genes that may confer a growth or survival advantage, including Rv0355c (*PPE8*) that has been associated with adaptation to host defense mechanisms30, *mmaA4* previously implicated in drug resistance through inhibition of the mycolic acid biosynthetic pathway31, and Rv1835c that appears to confer a growth advantage in *Mtb* when disrupted based on transposon mutagenesis32. Furthermore, we found five indels at a high frequency in MDR Ural strains that were not present in other isolates, including 2bp frameshift insertion in the *ndhA* gene, though this has been previously reported to be a non-essential gene for *Mtb* survival33. Finally, we performed a genome-wide association study (GWAS) to identify SNPs associated with multiple drug resistance in both the Ural and Beijing sublineages (Supplementary table 3). In the Ural strains, we identified only the *rpoB* S450L mutation as significantly associated with MDR strains. We found the same mutation in the Beijing strains, along with the *katG* S315R mutation, which confers isoniazid resistance34, and the *rpsL* K43R mutation, which confers streptomycin resistance22. ### Evidence of highly transmissible Ural MDR-TB strains outside of Moldova The Ural lineage 4.2 has been identified in other countries of the former Soviet Union10,35 and in light of this, we compared our strains to publicly available *Mtb* sequences from Georgia10. Although the fraction of MDR strains among Ural strains isolated in Georgia was far lower than in Moldova (∼6% MDR in Georgia compared to ∼61% MDR in Moldova), there was evidence that nine of the 27 MDR Ural isolates from Georgia had similar mutational profiles to our MDR Ural strains from Moldova (Supplementary table 2). Of the 70 MDR Ural strain-specific SNPs identified above in Moldovan isolates, 59 (84%) of these SNPs were found in approximately one third of MDR Ural strains from Georgia (9 strains), and four of these strains also harbored a further 7 SNPs (94% of the total SNPs). Two of the known AMR associated SNPs in *rpoB* and *rpsL* were found in higher frequencies in the Georgian MDR Ural strains (78% and 52% respectively) and only two SNPs were not present at all in the MDR Ural isolates from Georgia. All the 70 MDR Ural strain-specific SNPs identified in Moldovan samples were found at a low frequency in pan-susceptible Ural strains from Georgia (< 10%). ### Sensitivity Analyses We conducted several additional analyses to determine whether these results were sensitive to modeling assumptions. To account for the evidence that the average mutation rate in *Mtb* L4 is believed to be slower than in Beijing strains4, we reconstructed an additional timed phylogenetic tree for the Ural sub-lineage 4.2 with a fixed mutation rate of 0.3 SNPs/genome/year and recalculated effective population sizes. We found that tree could be subdivided into time-based clades that were similar in size and composition to the those identified in the main analysis, and that the majority MDR Ural clade was still predicted to have a faster growth rate than the non-MDR Ural clades and the Beijing clades in the main analysis (Supplementary figure 3). We also randomly resampled the Ural sub-lineage 4.2 sequences and re-ran the MTBD model. We found that Re estimates were consistent across the main analysis and the re-sampled run (Supplementary figure 4). Finally, we looked at the demographic composition of the timed clades used in this study to determine whether these factors may contribute to the observed differences in transmission. We found that clades had similar distributions of age, sex, and homelessness (Supplementary table 4). The Bejing_A clade had a higher proportion of formerly incarcerated individuals (24.5%) than the overall average (11.4%), though all other clades did not vary greatly in their fraction of formerly incarcerated individuals (Supplementary table 4). ![Supplementary figure 3.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/11/10/2023.11.10.23298377/F7.medium.gif) [Supplementary figure 3.](http://medrxiv.org/content/early/2023/11/10/2023.11.10.23298377/F7) Supplementary figure 3. (A) Timed phylogenetic tree of Ural strains reconstructed using a lower mutation rate of 0.3 SNPs/genome/year, compared to 0.5 SNPs/genome/year in the main analysis, with MDR status and new time-based clade designation shown. (B) Effective population size estimates for each new time-based clade, inferred using *SkyGrowth*, showing the rapid growth of a majority MDR Ural clade (Ural_26), with corresponds to the Ural\_A clade in the main analysis. ![Supplementary figure 4.](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/11/10/2023.11.10.23298377/F8.medium.gif) [Supplementary figure 4.](http://medrxiv.org/content/early/2023/11/10/2023.11.10.23298377/F8) Supplementary figure 4. The effective reproduction number (Re) of MDR (light red) and non-MDR (dark red) randomly down-sampled Ural strains inferred using a multi-type birth death model in *BEAST2*. ## Discussion We investigated the reproductive fitness of the *Mtb* Ural 4.2 and Beijing 2.2.1 sub-lineage strains responsible for most MDR-TB circulating in the Republic of Moldova. We found that the fitness of MDR Ural strains is high relative to pan-susceptible strains of the same lineage, and even higher than strains of the Beijing sub-lineage commonly associated with successful transmission. We also found strong evidence that there has been a recent, rapid expansion of a large MDR clade of Ural lineage 4.2 strains in Moldova in the past 10-15 years. Given the high burden of MDR-TB in Eastern Europe, these results suggest that MDR-TB of the Ural 4.2 sub-lineage should be closely monitored in Moldova and the surrounding region to enhance TB control efforts and prevent widespread transmission. Our findings are consistent with our earlier work in the Republic of Moldova, which highlighted the presence of highly drug resistant Ural lineage 4.2 strains that appeared to be readily transmitting within the population13,36–38. However, the Re estimates for Ural 4.2 strains reported here contrast with a recently published study that found that there was a reduced transmission fitness associated with multi-drug resistance in L4 strains in another former Soviet Union country10. This Georgian study, conducted on strains collected two years prior to our work in Moldova, identified only 27 MDR Ural 4.2 strains, and only 14% of the MDR-TB in the country was in L4 strains. Notably, we found that a small number of the Georgian MDR Ural strains had similar mutational profiles to the Moldovan MDR Ural strains identified in our study. Given the evidence of the recent expansion of MDR Ural strains in Moldova in the past decade, it is possible that the MDR Ural sub-lineage may now be circulating more widely in Georgia and other countries of the former Soviet Union. We identified several genes associated with anti-tuberculosis drug resistance in the MDR Ural strains that were not present in high proportions in MDR Beijing strains, including genes conferring resistance to second-line drugs streptomycin and ethionamide. In addition, we identified several mutations that are associated with improved bacterial survival, including in ESX family of genes which are involved in host-pathogen interactions39. While these findings provide potential mechanisms by which the relative fitness of the MDR Ural strains is maintained or increased, we only looked for the presence of single variants that were at a high frequency in our MDR strains. A more comprehensive investigation could provide further insights, including the specific fitness effects conferred by individual drug-resistant mutations and putative compensatory mutations, along with an analysis of epistatic interactions between multiple variants present concurrently in MDR strains and the functional modelling of protein changes. A limitation of this study was that our dataset included *Mtb* isolates collected over a two-year period, making it difficult to observe longer term trends in strain evolution. For this reason, our data are not sufficient to estimate a mutation rate; we used a fixed mutation rate to estimate time-resolved phylogenies. This choice may impact estimates of the effective population size. We reconstructed the timed phylogeny with the Ural strains using a lower fixed mutation rate and found the same pattern of rapid expansion in the majority MDR Ural clade. Furthermore, the reported clade expansion estimates were inferred using clades that were defined by applying a cut-off to the timed phylogeny of 120 years before the study period, removing any tips with a coalescent time to the clade beyond this threshold. As such, these results may be biased by excluding MDR strains with low transmission fitness and by the specific choice of threshold. Even so, the final clades include the majority of Ural and Beijing strains from our dataset and the overall finding, that MDR Ural strains appear to have a high transmission fitness than susceptible Ural strains, is supported by the LBI estimates that included all strains in the collection. In conclusion, we found that MDR-TB Ural lineage 4.2 strains have significantly higher effective reproductive fitness than non-MDR Ural strains circulating in the Republic of Moldova. There is also evidence of a recent, rapid expansion of a large clade of MDR Ural strains, in contrast to non-MDR Ural and both MDR and non-MDR Beijing clades that do not show the same concerning trajectory. Given that we find that approximately 18% of Ural strains in our study were resistant to fluroquinolones, there is also the risk that the MDR-TB Ural lineage 4.2 strain develops further drug resistance and could evolve to become extensively drug resistant40,41. Tracking the evolution and transmission of MDR-TB Ural strains in countries of the former Soviet Union is vital to reduce the burden of MDR-TB in these regions and to implement effective TB control strategies. ## Online Methods ### Study Population and Sequence Analysis Between January 2018 and December 2019, TB culture-positive individuals in the Republic of Moldova were asked to consent for the use of routinely collected demographic and epidemiological data, and for whole genome sequencing of bacterial isolates. Ethical approval was obtained for this study from the Ethics Committee of Research of the Phthisiopneumology Institute in Moldova and the Yale University Human Investigation Committee (Number 2000023071). Next generation sequencing was carried out using the Illumina MiSeq platform, with the resulting raw reads mapped to the H37Rv reference strain with *BWA* ‘mem’42 and single nucleotide polymorphisms (SNPs) called using the GATK software package43. We identified polyclonal infections using *MixInfect*44 and excluded these samples from further analysis. Lineage typing and *in silico* drug resistance predictions were carried out with *TB-Profiler* v2.8.1416. Full details on the study enrollment, specimen and individual metadata collection, and whole genome sequencing have been described previously13. ### Phylogenetic tree construction and clade identification A maximum likelihood (ML) phylogeny comprising all isolates was produced from a multiple sequence alignment of concatenated SNPs using *RAxML-NG*45 with the ‘GTRGAMMA’ substitution model. We next constructed time-calibrated phylogenies separately for Beijing lineage 2.2 and Ural lineage 4.2 strains using the R package *BactDating*46 and the well-resolved ML tree as input. Timed phylogenies were built using the ‘strictgamma’ clock model, a root finding algorithm,and a fixed mutation rate 0.5 SNPs/genome/year, which corresponds to previously reported estimates for the *M. tuberculosis* mutation rate4. We also conducted a sensitivity analysis by repeating the analysis on the Ural lineage 4.2 strain using a fixed mutation rate of 0.3 SNPs/genome/year to reflect the potentially slower rate reported for Mtb lineage 4. We ran *BactDating* for 1×107 Markov Chain Monte Carlo (MCMC) iterations, thinning by a factor of 10,000, without a coalescent prior. We assigned taxa to clades with a shared evolutionary history using two different phylogenetic methods. First, we used a time-based method to identify clades with 20 or more tips that emerged within 120 years start of the study period. Second, we assigned taxa to clades based on evidence of an underlying population structure in the phylogeny using the R package *treestructure*17. ### Phylodynamic Analyses We estimated the Local Branching Index (LBI) of each tip in the ML phylogeny using the R package *TreeImbalance*47. LBI is a quantitative method to estimate fitness from the shape of phylogenetic tree. LBI estimates are based on the total branch length surrounding each internal and terminal node in a phylogeny, discounting exponentially with increasing distance from the node. Higher values of LBI are associated with rapid expansion and increased fitness. The scale of the exponential discounting, τ, is a function of the relevant neighborhood size around a node. Neher et al. report that τ should be the average pairwise distance within the tree scaled by a factor of 0.062518. We calculate the LBI for each internal and terminal node based on the ML tree and the scaling factor τ and summarize results for each terminal node by drug resistant genotype and by clade. Next, we estimated effective population sizes (Ne) using the R package *skygrowth*19, a nonparametric autoregressive model to estimate effective population sizes through time. We fitted the model to each time clade separately with a Bayesian MCMC approach and ran the model for 100,000 iterations, allowing estimation of Ne over 50 timepoints prior to the beginning of sampling and a precision factor of 1 standard deviation. We directly inferred the effective reproduction numbers (Re) for MDR and non-MDR Ural lineage 4.2 strains using a multi-type birth death (MTBD) model 48. We sampled 200 taxa from the Ural lineage 4.2, maintaining the same proportion of MDR and non-MDR taxa as in the full dataset. The *bdmm* package49 in BEAST220 version 2.7.4 was used to run the MTBD model using the multisequence alignment of concatenated SNPs and collection dates, accounting for invariant sites. We ran the model for 100 million MCMC iterations or until all parameters had an effective sample size greater than 200. ### Genomic Analysis We identified mutational differences between the MDR and non-MDR strains Ural strains and the MDR Ural and MDR Beijing strains by comparing the allele frequency of all SNPs and short insertions and deletions in these groups. We examined genomic differences between the MDR Ural strains from Moldova and known antimicrobial resistance mutations were annotated using in-house lists based on associated variants described previously16,50 (Supplemental table S5). Finally, we conducted a GWAS to find SNPs associated with MDR in both the Ural and Beijing strains using *treeWAS*51. ## Supporting information Supplementary Tables [[supplements/298377_file02.xlsx]](pending:yes) ## Data Availability All data produced in the present study are available upon reasonable request to the authors ## Footnotes * * These authors are joint senior authors on this work * Received November 10, 2023. * Revision received November 10, 2023. * Accepted November 10, 2023. * © 2023, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution-NonCommercial-NoDerivs 4.0 International), CC BY-NC-ND 4.0, as described at [http://creativecommons.org/licenses/by-nc-nd/4.0/](http://creativecommons.org/licenses/by-nc-nd/4.0/) ## References 1. 1.World Health Organization. Annual Report of Tuberculosis. Annual Global TB Report of WHO 8, (2022). 2. 2.Kendall, E. A., Fofana, M. O. & Dowdy, D. W. Burden of transmitted multidrug resistance in epidemics of tuberculosis: a transmission modelling analysis. Lancet Respir. Med. 3, 963–972 (2015). 3. 3.Cohen, T., Sommers, B. & Murray, M. The effect of drug resistance on the fitness of Mycobacterium tuberculosis. Lancet Infect. Dis. 3, 13–21 (2003). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S1473-3099(03)00483-3&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12505028&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000180201400020&link_type=ISI) 4. 4.Ford, C. B. et al. Mycobacterium tuberculosis mutation rate estimates from different lineages predict substantial differences in the emergence of drug-resistant tuberculosis. Nat. Genet. 45, 784–790 (2013). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.2656&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=23749189&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 5. 5.Grandjean, L. et al. Transmission of Multidrug-Resistant and Drug-Susceptible Tuberculosis within Households: A Prospective Cohort Study. PLoS Med. 12, 1–22 (2015). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pmed.1001878&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25710373&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 6. 6.Fenner, L. et al. Effect of mutation and genetic background on drug resistance in Mycobacterium tuberculosis. Antimicrob. Agents Chemother. 56, 3047–3053 (2012). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiYWFjIjtzOjU6InJlc2lkIjtzOjk6IjU2LzYvMzA0NyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIzLzExLzEwLzIwMjMuMTEuMTAuMjMyOTgzNzcuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 7. 7.Borrell, S. & Gagneux, S. Infectiousness, reproductive fi tness and evolution of drug-resistant. Int J Tuberc Lung Dis 13, 1456–1466 (2009). [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19919762&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 8. 8.Coll, F. et al. Genome-wide analysis of multi- and extensively drug-resistant Mycobacterium tuberculosis. Nat. Genet. 50, 307–316 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41588-017-0029-0&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29358649&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 9. 9.Comas, I. et al. Whole-genome sequencing of rifampicin-resistant Mycobacterium tuberculosis strains identifies compensatory mutations in RNA polymerase genes. Nat. Genet. 44, 106–10 (2012). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.1038&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=22179134&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 10. 10.Loiseau, C. et al. The relative transmission fitness of multidrug-resistant Mycobacterium tuberculosis in a drug resistance hotspot. Nat. Commun. 14, (2023). 11. 11.Merker, M. et al. Transcontinental spread and evolution of Mycobacterium tuberculosis W148 European/Russian clade toward extensively drug resistant tuberculosis. Nat. Commun. 13, (2022). 12. 12.Stucki, D. et al. Mycobacterium tuberculosis lineage 4 comprises globally distributed and geographically restricted sublineages. Nat. Genet. 48, 1535–1543 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/ng.3704&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27798628&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 13. 13.Yang, C. et al. Phylogeography and transmission of M. tuberculosis in Moldova: A prospective genomic analysis. PLOS Med. 19, e1003933 (2022). 14. 14.Mokrousov, I. Mycobacterium tuberculosis phylogeography in the context of human migration and pathogen’s pathobiology: Insights from Beijing and Ural families. Tuberculosis 95, S167–S176 (2015). 15. 15.Sinkov, V. et al. New epidemic cluster of pre-extensively drug resistant isolates of Mycobacterium tuberculosis Ural family emerging in Eastern Europe. BMC Genomics 19, 1–9 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s12864-017-4368-0&link_type=DOI) 16. 16.Phelan, J. E. et al. Integrating informatics tools and portable sequencing technology for rapid detection of resistance to anti-tuberculous drugs. Genome Med. 11, 1–7 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s13073-019-0650-x&link_type=DOI) 17. 17.Volz, E. M. et al. Identification of Hidden Population Structure in Time-Scaled Phylogenies. Syst. Biol. 69, 884–896 (2020). 18. 18.Neher, R. A., Russell, C. A. & Shraiman, B. I. Predicting evolution from the shape of genealogical trees. Elife 3, 1–18 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.7554/eLife.03671&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25049222&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 19. 19.Volz, E. M. & Didelot, X. Modeling the growth and decline of pathogen effective population size provides insight into epidemic dynamics and drivers of antimicrobial resistance. Syst. Biol. 67, 719–728 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/sysbio/syy007&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29432602&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 20. 20.Bouckaert, R. et al. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 15, 1–28 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pcbi.1006858&link_type=DOI) 21. 21.Kumar, S. & Jena, L. Understanding Rifampicin Resistance in Tuberculosis through a Computational Approach. Genomics Inform. 12, 276 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.5808/GI.2014.12.4.276&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25705170&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 22. 22.Jagielski, T. et al. Screening for streptomycin resistance-conferring mutations in Mycobacterium tuberculosis clinical isolates from Poland. PLoS One 9, (2014). 23. 23.Sun, H. et al. Interaction between rpsL and gyrA mutations affects the fitness and dual resistance of mycobacterium tuberculosis clinical isolates against streptomycin and fluoroquinolones. Infect. Drug Resist. 11, 431–440 (2018). 24. 24.Chesov, E. et al. Emergence of bedaquiline-resistance in a high-burden country of tuberculosis. Eur. Respir. J. 59, 2100621 (2021). 25. 25.Ang, M. L. T. et al. EthA/R-independent killing of Mycobacterium tuberculosis by ethionamide. Front. Microbiol. 8, 1–12 (2017). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3389/fmicb.2017.01020&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28197127&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 26. 26.Kieser, K. J. & Rubin, E. J. How sisters grow apart: Mycobacterial growth and division. Nat. Rev. Microbiol. 12, 550–562 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nrmicro3299&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=24998739&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 27. 27.Patel, Y., Soni, V., Rhee, K. Y. & Helmann, J. D. Mutations in rpoB That Confer Rifampicin Resistance Can Alter Levels of Peptidoglycan Precursors and Affect β-Lactam Susceptibility. MBio 14, (2023). 28. 28.Uplekar, S., Heym, B., Friocourt, V., Rougemont, J. & Cole, S. T. Comparative genomics of ESX genes from clinical isolates of Mycobacterium tuberculosis provides evidence for gene conversion and epitope variation. Infect. Immun. 79, 4042–4049 (2011). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiaWFpIjtzOjU6InJlc2lkIjtzOjEwOiI3OS8xMC80MDQyIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjMvMTEvMTAvMjAyMy4xMS4xMC4yMzI5ODM3Ny5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 29. 29.Grandjean, L. et al. Convergent evolution and topologically disruptive polymorphisms among multidrug-resistant tuberculosis in Peru. PLoS One 12, 1–16 (2017). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pone.0179159&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=28957441&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 30. 30.Sassetti, C. M., Boyd, D. H. & Rubin, E. J. Genes required for mycobacterial growth defined by high density mutagenesis. Mol. Microbiol. 48, 77–84 (2003). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1046/j.1365-2958.2003.03425.x&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12657046&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000181802600007&link_type=ISI) 31. 31.Grzegorzewicz, A. E. et al. A common mechanism of inhibition of the Mycobacterium tuberculosis mycolic acid biosynthetic pathway by isoxyl and thiacetazone. J. Biol. Chem. 287, 38434–38441 (2012). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiamJjIjtzOjU6InJlc2lkIjtzOjEyOiIyODcvNDYvMzg0MzQiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMy8xMS8xMC8yMDIzLjExLjEwLjIzMjk4Mzc3LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 32. 32.Dejesus, M. A. et al. Comprehensive essentiality analysis of the Mycobacterium tuberculosis genome via saturating transposon mutagenesis. MBio 8, (2017). 33. 33.Awasthy, D., Ambady, A., Narayana, A., Morayya, S. & Sharma, U. Roles of the two type II NADH dehydrogenases in the survival of Mycobacterium tuberculosis in vitro. Gene 550, 110–116 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.gene.2014.08.024&link_type=DOI) 34. 34.Ando, H. et al. Identification of katG mutations associated with high-level isoniazid resistance in Mycobacterium tuberculosis. Antimicrob. Agents Chemother. 54, 1793–1799 (2010). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiYWFjIjtzOjU6InJlc2lkIjtzOjk6IjU0LzUvMTc5MyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIzLzExLzEwLzIwMjMuMTEuMTAuMjMyOTgzNzcuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 35. 35.Chernyaeva, E. et al. Genomic variations in drug resistant Mycobacterium tuberculosis strains collected from patients with different localization of infection. Antibiotics 10, 1–12 (2021). 36. 36.Brown, T. S. et al. Evolution and emergence of multidrug-resistant mycobacterium tuberculosis in chisinau, moldova. *Microb*. Genomics 7, (2021). 37. 37.Crudu, V. et al. Nosocomial transmission of multidrug-resistant tuberculosis. Int. J. Tuberc. Lung Dis. 19, 1520–1523 (2015). 38. 38.Wollenberg, K. et al. A retrospective genomic analysis of drug-resistant strains of M. tuberculosis in a high-burden setting, with an emphasis on comparative diagnostics and reactivation and reinfection status. BMC Infect. Dis. 20, 1–12 (2020). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s12879-020-05251-9&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 39. 39.Vaziri, F. & Brosch, R. ESX / Type VII Secretion Systems — An Important Way Out for Mycobacterial Proteins. Microbiol. Spectr. 7, (2019). 40. 40.Cegielski, J. P. et al. Extensive drug resistance acquired during treatment of multidrug-resistant tuberculosis. Clin. Infect. Dis. 59, 1049–1063 (2014). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/cid/ciu572&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25057101&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 41. 41.Shin, S. S. et al. Development of extensively drug-resistant tuberculosis during multidrug-resistant tuberculosis treatment. Am. J. Respir. Crit. Care Med. 182, 426–432 (2010). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1164/rccm.200911-1768OC&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20413630&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000280697400019&link_type=ISI) 42. 42.Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–60 (2009). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bioinformatics/btp324&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19451168&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000267665900006&link_type=ISI) 43. 43.McKenna, A. et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 20, 1297–1303 (2010). [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NjoiZ2Vub21lIjtzOjU6InJlc2lkIjtzOjk6IjIwLzkvMTI5NyI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIzLzExLzEwLzIwMjMuMTEuMTAuMjMyOTgzNzcuYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 44. 44.Sobkowiak, B. et al. Identifying mixed Mycobacterium tuberculosis infections from whole genome sequence data. BMC Genomics 19, 613 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/S12864-018-4988-Z&link_type=DOI) 45. 45.Kozlov, A. M., Darriba, D., Flouri, T., Morel, B. & Stamatakis, A. RAxML-NG: A fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference. Bioinformatics 35, 4453–4455 (2019). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/BIOINFORMATICS/BTZ305&link_type=DOI) 46. 46.Didelot, X., Croucher, N. J., Bentley, S. D., Harris, S. R. & Wilson, D. J. Bayesian inference of ancestral dates on bacterial phylogenetic trees. Nucleic Acids Res. 0044, 1–23 (2018). 47. 47.Fischer, M., Herbst, L., Kersting, S., Kühn, L. & Wicke, K. Tree balance indices: a comprehensive survey. (2021). 48. 48.Kühnert, D., Stadler, T., Vaughan, T. G. & Drummond, A. J. Phylodynamics with Migration: A Computational Framework to Quantify Population Structure from Genomic Data. Mol. Biol. Evol. 33, 2102–2116 (2016). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/molbev/msw064&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=27189573&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F11%2F10%2F2023.11.10.23298377.atom) 49. 49.Scire, J., Barido-Sottani, J., Kühnert, D., Vaughan, T. G. & Stadler, T. Improved multi-type birth-death phylodynamic inference in BEAST 2. bioRxiv 2020.01.06.895532 (2020). 50. 50.World Health Organization. Catalogue of mutations in Mycobacterium tuberculosis complex and their association with drug resistance. (2021). 51. 51.Collins, C. & Didelot, X. A phylogenetic method to perform genome-wide association studies in microbes that accounts for population structure and recombination. PLoS Comput. Biol. 14, 1–21 (2018). [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1371/journal.pcbi.1006223&link_type=DOI)