ARTICLE https://doi.org/10.1038/s41467-021-22705-z OPEN Genomic analyses of Mycobacterium tuberculosis from human lung resections reveal a high frequency of polyclonal infections Miguel Moreno-Molina 1, Natalia Shubladze 2, Iza Khurtsilava2, Zaza Avaliani2, Nino Bablishvili2, Manuela Torres-Puente1, Luis Villamayor3, Andrei Gabrielian4, Alex Rosenthal 4, Cristina Vilaplana 5,6,7, Sebastien Gagneux 8,9, Russell R. Kempker10, Sergo Vashakidze2 & Iñaki Comas 1,11 ✉ 1234567890():,; Polyclonal infections occur when at least two unrelated strains of the same pathogen are detected in an individual. This has been linked to worse clinical outcomes in tuberculosis, as undetected strains with different antibiotic resistance profiles can lead to treatment failure. Here, we examine the amount of polyclonal infections in sputum and surgical resections from patients with tuberculosis in the country of Georgia. For this purpose, we sequence and analyse the genomes of Mycobacterium tuberculosis isolated from the samples, acquired through an observational clinical study (NCT02715271). Access to the lung enhanced the detection of multiple strains (40% of surgery cases) as opposed to just using a sputum sample (0–5% in the general population). We show that polyclonal infections often involve genetically distant strains and can be associated with reversion of the patient’s drug sus- ceptibility profile over time. In addition, we find different patterns of genetic diversity within lesions and across patients, including mutational signatures known to be associated with oxidative damage; this suggests that reactive oxygen species may be acting as a selective pressure in the granuloma environment. Our results support the idea that the magnitude of polyclonal infections in high-burden tuberculosis settings is underestimated when only testing sputum samples. 1 Instituto de Biomedicina de Valencia IBV-CSIC, Valencia, Spain. 2 National Center for Tuberculosis and Lung Diseases of Georgia, Tbilisi, Georgia. 3 FISABIO Public Health, Valencia, Spain. 4 National Institute of Allergy and Infectious Diseases, National Institutes of Health, U.S. Department of Health and Human Services, Maryland, USA. 5 Fundació Institut Germans Trias i Pujol (IGTP), Barcelona, Spain. 6 Universitat Autònoma de Barcelona (UAB), Barcelona, Spain. 7 CIBER of Respiratory Diseases, Madrid, Spain. 8 Swiss Tropical and Public Health Institute, Basel, Switzerland. 9 University of Basel, Basel, Switzerland. 10 Department of Medicine, Division of Infectious Diseases, Emory University School of Medicine, Atlanta, USA. 11 CIBER in Epidemiology and Public Health, Madrid, Spain. ✉email: icomas@ibv.csic.es NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z H ow Mycobacterium tuberculosis evolves during the infec- higher than previously recognized5. This is probably also true for tion and treatment of a patient and transmits is key to superinfections with strains resistant to the treatment being used understand phenomena like the spread of drug resistance. against the first infection. So far, the diversity of the pathogen has been mostly studied from The difficulty in identifying polyclonal infections is particularly sputum samples except for one study in the context of HIV/TB pronounced in high-burden MDR-TB countries. Often, similar co-infection1. However, individual sputum samples may under- genotypes are responsible for a large proportion of recent TB estimate the true bacterial diversity within the lung as they are transmission, making it more difficult to distinguish between likely limited to reveal the coexistence of multiple M. tuberculosis these closely related strains when they coexist in the same strains in the same patient. Generally, infections involving two or patient6. In addition, in the context of drug resistance, patterns of more unrelated genotypes can be referred to as polyclonal (Fig. 1). within-host diversity in sputum cultures may be biased towards Polyclonal infections of M. tuberculosis complicate the diagnosis drug-resistant genotypes with high fitness. As a result, sputum- and treatment of tuberculosis (TB), particularly when the based cultures may not reflect the true extent of the pathogen infecting strains differ in their antibiotic susceptibility; which can diversity inside the patient’s lung. For obvious reasons, studying lead to the total replacement of the susceptibility profile during M. tuberculosis directly from the lung is not usually possible. treatment or to heteroresistance2. There are however a few studies, some of them based on post- Polyclonal infections are also relevant to evaluate if and how an mortem biopsies, which suggest that the diversity of M. tuber- initial TB infection protects from a second infection. Recent culosis within the host is higher than what can be detected in experiments in macaques suggest that an initial infection is highly sputa1,7. In addition, bacterial diversity within the lungs of TB protective against reinfection and disease caused by the same patients may affect clinical outcomes as TB drugs are known to M. tuberculosis strain3. However, whether the initial infection differ in their capacity to penetrate into the different types of lung protects against heterologous challenges with a different strain lesions8 or into the variable immune microenvironments of is unknown. Until the advent of molecular epidemiology, the individual granulomas9–11. Therefore, the role of within-host dominant view was that primary TB episodes were due to diversity in general and of polyclonal infections is pivotal to TB endogenous reactivation after years of latent infection. However, control as it has implications at many levels12. First, effective progression to active disease occurs mainly in the first 2 years treatment can be compromised13. Second, public health inter- after infection4, suggesting that many of the episodes of TB fol- ventions may need adjustments in settings with high rates of lowing prolonged exposures are not due to endogenous reacti- exogenous reinfection13. And third, there is a need to understand vation but to reinfection. Molecular epidemiological studies also the role of M. tuberculosis strain variation in the context of new show that in high-burden countries the rate of reinfection is vaccines and the factors behind the limited success of the existing BCG vaccine14,15. The country of Georgia has a yearly incidence of 80 TB cases per 100,000 population. Out of those, 14 are MDR-TB cases (17.5%). Importantly, 12% among all new cases are MDR-TB indicating that transmission plays an important role in the epi- demic of MDR-TB in Georgia16,17. Like in many countries of the former Soviet Union, adjunctive surgical resections are sometimes performed in patients not responding to treatment18,19. Utilizing these resected lung tissue samples, we studied the diversity of M. tuberculosis within different parts of cavitary lesions and com- pared it to the M. tuberculosis diversity seen in sputum samples of the same patients. We achieved this using culture-derived bulk sequencing of the bacteria and detecting minority variants down to 3%. In this work, we show a high frequency of polyclonal infections and important differences between patients in terms of genetic diversity in granulomatous lesions. Some of these differences are driven by an increased mutational supply mediated by host- derived reactive oxygen species (ROS). In most patients infected with multiple strains, these strains differ in their drug resistance profiles. Furthermore, the high genetic distance observed between two strains infecting a given patient suggests an important role of strain genetic diversity in establishing a polyclonal infection. Our results represent a challenge for treatment and control of TB in the setting and highlight a possible limitation of new vaccines against TB. Results Fig. 1 Theoretical scenarios for Mycobacterium tuberculosis polyclonal A high frequency of polyclonal infections in MDR-TB patients infections. Blue represents healthy patients and red represents infected from Georgia. A total of 370 M. tuberculosis cultured isolates from patients. a A transmission event of two strains from an infected individual 275 patients were included in this study (Fig. 2a). Lung surgical to another result in two different genotypes being present in the same samples from 18 patients were also available. For nine patients, two space or sample. b An infected patient on treatment clears infection and samples were analyzed from sputum and caseum, respectively, and gets infected again resulting in two different genotypes present over time. for nine others, we analyzed multiple samples from cavitary c An already infected patient get superinfected with a different genotype. granulomatous lesions, samples from remote tubercular foci, and The second genotype will either coexist with the first one or replace it. visually healthy lung tissue surrounding the cavity. In addition, for 2 NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z ARTICLE Fig. 2 Phylogenetic diversity in Georgia and identification of polyclonal infections. a Phylogeny of all Georgian M. tuberculosis isolates included in this study (L2 in blue, L4 in red). Surgery patients are highlighted in different colors. Drug resistance profiles are illustrated in the inner band (see legend) and transmission clusters represented on the outer band were estimated using a 10 SNP phylogenetic distance threshold. Branch lengths and bootstrap values not represented. b Curves connecting patient samples located in distant places of the phylogeny, suggesting polyclonal infections. The phylogeny and associated data can be browsed at the ITOL website (https://itol.embl.de/tree/161111218247381131580915038). 38 patients, we analyzed two consecutive cultures from sputum isolates in the phylogeny instead of with each other. In total, we samples. In terms of drug resistance profiles, patient samples were found 6 out of 7 polyclonal infection surgical patients involved in 16% pan-susceptible (44), 2% mono-resistant (3 mono-INH, 2 recent transmission clusters (Supplementary Fig. 1). These values mono-SM, 1 mono-EMB), 2% poly-resistant (5), 41% MDR (112), suggest that superinfections are common among the MDR-TB 24% pre-XDR (66) and 15% XDR (42). population of this setting. In regards to the sputum and lung tissue M. tuberculosis clinical isolates from patients undergoing adjunctive surgical TB patients undergoing lung surgery show complex infection resection, a genome-based phylogeny of all bulk-sequenced scenarios. For nine surgery patients, several bacterial cultures cultures revealed that for seven out of eighteen cases, samples from different parts of the cavitary lesion could be analyzed in from the same patient did not cluster together (Fig. 2b). This detail, including caseum (C), inner wall (I), external wall (E), strongly suggested the presence of polyclonal infections. Patients remote nodule (N), and surrounding healthy tissue (H), in harboring different strains in different samples become evident addition to the diagnostic sputum culture (S). Analysis of the M. when considering pairwise genetic distances between these tuberculosis genomic diversity within and around lesions showed strains. It has been suggested that >10 SNPs between isolates of very different patterns across patients (Supplementary Figs. 2 and the same patient are unlikely to emerge by clonal diversification 3). Detailed descriptions of the surgical cohort are provided in from a single M. tuberculosis genotype20. For all seven patients in Supplementary Data 1. which the isolates did not cluster together, the distance between For three surgery patients, the situation was even more the genotypes was much higher than 10 SNPs (range 105 to 1132 complicated. For those patients, we found instances of polyclonal SNPs, Table 1). By contrast, in the eleven patients whose isolates infections of two unrelated genotypes not only between isolates of clustered together, the pairwise SNP distance was between 0 and 2 the same individual but also within the same lesion. For instance, SNPs (Wilcoxon Rank Sum Test; W 0, p-value 0.001). For patient G039 harbored two genotypes in different parts of the example, the caseum isolate from patient G019 clustered 120 same lesion (C, I, E, H), while the nodule (N) and sputum (S) SNPs apart from the sputum isolate of the same patient. Extreme samples harbored a single genotype. Deconvolution of the two examples were patient G036 or G327, in which the sputum and genotypes present in the caseum isolate (see “Methods”) revealed the caseum isolates belonged to different lineages (L2 and L4, that the genotype absent from the nodule and sputum belonged more than 1100 SNPs apart). For three patients, the isolates were to a different sublineage and was 105 SNPs apart. A similar not only located in different parts of the phylogeny but were also phenomenon was observed in patient G240, where two genotypes assigned to different transmission clusters (Supplementary Fig. 1). coexisted in the granuloma center, but only one was present in For patient G324, whose surgical isolates clustered with the the sputum isolate and clustered with another patient’s sputum caseum isolate from patient G327, the sputum isolates clustered isolate. Finally, patient G033 also showed evidence of two with patient G005 and G312 sputum isolates. Patients G019 and genotypes belonging to different L4 sublineages, coexisting in the G036 also had paired isolates clustering independently with other nodule and healthy tissue samples. NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z Table 1 Summary of M. tuberculosis genotypes found in the surgical cohort. Patient Infection G1 Lin. G1 DRP G2 Lin. G2 DRP Dist. C I E H N S G018 Clonal L2.2.10 XDR – – 0 G1 G1 G1 G1 G1 G1 G021 Clonal L2.2.10 Pre-XDR – – 0 G1 – – – – G1 G023 Clonal L4.2.1 Sus. – – 1 G1 – – – – G1 G025 Clonal L2.2.10 Pre-XDR – – 0 G1 – – – – G1 G031 Clonal L2.2.10 Pre-XDR – – 2 G1 G1 G1 – – G1 G032 Clonal L2.2.9 MDR – – 0 G1 – – – – G1 G034 Clonal L2.2.10 Pre-XDR – – 0 G1 G1 G1 G1 G1 G1 G035 Clonal L.2.2.9 Pre-XDR – – 0 G1 G1 G1 G1 G1 G1 G085 Clonal L2.2.10 Pre-XDR – – 2 G1 G1 G1 G1 G1 G1 G330 Clonal L2.2.10 Sus. – – 0 G1 – – – – G1 G335 Clonal L2.2.10 Sus. – – 0 G1 G1 G1 G1 G1 G1 G033 Polyclonal L4.3.3 Pre-XDR L4.2.1 Sus. 611 G1 G1 G1 M M G1 G039 Polyclonal L2.2.9 Pre-XDR L2.2.10 Pre-XDR 105 M M M M G2 G2 G240 Polyclonal L4.8 Sus. L4.8 Sus. 162 M – – – – G1 G019 Polyclonal L2.2.10 MDR L2.2.10 Pre-XDR 120 G1 – – – – G2 G036 Polyclonal L2.2.9 MDR L4.8 Sus. 1132 G1 – – – – G2 G324 Polyclonal L2.2.10 Sus. L2.2.10 Pre-XDR 128 G1 G1 G1 G1 G1 G2 G327 Polyclonal L4.3.3 Sus. L2.2.10 Pre-XDR 1116 G1 – – – – G2 Lineage (Lin.) and drug resistance profile (DRP) of the different genotypes (G1, G2), as well as their genetic distances (Dist.) are shown. Genotype locations are indicated with: G1 genotype 1, G2 genotype 2, M mix, – no sample available. C caseum, I inner wall, E external wall, H healthy tissue, N nodule, S sputum. Taken together, our results show that seven out of the proportion of genotype 1 decreased as we moved out of the eighteen surgery patients (39%) showed evidence of infection by caseum until reaching a ratio of nearly 50:50 in the external two phylogenetically unrelated genotypes, either in the same walls and healthy tissue. The situation was reversed in the sample or in separate samples (Table 1). To compare the high remote nodule where genotype 2 was fixed at 100% frequency. percentage of polyclonal infections in our surgical dataset with The sputum of the patient also contained genotype 2 at a 100% the percentage expected when using only a single sputum frequency. Deconvoluting these two genotypes were challenging diagnostic sample, we scanned for possible polyclonal infections as both belonged to the same lineage 2 sublineages. We could in 218 sputum culture-positive patients from the country of extract the genetic base of genotype 2 from the nodule as it Georgia. As only one sample per case was available we could represented 100% of the culture but genotype 1 always existed as only scan for cases of mixed infection. We used two different a mixture with genotype 2 in the other samples. To deconvolute methods to identify the co-existence of two different genotypes genotype 1, we assigned as a fixed SNP any position that was in in single sputum samples, including a phylogenetic-based the caseum sample at above 75% frequency, thus reflecting the approach specifically designed for this analysis (see “Methods”). expected frequency of SNPs associated with genotype 1. In total, we identified 11/218 (5%) likely cases of co-existence in Comparison with genotype 2 corroborated that both genotypes these patients based on a single sputum sample. Supporting this belonged to the same lineage but different sublineages (2.2.9 and result, when we only took into consideration the sputum isolates 2.2.10), being only 105 SNPs apart. The frequencies of genotype from surgery patients no mixed infections were detected. This is 1 and genotype 2 correlated with resistance mutations to several in contrast with findings in surgical samples where we identify antibiotics at the corresponding frequencies (Fig. 3b, c). Several three cases in which two genotypes co-existed (17% in surgical of the resistance mutations were the same in both genotypes, samples vs. 5% in sputum samples from the general population; including those for isoniazid, rifampicin, ofloxacin, and chi-square 4.39, p-value = 0.036). To increase the power to kanamycin. The combined frequency of different resistance detect additional polyclonal infections in the general population, mutations to ethambutol and streptomycin across the two we analyzed a set of 38 patients with two consecutive sputum genotypes resulted in 100% of the culture is resistant to those samples (Supplementary Table 1) collected from days to months drugs. For other drugs like pyrazinamide, resistance mutations apart (range 1–595 days). In total, we identified seven patients were only found in genotype 1. out of the 38 (18.4%) with evidence of polyclonal infection at some point during their TB disease. This was lower than the Patient follow-up treatment can be compromised by the frequency of polyclonal infections in the surgical samples (39%) polyclonal infection. To analyze the role of polyclonal infections but higher than in the single diagnostic sputum samples (5%). in drug resistance, we used the available DST results for the Despite the study limitations, this suggests that the real extent of surgical cohort’s samples. First, we explored whether there was polyclonal infections cannot be accurately estimated from a variability within a patient in phenotypic DST results. In 13 out of single sputum sample. the 18 surgical patients, the DST results did not change across We chose patient G039 to illustrate how complex infection sites even when the sputum sample was taken into account. In and transmission patterns can be in a high-burden MDR-TB five patients, the DST profiles differed across samples. Sometimes setting. This patient showed a polyclonal infection of two the difference involved only one drug like OFX in G019 but on genotypes in the caseum, the inner and external wall as well as other occasions like patient G327, the DST profile was fully in the healthy tissue, but not in the remote nodule or sputum reversed when comparing the caseum and the sputum sample (Fig. 3a). A detailed analysis of the SNP frequencies suggested (from pre-XDR to pansusceptible in this case). In other instances, that there were two genotypes coexisting in the caseum at a ratio we had polyclonal infections like patient G033, in which the mix of 80:20 (genotype1:genotype2). Tracking of variable positions between the two genotypes in H and N samples (susceptible + across samples of the patient revealed that this ratio varied. The pre-XDR) resulted in an overall profile of pre-XDR masking the 4 NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z ARTICLE Fig. 3 Dynamics of polyclonal infection with two different strains across the granuloma. a Descriptive summary of patient G039 genotypes frequencies and distribution across the lesion. The PCA graph illustrates their separation: while the X-axis explains ~90% of the variance, leaving N and S apart from the rest, the Y-axis only explains ~5%, mainly defined by low-frequency variants not shared by N and S (due to reading depth differences). b Antibiotic resistance-associated mutations across surgical and sputum samples for the two co-existing genotypes. c Deconvolution of the two genotypes by frequency patterns clustering. Resistance-associated variants from (b) are represented in purple. Unassigned subpopulations (in orange) are low-frequency variants that we cannot assign to any of the two genotypes and shared fixed variants (in green) are common to both genotypes. C caseum, I inner wall, E external wall, H healthy tissue, N nodule, S sputum. Source data are provided as a Source Data file. second genotype. Although the number of patients was low, three excluding G021 which showed a high diversity only in sputum). out of five patients with discrepancies in DST vs. genomic pre- Thus, the general assumption that the M. tuberculosis diversity in diction had a polyclonal infection, suggesting that these infections sputum samples is only a subset of the within-patient diversity is can mislead individual treatment of MDR-TB patients. For a not always true. Indeed, we saw patients in which the sputum comparison of phenotypic and genotypic results as well as harbored exclusive SNPs not found in the surgical samples and identification of novel drug resistance markers see Supplementary vice versa. TB patients widely differed in the amount of within- Note 2. host MTB genomic diversity, with approximately half of the Another complication might result from the replacement over patients harboring almost no MTB diversity (e.g., G031, G034, time of one strain by another, which may compromise the G035 with no low-frequency SNPs; Fig. 4a and Supplementary treatment if no additional DST is done during follow-up. In Fig. 3), while others showed a large diversity across sputum and support of this notion, in all four patients harboring different surgical samples (e.g., G023, G025 with >600 low-frequency strains in temporally separated samples (Table 1), the DST profile SNPs; Fig. 4a and Supplementary Fig. 3). However, the com- changed for several drugs. This included cases that were parison between sputum and surgery samples must be taken with susceptible in the first isolate but MDR/XDR in the second caution as the sputum and surgical samples were obtained at isolate or vice versa, like patients G324 or G327. This data shows different time points and may also reflect the sampling of dif- that treatment based on the results from one diagnostic sample ferent lesions. can be misleading and suggests the need for sequential testing of Several selective forces, including the host response and isolates under programmatic conditions, especially in settings antibiotic treatment, may affect the diversity within granulomas with a high MDR-TB burden such as the country of Georgia. from the time of diagnosis to the time of surgery. It has been suggested that the host immune system may exert a mutational Diversity within and around lesions reflects host micro- pressure on the infecting bacteria through the production of ROS. environment pressures. For those patients with both surgical and This has been demonstrated in single colony analyses of cultured sputum samples and with no polyclonal infection detected, we sputum samples showing that an elevated mutational supply can analyzed the M. tuberculosis genetic diversity across their sam- be identified in immunocompetent individuals but not in HIV- ples. The overall diversity of M. tuberculosis within patients was positive individuals21. A hallmark of an elevated mutational rate well represented in sputum samples when compared to the due to ROS is a mutational signature associated with oxidative caseum sample (correlations >90%). We then calculated the damage (increased changes C > T and G > A)22. It could be number of bacterial SNPs exclusive to each sample. These are expected that such a ROS mutational signature was amplified in SNPs that have accumulated since the divergence of the sample surgical samples compared to sputum. As our study was not from the closest isolate in the phylogeny. In general, patients with based on single-colony sequencing but on bulk sequencing, we a high number of exclusive SNPs in the lesions also had a high reasoned that identifying variants present in a few colonies in number of SNPs in the sputum (R2 = 0.9623 between cavity culture should roughly correlate with variants at very low centers and corresponding sputa using exclusive SNPs and frequency in bulk culture sequencing. Thus for this analysis, we NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z Fig. 4 Impact of the granuloma microenvironment on within-host diversity. a Average diversity of each patient’s surgical samples (measured by a number of vSNPs and excluding polyclonal infections, n = number of independent surgical samples available for each patient, range 1–3). Data are presented as mean values ± SD. b Pooled comparisons of ROS signature in the different datasets. Categories include surgery specimens, sputum samples from surgery patients, and serial sputum samples from the same patient. c Individual values of caseum and sputum from surgery patients. Purple lines connect those patients whose differences are statistically significant by the two-tailed X2 test. Source data are provided as a Source Data file. focused on variants with frequency lower than 5% for samples The surgical cohort showed seven out of the eighteen patients with enough depth of coverage and excluding patients showing (39%) with evidence of infection by two phylogenetically unre- polyclonal infection as frequencies in those cases are not lated genotypes, either in the same sample or in separate samples. straightforward to interpret (total eligible patients n = 10). Our This represents a high percentage when compared to the 5% analysis revealed that, as hypothesized, surgical samples showed a found in single sputum patients of our dataset and percentages stronger ROS mutational signature compared to the sputum described by others2. Three surgical patients showed true mixed samples from the same patients (21% vs. 8.6%, pooled analysis, infections (multiple genotypes co-existing) and four are likely the Fig. 4b; chi-square 14.5, p-value < 0.0001). As the sputum and results of superinfections due to multiple transmission events. surgical samples were obtained at different time points, we also The high rate of polyclonal infections in these patients suggests analyzed patients with serial sputum samples to determine the effect deficiencies in infection control in the setting23 and agrees with of time on the diversity observed in sputum. While there was an epidemiological and model data showing that repeated exposure effect of time (S1 3% vs. S2 5.9%; chi-square 22.3 p-value < 0.0001), to infection as seen in high-burden settings increases the risk of the magnitude of the ROS signature in serial samples was similar to reinfections, disease progression, and outbreaks23,24. In addition, that seen in sputum samples from surgical patients (Fig. 4b) but polyclonal infections in the surgical cohort usually involve strains significantly lower than in the surgical samples (chi-square p-value with different drug resistance patterns (5 out of 7 cases), a fact < 0.0001 for all comparisons against surgical samples). Some that can hamper successful treatment if not assessed correctly25. patients had stronger ROS signatures than others (Fig. 4c; 5 out of There are several studies showing that infection with multiple 10 significant at the individual level, G021, G023, G025, G031, strains resulted in a poor treatment outcome since the presence of G032; chi-square p-value < 0.05). This is in accordance with results any undetected resistance during standard treatment can propel by Liu and colleagues21, where 4 out of 18 cases showed a the acquisition of further drug resistance13. Importantly, in our significant increase in ROS-associated mutations. Thus, our results results, polyclonal infections change, and many times fully revert, suggest that the granuloma microenvironment can increase the the drug resistance profile of the patients. Thus, follow-up DST mutational supply for transitions during the within-host evolution samples should be implemented in these settings for better of M. tuberculosis. patient management. Figure 1 shows the theoretical scenarios to explain the natural Discussion history of infection for these patients. The fact that in several Studies on within-host diversity of M. tuberculosis are based on cases the second isolate is in a transmission cluster suggests that cultured samples and very few on non-sputum samples, limiting superinfection is one of the most common mechanisms. How- our capacity to understand diversity patterns at the site of ever, we cannot be sure if two strains were already infecting at infection. In this study, we analyzed an 18-patient cohort, most of baseline while only one is detected in sputum as we do not have the MDR to XDR-TB cases, with available lung surgical samples access to lung samples before treatment. Similarly, we cannot test from Tbilisi (Georgia). We performed bulk sequencing of cul- the hypothesis suggested from mice experiments that a secondary tured bacteria from the surgical resections and sputum isolates of infection drives the progression of an asymptomatic primary every patient, finding an elevated percentage of polyclonal infection via immune response or by expression of resuscitation- infections with an impact on drug resistance in the surgical promoting factors13,26. However, our results contribute to the cohort (39%), which cannot be attained by single sputum sam- evidence from mice data26 (but not from macaques3) that rein- ples. In addition, we demonstrate that patterns of genetic diversity fection with a second strain can be common after a primary differ within lesions and across patients suggesting a role of host infection under certain circumstances and calls to re-evaluate the immune microenvironments. natural history of TB in settings with the rampant transmission. 6 NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z ARTICLE The extent of polyclonal infections also has consequences to can increase the risk of superinfection. In addition, diagnostic understand the protective potential of vaccines. Recent in vivo sputum and surgical samples have several months time difference, reinfection experiments suggested that the first episode of TB was which limited some of our analyses, but we tried to assess the more protective against a second episode than previously effect of sampling time by comparing our surgical dataset with thought6. However, our results suggest that in a high-burden the serial sputa dataset (Supplementary Table 1) in which the setting this might not be the case, and future studies designed to time difference between samples is similar overall. Also, even corroborate this observation will be needed. Differences may be though we had access to surgery samples, analyses were done on related to how well the experimental set-up recapitulates natural cultured samples and not directly on the surgery or sputum infection in such settings. In the experiments on macaques, the sample which can impose different biases34. Finally, a possible same strain was used for re-challenge3. By contrast, in our dataset limitation of this study would be cross-contamination as an we frequently observe different strains, usually from different explanation for the high frequency of polyclonal infections, (sub-)lineages, infecting the same patient. In fact, we did not see although we think it is highly unlikely for a number of reasons. polyclonal infections with strains less than 100 SNP apart despite First, sample collection dates and processing dates from all our phylogenetic approach being designed to identify cases patients are not close in time thus they have not shared the same up to 20 SNP apart. This raises the possibility that previous space. Second, sample homogenization is carried out in a closed infection might indeed protect against superinfection with the special, disposable tubes so that samples cannot be mixed. Third, same or a very similar strain, as observed in the macaque model. genotypes match their DST phenotypes in nearly all cases, The differences with animal experiments may also be related to arguing against a general contamination problem. Lastly, not all the nature of the re-challenge, as in TB transmission hotspots, polyclonal infections are in a transmission cluster and thus they the exposure to second infections is recurrent27. However, co- do not match any other strain processed in the laboratory. morbidities and social determinants may also play an important In summary, we have shown that surgical lung resections part in the susceptibility of the host to second infections. Lastly, from TB patients reveal a more complete picture of the within- many clinical trials and vaccine studies use sputum as their basis patient diversity of M. tuberculosis compared to sputum sam- for distinguishing between relapse and reinfection and our results ples. The high frequency of polyclonal infections found may be suggest that it may be a poor correlate to ascertain lung poly- related to the nature of our patient population, yet this allowed clonal infection. us to study complex infection scenarios with the potential to Beyond multiple genotypes, analysis of the genomic diversity confound diagnosis and/or DST results, as in many cases the within and around lesions showed very different patterns across two genotypes involved in the infection had different drug patients. Some showed almost no diversity while others were resistance profiles. The surgery patients are often in a trans- highly diverse. In some cases (notably G-018) the distribution of mission cluster suggesting that either there are uncontrolled drug resistance variants suggests gradients of antibiotic con- hotspots of transmission shared by MDR-TB patients, there is centration as having been described in recent literature8,28,29. The increased host susceptibility to reinfection by different geno- fact that this pattern is not observed in all patients, may be types, or both. The fact that polyclonal infections usually involve explained by different pharmacodynamics, but also by the time strains of different lineages or sub-lineages suggests that vaccine difference until resection between patients as the granuloma walls preclinical models must take into account the genetic diversity get hardened over time. This would limit the penetration of drugs of the bacteria to assess protection. Finally, as observed in inside, which could affect selective pressure and thus the diversity culture-based sputum studies1,25,35, there are profound differ- we are able to recover from genomic data8. But drug pressure is ences in diversity between patients and here we show that those not the only selective force that the pathogen encounters during differences can also be seen across regions of the cavitary lesions infection. It is known that different immune microenvironments and between patients. In some patients, differences are driven by exist within granulomas30. Chief among macrophage defense a mutational signature associated with ROS and by treatment mechanisms is the production of reactive oxygen and nitrogen and suggest a link between immune and drug resistance selective species31. Looking at sputum cultures, it has been proposed that pressures. Overall, our results exemplify a challenge for tuber- an unintended effect of ROS is to increase the mutation rate of culosis treatment and control, highlighting our knowledge gap the bacteria21. Our data corroborate those results showing that on the natural history of the disease in certain settings and the the mutational signature of ROS, i. e., increased number of need for better patient follow-up in high burden MDR-TB areas transitions, was overall enriched in surgery samples compared to to halt the spread of resistant strains. sputum samples. A recent study suggests that transition bias is linked to the acquisition of drug resistance variants32. It is thus tempting to link ROS-increased mutational supply to an Methods Sampling. We included 275 TB patients in this study (257 with diagnostic sputum accelerated acquisition of drug resistance mutations. Combined and 18 with paired sputum and lung tissue samples, Supplementary Data 1) with a analyses of the bacterial population in the lung during infection in total of 370 M. tuberculosis clinical isolates from between 9 January 2013 and 20 humans and relevant animal models33, drug penetration and March 2018 (see Fig. 5 for details on the surgical cohort). Sputum samples were pharmacodynamics28 and lesion imaging29 will help to better collected according to standard clinical practice. The surgical samples used in this understand the bacterial population diversity driven by host study were acquired within the framework of observational clinical study NCT02715271, which is ongoing as of this publication date but no longer pressures, how it is linked to the emergence and selection of drug- recruiting, and results are reporting on primary outcome #3. Consent to publish resistant subpopulations and, ultimately, to relapses and treat- clinical information potentially identifying individuals was obtained. The main ment failure in clinical settings. indication for lung surgery was the persistence of abnormal lung lesions (pre- Our results are necessarily limited by the characteristics of the dominantly cavities) identified on chest X-ray (13/18), both in DS (3/13) and MDR/XDR-TB (10/13), despite good treatment adherence. Additional indications epidemiological setting and the patients undergoing surgery. for surgery included treatment failure (2.9%), complications of TB related to Most of the patients who underwent surgery had already been pulmonary hemorrhage, spontaneous pneumothorax, or empyema36. Immediately diagnosed at least as MDR-TB and did not respond to treatment. after lung resection, samples were removed from the following zone of obtention: This population is characterized by substantial exposure to cavity center (C), cavity internal wall (I), cavity external wall (E), visually healthy tissue around the cavity (H) and nodule (N), placed in sterile tubes and sent to the multiple antibiotics, a higher frequency of prior TB disease, the National Center for Tuberculosis and Lung Disease microbiological laboratory for frequent presence of cavitary disease (associated with poor anti- processing. In total, we collected 219 single sputa, 38 pairs S1–S2 (serial sputa), 9 biotic penetration), and sometimes prolonged hospitalization that pairs C–S (caseum–sputum), and 9 complete surgical sets (C, I, E, H, N, S). NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z Fig. 5 Details on the surgical cohort. a Sampling sites for the surgical cohort. b Timeline of key events for the surgical cohort. Treatment periods are depicted in red, sputum samples date in blue points, and surgical samples date in green points. Starting case definition and treatment outcome are also represented. Source data are provided as a Source Data file. Culture and drug susceptibility testing. Resected tissue samples were homo- After that, we scanned for optical and PCR duplicates to remove them, as this helps genized using Minilys homogenizer, placed in both LJ and BACTEC 960 MGIT™ to reduce the number of artifactual variants in low-frequency ranges. Variant culture, and processed as per the manufacturer’s instructions. Sputum samples calling for sputum samples was carried out using the software and parameters from were also cultured in both media. All clinical isolates were confirmed for M. the calling module of the pipeline validated for sputum sample cultures at the IBV- tuberculosis complex (MTBC) using the standard microbiological method37,38. All CSIC (available at https://gitlab.com/tbgenomicsunit/ThePipeline). isolates underwent drug susceptibility testing (DST) for first and second-line drugs For a robust variant calling in surgery samples, we used three different variant in both solid (LJ) with increasing drug concentrations and liquid media BBL® callers (VarScan246, GATK’s HaplotypeCaller46,47, LoFreq48) and integrated SNPs MGIT™, as recommended by the manufacturer. reported by at least two of them to get a high-confidence list of low-frequency variants. VarScan2 was run with parameters “pileup2snp sample.pileup—min- coverage 20—min-reads2 4—min-avg-qual 20—min-var-freq 0.01—min-freq-for- DNA extraction and sequencing. For both sputum and resected tissue cultures, hom 0.9—strand-filter 1”, GATK was run with parameters “-T HaplotypeCaller -R after reaching sufficient bacterial growth in the LJ medium, plates were thoroughly ref.fasta -I sample.bam -o sample.vcf—min-base-quality-score 20 -ploidy 1” and scraped to maximize diversity recovery for bulk sequencing. DNA extraction was LoFreq was run with parameters “call-parallel—pp-threads 12 -f ref.fasta -o sample. performed by the standardized protocol in the BSL-2+ laboratory from the AFB vcf sample.bam” and “filter -i sample.vcf -v 20 -A 0.01 -Q 20 -o filtered.vcf”. positive culture. In short, bacteria were heated at 80 °C for 20 min, centrifuged, and From the initial list of variants, we applied a mapping filter that discarded resuspended in TE buffer. Lysozyme, SDS, and proteinase K were added in variants that arose in repetitive genomic regions like the PE/PPE families or sequential incubations. DNA was precipitated using chloroform-isoamyl alcohol phages49. To establish a threshold that discarded additional false variants, we and isopropanol and resuspended in TE buffer39. When there was no sufficient performed a synthetic read simulation. Using the ART software package50, we got growth in LJ, the CTAB/chloroform method was used for DNA extraction from the quality distribution and error profiles of our samples (using 1 mL of MGIT culture. The quality and quantity of DNA were analyzed with the art_profiler_illumina with default parameters), and simulated 100 sequencing runs help of spectrophotometers (Thermo Scientific NanoDrop 2000 and Qubit 3.0). with the data (art_illumina -p −1 profile.txt -2 profile.txt -na -iref.fasta -l 150 -f Extracted DNAs were sent for whole-genome sequencing to Broad Institute, C- 1000 -m 280 -s 137 -o out.fas). By analyzing simulations with the same pipeline, we PATH, TGen North, or Tuberculosis Genomics Unit (IBV-CSIC). The average defined a ~3% minimum frequency threshold to validate a variant in surgery sequencing depth for this study was 141X. The average sequencing depth for samples. In addition, we extended our mapping filter to new regions that showed surgical samples was 267X. high-frequency SNPs in the simulations and were due to systematic mapping errors to the predicted ancestor, especially in Lineage 2 strains (see Supplementary Data 2 Bioinformatic analysis. Read preprocessing was done using fastp40 to scan reads for a list of discarded genomic features). and trim low-quality ends with a mean window quality <20. We then used Kraken41 to taxonomically classify reads by means of a custom database, only keeping MTBC sequences to avoid false variants arising due to contaminant DNA. Phylogenetic and population genetics analyses. A maximum-likelihood phy- Filtered reads were mapped with BWA42 to a predicted MTBC ancestor reference logeny of the 370 samples in the dataset was constructed using IQ-TREE 251 from sequence43 using default parameters, and processed using samtools44 and Picard45. an alignment with gaps and no resistance SNP positions (run parameters: “-m GTR 8 NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z ARTICLE -bb 1000”, meaning the use of a general time-reversible model with unequal rates calculate their genetic distance. For this, we developed an in-house Python script and unequal base frequencies, plus 1000 ultrafast bootstrap replicates). In addition, that uses a phylogenetic marker database from the literature58,59 and takes the SNP an NJ phylogeny was constructed from the same alignment using MEGA52 and file from the sputum sample to perform the detection and separation of both phylogenetic Hamming distances between pairs of samples were extracted by strains, computing their genetic distance solely based on the number of differences. parsing the branch lengths from the tree using Python’s ete3 toolkit53. We parsed SNP files by means of a custom Python script to obtain summary Identification of polyclonal infections in multiple samples tables for each patient and a general one collecting information about the different Phylogeny manual inspection. To identify polyclonal infections happening in two numbers of fixed (fSNPs) and variable SNPs (vSNPs). All figures illustrating different samples from the same patient, we manually analyzed the phylogeny diversity across sites by the patient were produced using R54 and the ggplot2 looking for isolates of the same patient that seat in different parts of the phylogeny. package55. For comparisons within a patient’s samples, we analyzed every SNP in In addition, we recorded if samples were placed close enough to other patients their samples and calculated the differences. This was performed by means of an samples in the phylogeny to suggest a recent reinfection event involving those two in-house Python script that works with the SNP files for every patient, analyzing patients. them sequentially and establishing the population dynamic for every variant based on their presence/absence and frequency, obtaining stable, sporadic, ascending, or descending dynamics. The data were plotted using R and the ggplot2 package as Analysis of frequency spectra differences. To make the detection process systematic well. For multi-sample patients, PCAs were produced with a matrix of SNP for those cases where we had at least two samples (e.g., sputum1-sputum2 or frequencies across all samples and the R’s affycoretools package using the function sputum-caseum), we developed Python and R scripts to perform pairwise com- “plotPCA” for the first two principal components. We also calculated the parisons using SNP frequency differences obtained from the sample’s genomic percentage of explained variation of these two principal components using the data, generating differences profiles that were plotted using ggplot2 package. We “prcomp” function. generated simple XY plots in which every point is a single variant determined by the frequencies in each of the samples involved in the comparison, and calculated R2 as a measure of correlation (Supplementary Fig. 5). Polyclonal infections were Prediction of WGS-DST profiles and comparison with culture DST. Using the almost always associated with low values of this index, except in cases where they SNPs obtained from the analysis of genomic data, we performed an antibiotic were very subtle. We then proceeded to expand on the analysis generating a density resistance prediction for every sample. For this, reliable catalogs of resistance- graph that plots the distribution of SNP differences between the profiles. That is, in associated variants were used, namely PhyResSE56 and ReSeqTB57 databases. We which range of frequencies is the comparison more enriched. This procedure also considered as likely resistance variants any small INDEL present in genes allowed us to better define clonal infections when the enriched range was at the commonly associated with antibiotic resistance. Once the predictions were very low frequencies, and polyclonal infections when it was in the intermediate obtained, we systematically compared them to the available phenotypic DST results variable frequencies or at high fixed frequencies (Supplementary Fig. 6). to calculate matches and mismatches in the surgery patients dataset. We then computed sensitivity and specificity based on these coincidences and discrepancies. Ethical approval. Ethical approval for the study (#892/01-17) was obtained from the Institutional Review Board (IRB) of the National Center for Tuberculosis and Identification of polyclonal infections. Identification of polyclonal infections Lung Diseases (NCTLD). The study design and conduct complied with all relevant depends on the distribution of the different genotypes among patients samples regulations regarding the use of human study participants and was conducted in (Fig. 1). When two genotypes are in two samples of the same patient, phylogenetic accordance with the criteria set by the Declaration of Helsinki. All enrolled patients and genetic distances can be used to differentiate them. This happens mostly in or their legal guardians provided written informed consent prior to the inclusion in superinfection cases where one strain prevails over the other. In contrast, when the the study. All patient data were de-identified before final data analysis. All study different genotypes are in one sample, as it happens when both genotypes coexist, staff during the study period had up-to-date GCP/GCLP certificates. then deconvolution methods have to be applied to separate them before calculating their genetic distance. Consequently, we apply different methods to identify mul- Reporting summary. Further information on research design is available in the Nature tiple strains depending on the number of samples available from a patient. Research Reporting Summary linked to this article. Identification of polyclonal infections in single samples Data availability Phylogenetic identification. We reasoned that a single sample that shows evidence Sequencing files for all patients in the study are available at Bioproject accession codes of two different genotypes in the same isolate should have two characteristics in a PRJNA480888 and PRJNA318002. Individual accession codes for every sample are phylogeny. First, the sample should show no terminal branch length. Terminal provided in Supplementary Data 1. All relevant data related to this work is available from branch lengths represent private fixed SNPs only present in the sample. When the authors. Additional details about patients are also publicly available at TB Portals60 there is polyclonal infection no private SNPs are seen as any SNP not shared by the (https://data.tbportals.niaid.nih.gov). Resistance prediction databases are available at two genotypes will be at intermediate frequencies. The second rule is that they PhyResSE (https://bioinf.fz-borstel.de/mchips/phyresse) and ReSeqTB (https://platform. should not be part of a transmission cluster. It is known that transmission clusters are enriched in strains that are zero SNPs apart from other strains in the same reseqtb.org). Source data are provided with this paper. cluster. Thus we required the candidate polyclonal infection to be at least 20 SNP apart from another sample in the dataset. We developed a Python script able to Code availability analyze a phylogeny, extract branch lengths using Python’s ete3 toolkit, calculate Custom code used in this study, along with example inputs and outputs are available at phylogenetic distances between all samples, and identify terminal branches with 0 https://gitlab.com/tbgenomicsunit/georgia-polyclonal (https://doi.org/10.5281/ SNPs and more than 20 SNPs apart from any other isolate (Supplementary Fig. 4). zenodo.4604579)61. The code is available under Creative Commons Attribution 4.0 For this, we generated a phylogeny using a neighbor-joining approach and International. Hamming distance, which represents branch length proportional to the absolute number of observed differences. Given the low genetic diversity of MTBC and that most positions are biallelic, Hamming distances reconstruct reasonably well the Received: 5 November 2020; Accepted: 22 March 2021; overall phylogeny with respect to maximum likelihood and it is easier to parse. To show the accuracy and the limit of detection of our approach we also generated in- silico mixes of strains with increasing genetic distance between the pairs selected (5, 10, 15, 20, 25 SNP) and added them to our phylogeny to test the script. Lineage markers identification. To test the coexistence of two strains in one sample, References we checked the appearance of any phylogenetic markers using a database from the 1. Lieberman, T. D. et al. Genomic diversity in autopsy samples reveals within- literature58,59 by means of a custom Python script. For samples that showed evi- host dissemination of HIV-associated Mycobacterium tuberculosis. Nat. Med. dence of more than one marker from different MTBC lineages at a significant 22, 1470–1474 (2016). frequency (>5%), a polyclonal infection was called by this method and the esti- 2. Tarashi, S., Fateh, A., Mirsaeidi, M., Siadat, S. D. & Vaziri, F. Mixed infections mated proportion of the involved genotypes was defined by the approximate fre- in tuberculosis: the missing part in a puzzle. Tuberculosis 107, 168–174 (2017). quencies of their markers. 3. Cadena, A. M. et al. Concurrent infection with Mycobacterium tuberculosis confers robust protection against secondary infection in macaques. PLoS Deconvolution of individual genotypes. For those cases in which we observed a Pathog. 14, e1007305 (2018). polyclonal infection in just one sample by either the phylogenetic or lineage 4. Behr, M. A., Edelstein, P. H. & Ramakrishnan, L. Revisiting the timetable of markers method, we established the proportions of the two strains in the sample tuberculosis. BMJ 362, k2738 (2018). and deconvoluted both the individual genotypes if their difference was big enough 5. Verver, S. et al. Rate of reinfection tuberculosis after successful treatment is (e.g., 80/20 proportions, obtained by phylogenetic markers). By clustering the higher than rate of new tuberculosis. Am. J. Respir. Crit. Care Med. 171, frequencies of all variants matching the phylogenetic markers’ frequencies and 1430–1435 (2005). assigning them to their corresponding genotype, we could isolate both and NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z 6. Merker, M. et al. Compensatory evolution drives multidrug-resistant 36. Vashakidze, S. et al. Retrospective study of clinical and lesion characteristics of tuberculosis in Central Asia. Elife 7, e38200 (2018). patients undergoing surgical treatment for pulmonary tuberculosis in Georgia. 7. Plazzotta, G., Cohen, T. & Colijn, C. Magnitude and sources of bias in the Int. J. Infect. Dis. 56, 200–207 (2017). detection of mixed strain M. tuberculosis infection. J. Theor. Biol. 368, 67–73 37. Siddiqi, S. et al. Direct drug susceptibility testing of mycobacterium (2015). tuberculosis for rapid detection of multidrug resistance using the bactec MGIT 8. Dheda, K. et al. Drug-penetration gradients associated with acquired drug 960 system: a multicenter study. J. Clin. Microbiol. 50, 435–440 (2012). resistance in patients with tuberculosis. Am. J. Respir. Crit. Care Med. 198, 38. Orikiriza, P. et al. Evaluation of the SD Bioline TB Ag MPT64 test for 1208–1219 (2018). identification of Mycobacterium tuberculosis complex from liquid cultures in 9. Dheda, K. et al. Spatial network mapping of pulmonary multidrug-resistant Southwestern Uganda. Afr. J. Lab. Med. 6 (2017). tuberculosis cavities using RNA sequencing. Am. J. Respir. Crit. Care Med. 39. van Soolingen, D., Hermans, P. W., de Haas, P. E., Soll, D. R. & van Embden, 200, 370–380 (2019). J. D. Occurrence and stability of insertion sequences in Mycobacterium 10. Cadena, A. M., Fortune, S. M. & Flynn, J. L. Heterogeneity in tuberculosis. tuberculosis complex strains: evaluation of an insertion sequence-dependent Nat. Rev. Immunol. 17, 691–702 (2017). DNA polymorphism as a tool in the epidemiology of tuberculosis. J. Clin. 11. Lin, P. L. et al. Sterilization of granulomas is common in active and latent Microbiol. 29, 2578–2586 (1991). tuberculosis despite within-host variability in bacterial killing. Nat. Med. 20, 40. Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ 75–79 (2014). preprocessor. Bioinformatics 34, i884–i890 (2018). 12. Ley, S. D., de Vos, M., Van Rie, A. & Warren, R. M. Deciphering within-host 41. Wood, D. E. & Salzberg, S. L. Kraken: ultrafast metagenomic sequence microevolution of through whole-genome sequencing: the phenotypic impact classification using exact alignments. Genome Biol. 15, R46 (2014). and way forward. Microbiol. Mol. Biol. Rev. 83, e00062–18 (2019). 42. Li, H. & Durbin, R. Fast and accurate short read alignment with 13. Cohen, T. et al. Mixed-strain mycobacterium tuberculosis infections and the Burrows–Wheeler transform. Bioinformatics 25, 1754–1760 (2009). implications for tuberculosis treatment and control. Clin. Microbiol. Rev. 25, 43. Comas, I. et al. Human T cell epitopes of Mycobacterium tuberculosis are 708–719 (2012). evolutionarily hyperconserved. Nat. Genet. 42, 498–503 (2010). 14. Andersen, P. & Scriba, T. J. Moving tuberculosis vaccines from theory to 44. Li, H. et al. The sequence alignment/map format and SAMtools. practice. Nat. Rev. Immunol. 19, 550–562 (2019). Bioinformatics 25, 2078–2079 (2009). 15. Sable, S. B., Posey, J. E. & Scriba, T. J. Tuberculosis vaccine development: 45. “Picard Toolkit.” Broad Institute, GitHub Repository. http://broadinstitute. progress in clinical evaluation. Clin. Microbiol. Rev. 33, (2019). github.io/picard/; Broad Institute (2019). 16. World Health Organization. Global tuberculosis report 2018. (World Health 46. Koboldt, D. C. et al. VarScan 2: somatic mutation and copy number alteration Organization, 2018). discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012). 17. Maghradze, N. et al. Classifying recurrent Mycobacterium tuberculosis 47. Auwera, G. A. et al. From FastQ data to high‐confidence variant calls: the cases in Georgia using MIRU-VNTR typing. PLoS ONE 14, e0223610 genome analysis toolkit best practices pipeline. Curr. Protoc. Bioinform. (2019). 43, 11.10.1–11.10.33 (2013). 18. Kempker, R. R., Vashakidze, S., Solomonia, N., Dzidzikashvili, N. & Blumberg, 48. Wilm, A. et al. LoFreq: a sequence-quality aware, ultra-sensitive variant caller H. M. Surgical treatment of drug-resistant tuberculosis. Lancet Infect. Dis. 12, for uncovering cell-population heterogeneity from high-throughput 157–166 (2012). sequencing datasets. Nucleic Acids Res. 40, 11189–11201 (2012). 19. Vashakidze, S. et al. Favorable outcomes for multidrug and extensively drug 49. Cole, S. T. et al. Deciphering the biology of Mycobacterium tuberculosis from resistant tuberculosis patients undergoing surgery. Ann. Thorac. Surg. 95, the complete genome sequence. Nature 393, 537–544 (1998). 1892–1898 (2013). 50. Huang, W., Li, L., Myers, J. R. & Marth, G. T. ART: a next-generation 20. Guerra-Assunção, J. A. et al. Large-scale whole genome sequencing of M. sequencing read simulator. Bioinformatics 28, 593–594 (2012). tuberculosis provides insights into transmission in a high prevalence area. Elife 51. Minh, B. Q. et al. IQ-TREE 2: new models and efficient methods for 4, e05166 (2015). phylogenetic inference in the genomic era. https://doi.org/10.1101/849372. 21. Liu, Q. et al. Mycobacterium tuberculosis clinical isolates carry mutational 52. Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: molecular signatures of host immune environments. Sci. Adv. 6, eaba4901 (2020). evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35, 22. Kreutzer, D. A. & Essigmann, J. M. Oxidized, deaminated cytosines are a 1547–1549 (2018). source of C–T transitions in vivo. Proc. Natl Acad. Sci. USA 95, 3578–3582 53. Huerta-Cepas, J., Serra, F. & Bork, P. ETE 3: reconstruction, analysis, (1998). and visualization of phylogenomic data. Mol. Biol. Evol. 33, 1635–1638 23. Ackley, S. F. et al. Multiple exposures, reinfection and risk of progression to (2016). active tuberculosis. R. Soc. Open Sci. 6, 180999 (2019). 54. Tierney, L. The R statistical computing environment. Lecture Notes in 24. Lee, R. S., Proulx, J.-F., Menzies, D. & Behr, M. A. Progression to tuberculosis Statistics 435–447 (2012) https://doi.org/10.1007/978-1-4614-3520-4_41. disease increases with multiple exposures. Eur. Respir. J. 48, 1682–1689 55. Wilkinson, L. ggplot2: elegant graphics for data analysis by WICKHAM, H. (2016). Biometrics 67, 678–679 (2011). 25. van Rie, A. et al. Reinfection and mixed infection cause changing 56. Feuerriegel, S. et al. PhyResSE: a web tool delineating Mycobacterium Mycobacterium tuberculosis drug-resistance pattern. Am. J. Respir. Crit. Care tuberculosis antibiotic resistance and lineage from whole-genome sequencing Med. 172, 636–642 (2005). data. J. Clin. Microbiol. 53, 1908–1914 (2015). 26. Mollenkopf, H.-J., Kursar, M. & Kaufmann, S. H. E. Immune response to 57. Cirillo, D. M., Miotto, P. & Tagliani, E., ReSeqTB Consortium. Reaching postprimary tuberculosis in mice: Mycobacterium tuberculosis and consensus on drug resistance conferring mutations. Int J. Mycobacteriol. 5, S33 Miycobacterium bovis bacille Calmette–Guérin induce equal protection. J. (2016). Infect. Dis. 190, 588–597 (2004). 58. Coll, F. et al. A robust SNP barcode for typing Mycobacterium tuberculosis 27. Yates, T. A. et al. The transmission of Mycobacterium tuberculosis in high complex strains. Nat. Commun. 5, 4812 (2014). burden settings. Lancet Infect. Dis. 16, 227–238 (2016). 59. Stucki, D. et al. Mycobacterium tuberculosis lineage 4 comprises globally 28. Prideaux, B. et al. The association between sterilizing activity and drug distributed and geographically restricted sublineages. Nat. Genet. 48, distribution into tuberculosis lesions. Nat. Med. 21, 1223–1227 (2015). 1535–1543 (2016). 29. Ordonez, A. A. et al. Dynamic imaging in patients with tuberculosis reveals 60. Rosenthal, A. et al. The TB portals: an open-access, web-based platform for heterogeneous drug exposures in pulmonary lesions. Nat. Med. (2020) https:// global drug-resistant-tuberculosis data sharing and analysis. J. Clin. Microbiol. doi.org/10.1038/s41591-020-0770-2. 55, 3267–3282 (2017). 30. Marakalala, M. J. et al. Inflammatory signaling in human tuberculosis 61. Moreno-Molina M. Genomic analyses of Mycobacterium tuberculosis from granulomas is spatially organized. Nat. Med. 22, 531–538 (2016). human lung resections reveal a high frequency of polyclonal infections. 31. Roca, F. J. & Ramakrishnan, L. TNF dually mediates resistance and Zenodo. https://doi.org/10.5281/zenodo.4604579 (2021). susceptibility to mycobacteria via mitochondrial reactive oxygen species. Cell 153, 521–534 (2013). 32. Payne, J. L. et al. Transition bias influences the evolution of antibiotic resistance in Mycobacterium tuberculosis. PLoS Biol. 17, e3000265 (2019). Acknowledgements 33. Martin, C. J. et al. Digitally barcoding reveals infection dynamics in the The authors were supported by projects SAF2016-77346-R and PID2019-104477RB-I00 macaque model of tuberculosis. MBio 8, e00312–17 (2017). awarded to IC and the grant BES-2017-079656 awarded to MM by the Spanish Ministry 34. Shockey, A. C., Dabney, J. & Pepperell, C. S. Effects of host, sample, and in of Economy and Competitiveness and Ministry of Science, the ERC project 638553-TB- vitro culture on genomic diversity of pathogenic mycobacteria. Front. Genet. ACCELERATE awarded to IC, Spanish Government-FEDER Funds through CV contract 10, 477 (2019). CPII18/00031 and grant PI16/01511, and Generalitat Valencia Grant to I.C. (code 35. Trauner, A. et al. The within-host population dynamics of Mycobacterium PROMETEO/2020/012). The grant providers played no part in study design, data col- tuberculosis vary with treatment efficacy. Genome Biol. 18, 71 (2017). lection, and analysis, or the preparation of the manuscript. 10 NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22705-z ARTICLE Author contributions Reprints and permission information is available at http://www.nature.com/reprints I.C., M.M., S.V., N.S., Z.A., A.R., A.G., C.V., S.G., and R.K. contributed to the study design. N.S., I.K., Z.A., A.R., A.G., and S.V. collected or processed the data. SV performed Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in the surgical resections. N.S., I.K., N.B., M.T., and L.V. performed microbiological or published maps and institutional affiliations. molecular biology work. I.C., M.M., N.S., S.G., R.K., and S.V. analyzed the data, which was interpreted by all authors. I.C., M.M., N.S., S.G., and R.K. drafted the manuscript with input from all authors. All authors have read and approved the final manuscript. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, Competing interests adaptation, distribution and reproduction in any medium or format, as long as you give The authors declare no competing interests. appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless Additional information indicated otherwise in a credit line to the material. If material is not included in the Supplementary information The online version contains supplementary material article’s Creative Commons license and your intended use is not permitted by statutory available at https://doi.org/10.1038/s41467-021-22705-z. regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ Correspondence and requests for materials should be addressed to I.C. licenses/by/4.0/. Peer review information Nature Communications thanks Gregory Fox, Bouke de Jong and the other, anonymous, reviewer for their contribution to the peer review of this © The Author(s) 2021 work. Peer reviewer reports are available. NATURE COMMUNICATIONS | (2021)12:2716 | https://doi.org/10.1038/s41467-021-22705-z | www.nature.com/naturecommunications 11