Genetic Characterization and Population Structure of Drug Resistant Mycobacterium tuberculosis Isolated from Brazilian Patients Using Whole-Genome Sequencing

1. Introduction

Tuberculosis (TB) accounted for 1.4 million of deaths and an estimated 10.6 million of new cases in 2023, and was the leading cause of death from a single infectious agent worldwide until COVID-19 pandemic [1]. The emergence of drug resistant (DR) TB is a global threat that hinders successful TB treatment. TB multidrug resistant (MDR-TB), defined as simultaneous resistance to rifampicin (RIF) and isoniazid (INH), results in a worse prognosis, prolonged TB treatment with second-line drugs that are more toxic, more expensive and with the possibility of evolving to pre-extensive drug resistant (p-XDR-TB). The latter is defined as strains that fulfil the definition of multidrug resistant or rifampicin resistant (MDR/RR-TB) plus resistance to any fluoroquinolone, and subsequently can evolve to extensive drug resistant (XDR-TB), defined as strains p-XDR-TB plus resistance to bedaquiline and/or linezolid.

Brazil is one of the thirty countries with highest TB burden, and although the incidence rate decreased until 2014, during the period of 2015-2019 a notified increased incidence of 34.3 to 37.4 cases/100.000 inhabitants was observed [2]. In addition, due to the COVID-19 pandemic, it is estimated that after a decade of decline, TB mortality has increased in Brazil and globally [3]. An additional concern for TB control in Brazil is the considerable number of DR-TB patients; in 2018 around 2.500 MDR/RR-TB cases were estimated [4], including an increase of MDR-TB among patients that had not previously been treated in Rio de Janeiro State [5].

Rapid DR-TB detection and epidemiological surveillance, as well as knowledge about the genetic diversity of isolates of the Mycobacterium tuberculosis complex (MTBC) in different settings are factors that may contribute to DR-TB elimination. In this scenario, molecular tools have become important and quite recently, next generation sequencing (NGS) made it possible to quickly characterize the whole genome of MTBC strains, enabling both identification of resistance-related genetic variants and lineages identification involved in ongoing transmission [6,7]. Culture-based phenotypic drug susceptibility test (DST), although still the current gold standard, has limitations due to the slow growth rate of the MTBC organisms. Thus, molecular methods for drug resistance prediction are being steadily introduced as a routine in low TB incidence countries [8].

WGS is a promising tool and an approach for DR/MDR-TB detection, since it provides detailed sequence information from different genomic regions, thus enabling drug resistance prediction [9]. However, the high amount of sequencing data generated by WGS created the challenge to develop bioinformatics tools to translate the data into information of clinical and laboratory interest [8]. To permit the use of WGS by professionals with little or no bioinformatics skills, user friendly tools for analysis and interpretation of WGS data have been developed and implemented, permitting accessibility and broad implementation of NGS-based approaches [10,11,12,13]. Nonetheless, due to the complexity of large scale data analyses, some bioinformatics command line skills are still required and sometimes, user a friendly graphical interface is not available [14,15].

Due to global variability of prevalence of MTBC lineages and evidence of differential association with drug resistance [16,17], it is evident that combined detection of both characteristics in different countries and regions may interfere with performance and management of DR-TB detection. Therefore, genomic information of MTBC strains together with conventional phenotypic drug susceptibility testing (DST) and clinical outcome is urgently needed. In addition, evaluation of the different available tools for extraction of DST profile from WGS became important to evaluate the regional differences in DR-TB surveillance and to understand local TB transmission.

Here we conducted a genetic diversity study on genomes from a large collection of Mtb isolates from several states of Brazil that mostly had phenotypic DST for primary and secondary drugs. WGS data were evaluated by an in-house WGS pipeline and different online available pipelines, including Mykrobe Predictor, KvarQ, TB Profiler version 0.3.4 and the more recent TB Profiler 5.0 to predict drug susceptibility and the in-house WGS pipeline, KvarQ, TB Profiler 5.0 and RD-Analyzer to predict the genotype (TB lineages).

2. Results

2.1. Study Samples, Patient Data and Phenotypic Drug Susceptibility Testing

Among 298 Mtb isolates from presumed DR-TB patients with high quality WGS data, 294 had conventional phenotypic DST data for at least one drug. Phenotypic DST resulted in 64 (21.8%) pan-susceptible; 22 (7.5%) isoniazid monoresistant (IMR); nine (3.1%) rifampicin monoresistant (RMR); one (0.3%) ethambutol monoresistant (EMR); 10 (3.4%) poly-resistant (Poly-R); 161 (54.8%) MDR and 27(9.2%) p-XDR (Supplementary Table S1).

Among 170 patients for which we had identification, 112 had sociodemographic, clinical, radiographic and treatment outcome data available (Supplementary Table S2). Among them, 65 (57.5%) were male, seven (6.3%) were smokers, eight (7.1%) were illicit drug users, 16 (14.2%) had diabetes, 12 (10.7%) were alcohol users and seven (6.3%) were HIV infected. In relation to chest radiographic images, 72 patients (64.3%) had bilateral cavitary disease, 15 (13.4%) unilateral and 16 (14.3%) had no image available. Regarding treatment outcome, 50 (44.6%) cured or had treatment completion, 32 (28.6%) died due to TB, 16 (14.3%) were lost to follow-up, five (4.5%) were still under treatment, five (4.5%) died from another cause and four (3.6%) patients had no clinical data available.

2.2. Novel Mutations, In Silico Drug Susceptibility and Resistance Prediction Using Different Pipelines

Considering the in-house WGS pipeline, the 298 sequenced genomes were classified as follows: 80 (26.8%) pan-susceptible; 11 (3.7%) IMR; 4 (1.3%) RMR; 1 (0.3%) EMR; 1 (0.3%) SMR; 10 (3.4%) Poly-R; 124 (41.6%) MDR and 67 (22.5%) p-XDR. The nature and frequencies of the mutations detected with this pipeline are presented in Figure 1 and in Supplementary Table S3.

Using the 294 genomes that had DST available for at least one drug, we compared the performance for drug resistance inference of the five predictive drug resistance pipelines against phenotypic DST (three different pipelines and two versions of TB Profiler). The in-house WGS pipeline showed slightly better overall results when considering sensitivity, specificity, positive and negative predictive values, and accuracy when comparing with Mykrobe Predictor, both versions of TB Profiler and kvarQ (Table 1). For prediction of resistance to RIF and INH, in all tools, the sensitivity was higher than 80% but the specificity was lower than 95%. Regarding the prediction to the other first line drugs, the in-house WGS pipeline showed sensitivity higher than 80% to EMB, while high specificity to STR was observed using the in-house WGS pipeline and Mykrobe Predictor. With respect to second line drugs, highest sensitivity to FQ ofloxacin (OFL) was obtained using TB Profiler, and over 90% specificity to AMK, KAN and CAP was observed in all pipelines. Highest sensitivity to PZA was obtained in the most updated version of TB Profiler (68%).

The agreement among phenotypic DST and the in-house WGS pipeline, KvarQ, Mykrobe Predictor, TB Profiler 0.3.4 and TB Profiler 5.0 was assessed for all drugs considering resistance and susceptibility prediction (regardless mutation type in pipelines). The highest agreement level among the six methods was found in RIF (k = 0.89), and the worst agreement level was in PZA (k = 0.32). An excellent level of agreement was also detected to INH, AMK, KAN, CAP (Supplementary Table S4).

Using TB Profiler 5.0, that is able to detect and report unique mutations, i.e. not yet reported, 20 novel mutations in three resistance genes (katG, pncA and ethA) were identified, which six were present in eight phenotypically resistant isolates (gene-drug related) as shown in Table 2. Unfortunately, ethionamide susceptibility had not been phenotypically characterized in the present study. Mostly isolates with novel mutations were phenotypically MDR (68.7%) and genotypically p-XDR (46.8%). All mutations observed using TB Profiler 5.0 are shown in Supplementary Table S5.

2.3. Genomic Diversity, Phylogenetic Analysis and Lineage Classification Using Different Pipelines

Among the 298 genomes analyzed by the in-house WGS pipeline, 99% belonged to Lineage 4, of which 155 (52%) to sub-lineage 4.3 [LAM], 42 (14.1%) to sub-lineage 4.10 [PGG3], 40 (13.4%) to sub-lineage 4.1.2 [Haarlem], 35 (11.7%) to sub-lineage 4.1.1 [X], 17 (5.7%) to sub-lineage 4.4 [Vietnam], two (0.6%) had not been sub-lineage detected and four (1.3%) demonstrated a mixed classification of L4.3/L4.1.1 (n = 1), L4.3/L4.1.2 (n = 1), L4.3/L4.10 (n = 1) and L4.1.1/L4.4 (n = 1), indicating mixed infections or contaminations; three (1%) were classified as Lineage 1.

A phylogenetic tree (Figure 2) based on 38,563 genome-wide SNPs was constructed using 293 genomes. Five genomes with high count of mixed SNP calls were excluded (four with mixed sub-lineage classification and one L4.1.2 [Haarlem]). The topology was congruent with different lineage and sub-lineages classifications assigned by all SNP-based pipelines as previously described [33,34,35].

Upon comparing SNP-based lineage classification pipelines (in-house WGS pipeline, kvarQ-barcode-Coll14 and TB Profiler) and a RD-based analysis pipeline (RD-Analyzer), the L4.3 [LAM] (and higher resolutions classifications when available, e.g., L4.3.2, L4.3.3, L4.3.4, L4.3.4.1, L4.3.4.2, L4.3.4.2.1) was predominant with all tools. Classification results and lineage proportions are presented in Supplementary Figure S1.

Some discordant genotypic classifications were encountered and described in the Supplementary Figure S2, mainly among RD-based pipeline analyzes and SNP-based pipelines. Nevertheless, the free-marginal kappa coefficient obtained presented an excellent agreement level (k = 0.89 [95%CI: 0.86-0.93]) and the overall agreement was 92.06%.

With respect to the sample origin (Figure 3), 212 (71.2%) were isolates from patients residents from the Southeast region including 173 (58,1%) isolates from Rio de Janeiro and 39 (13.1%) from São Paulo; 46 (15.4%) were from the Midwest region, including 42 (14,1%) from the Distrito Federal, two (0,7%) from Goiás and two (0,7%) from Mato Grosso; 15 (4.9%) were form the northeast including six (2%) from Pernambuco, four (1,3%) from Ceará, three (1%) from Maranhão, one (0,3%) from Piauí and one (0,3%) from Sergipe; nine 9 (3%) were from the South Region including six (2%) from Santa Catarina and three (1%) from Paraná and finally, six (2%) were from the North Region, including three (1%) from Amazonas, two (0,7%) from Tocantins and one (0,3%) from Acre. From eight (2.7%) isolates we had not data about state origin.

2.4. Genomic Clusters Analysis

Among the 170 isolates with patient identification available, 74 (43.5%) belonged to 24 genomic clusters with sizes ranging from 2 to 10 isolates, therefore characterized by a clustering rate of 0.294. The most frequent genotypic resistance profiles in this sampling was MDR-TB (53.5%), followed by p-XDR-TB (32.9%) and a similar distribution was observed among the clustered population, composed by 63.5% MDR-TB and 33.7% p-XDR-TB isolates (Figure 4). Importantly, a significantly association between being part of a cluster and having an MDR-TB genotype was observed (Chi-square = 5.2512; p=0.02).

Upon analysis of genotypes of these 170 isolates using the in-house WGS pipeline, 95 (55.8%) belonged to sub-lineage 4.3 [LAM], 26 (15.2%) were 4.1.2 [Haarlem], 20 (11.7%) were 4.10 [PGG3], 18 (10.5%) were 4.1.1 [X], nine (5.3%) were 4.4 [Vietnam], one (0.6%) was classified as Lineage 1, and one (0.6%) as Lineage 4 without sub-lineage detection. Among the isolates in clusters, L4.3 LAM also was observed most frequently (n = 39; 52.7%).

When concentrating on the origin of these samples, 138 were from Rio de Janeiro (81.1%), six from Pernambuco (3.5%;) and Santa Catarina (3.5%) each, four from Ceará (2.3%), three from Amazonas (1.7%), Maranhão (1.7%) and Paraná (1.7%) each, two from Minas Gerais (1.1%) and Goiás (1.1%) each and one from Tocantins (0.6%), Acre (0.6%) and Mato Grosso (0.6%) each. An even more pronounced amount of the clustered isolates (95.9%) were from Rio de Janeiro while only three isolates were from other states (Ceará [GC8], Goiás and Tocantins [GC15]). The latter explains why we observed two genomic clusters with isolates belonging to patients from different states, GC8 with one from Ceará and three from Rio de Janeiro, and GC15 with one from Goiás and one from Tocantins.

2.5. Evolution of Drug Susceptibility Patterns in Patients with Multiple Isolates of Mycobacterium Tuberculosis

Included were 133 genomes from 112 patients from whom clinical and epidemiologic data were available and from which 15 had been followed-up (Table 3 and Supplementary Table S6). These 15 patients had at least two isolates with WGS resulting in a total of 36 genomes: two patients with four isolates, two patients with three isolates and 11 patients with two isolates each. Only two patients presented a change in phenotypical drug resistance: one was susceptible to AMK and changed to resistant without changing the mutational profile (patient P39). Another was phenotypically susceptible to KAN and CAP and changed to resistant (patient P64) but again, without changing the mutational profile.

Regarding the mutational profile changes of those 15 patients, seven patients (two with three isolates each and five with two isolates each) had their genotypic drug resistance profile changed (Table 3). The other eight patients presented no change in the genomes of their respective isolates (genetic distance = 0 SNPs) (Supplementary Table S6).

Emergence of resistance mutations for FQ was observed in five patients in gyrA gene (including P28 = S91P [T→C], P27 = A90V [C→T], P63 = D94H [G→C], P10 = D94A [A→C], P39 = D94G [A→G] / A90V [C→T]); one patient cured from disease while the other four died of tuberculosis. Two patients had an EMB mutation resistance emergence in embB gene (P54 = G406D [G→A] and P35 = Q497R [A→G]); the first was cured while the second patient abandoned treatment with no further information obtained.

2.6. Treatment Outcome, Risk Factors and Lineage Associations with Genotypic and Phenotypic Resistance

The association of TB treatment outcome and genotypic and phenotypic drug resistance was evaluated. Genotypic p-XDR-TB was significantly associated with mortality (Chi-square = 4.6231; p = 0.03). Both phenotypic resistance (Chi-square = 7.3445; p = 0.006) and genotypic resistance to EMB (Fisher’s exact test p = 0.0007) were significantly associated with mortality. Another significant association was EMB genotypic resistance (Chi-square = 24.364; p = 0.00001) with the p-XDR-TB genotype. All isolates with phenotypic EMB available are summarized in Supplementary Table S7. In short, the table shown EMB test available isolates divided in four groups: EMB genotypically susceptible and resistant and phenotypically susceptible and resistant; each of these groups was also divided in two groups: cured and deceased.

Comparing patients who had any risk factors to develop TB (HIV, diabetes, alcoholism, smoking and/or illicit drugs abuse) (Supplementary Table S2) among the two outcome groups, deceased and cured, a significant association between presence of any risk factor and deceased outcome was observed (Chi-square = 5.7672; p = 0.01).

Upon comparing of the frequency of the presence drug resistance related mutations in Mtb Isolates classified as L4.3 (LAM) when compared to those of other sub-lineages, a statistically significant lower frequency of mutation related with resistance to STR (Chi-square = 5.3496; p = 0.02), CAP (Chi-square = 6.4498; p = 0.01), ETH (Chi-square = 5.7636; p = 0.01) and AMK (Fisher’s exact test p = 0.001) were associated with LAM sub-lineage in comparison to non-LAM isolates.

3. Discussion

In this study we performed a genome-based characterization of the genetic diversity of mainly drug resistant isolates of Mtb in different regions in Brazil. The analyzed Mtb strains predominantly belonged to the sub-lineage 4.3 (LAM) (52%). This is in accordance with previous studies performed in various Brazilian regions in both susceptible and DR-TB samples [18,19,20,21,22,23,24]. L4.3 was predominant in the South, North, Southeast and Midwest region. In the Northeast, this predominance by L4.3 was not observed but the number of samples from this region was low. Three isolates (from São Paulo, Rio de Janeiro and Distrito Federal) belonged to Lineage 1, known to be restricted mostly to Eastern Africa and South and East of Asia [16], but apparently described with a certain frequency in the Northern Brazil region [24,25,26], probably imported trough slave trade from East Africa [27]. Although L1 is usually not associated with DR, two of three isolates in this study presented some resistance (one IMR and one Poly-R). This might be due to sampling bias as the CRPHF mainly receives samples from patients suspicious for being DR, either due to treatment failure, treatment abandonment, relapse TB or contact with TB resistant patients.

The four pipelines presented a high agreement level (k = 0.89), with almost complete agreement among SNP-based pipelines, but showing divergences mainly in lineage classifications by RD-Analyzer. The Mtb isolates with high counts of mixed SNP calls may be indicative of mixed infection as described by others [28,29,30]. Analysis among five isolates with a high proportion of mixed SNP calls, showed four mixed classifications by the in-house WGS pipeline only, while TB Profiler and KvarQ detected only one of the two lineages present in these isolates. The RD-Analyzer had some difficulties to classify lineage in comparison with SNP-based tools and presented the most divergences among all tools tested, including 31 isolates with mixed classification and six unidentified despite good genome quality and low mixed SNP call counts.

The majority of clusters were identified in Rio de Janeiro, and although sampling from this state was clearly over represented in our setting, this may represent higher levels of recent infection and drug resistant TB transmission, differing of TB cases due to reactivation of latent infection [31,32]. In this Brazilian State, an increase of primary MDR-TB has been described as associated with low TB control performance [5]. Studies in other Brazilian States also shown important ongoing transmission of MDR-TB, as observed in Sao Paulo [33,34] and Santa Catarina [35] and also, in a genome-based study [36], in Rio Grande do Sul. In the latter, as well as in our study, a significant association between clustering and genotypic MDR-TB and significant ongoing transmission in p-XDR-TB resistance profiles was observed, with practically half of the patients infected with genotypic p-XDR-TB (44.6%; n = 25) in clusters, reflecting an alarming scenario of insufficient DR-TB control [36].

In addition to clusters of isolates derived from residents of particular states, we also observed clusters suggestive of interstate transmission of particular genotypes that seem to have been circulating in the country for a considerable time. Examples are GC8 (three patients from Rio de Janeiro and one from Ceará) and GC15 (one from Goiás and one from Tocantins) presenting genetic distances ≤12 SNPs. Furthermore, we observed among isolates of unidentified patients, that were not included in clustering analysis, one isolate from Distrito Federal (isolated in the year 2007) that had genetic distance of ≤12 SNPs to the GC4 isolates that were from Rio de Janeiro (2008-2012), and also one isolate from São Paulo (2004) and one isolate from Distrito Federal (2007) that had genetic distance of 11 SNPs between them. These could point to other two interstate transmission events, requiring further investigation.

From 15 patients, we had several isolates collected during and after treatment. In most of these cases, the infecting lineage strain during follow up of the patients was the same and only in one patient (P35), reinfection with another strain and lineage was observed. In a retrospective cohort study from China [37], a MIRU-VNTR-based recurrence definition observed a lower frequency of reinfection (n = 21; 36.2%) in comparison to relapse (n = 37; 63.8%). In a population-based study from Malawi [38] using WGS, lower frequency of reinfection (n = 20; 14.3%) compared with relapse (n = 55; 39.5%) was also observed; 64 (46%) of recurrent patients remained unclassified, showing how challenging it remains to classify such cases.

We observed emergence of FQ resistance conferring mutations among one third of the patients that were followed-up during treatment. These bacterial populations had likely been selected due non-compliance of drug treatment [39,40]. In addition, emergence of phenotypic resistance to AMK in patient P39, and phenotypic resistance to KAN and CAP in patient P64 was observed, but without detection of known mutations in rrs and eis. This could be caused by not yet known resistant conferring mutations to those drugs and/or by efflux pump activity, not investigated here. Such pumps have been described to be induced by sub-inhibitory levels of antibiotics such as FQ, RIF and INH in inappropriate treatment regimens [41]. Additionally, the presence of a minor undetected population of drug resistant bacteria, so-called heteroresistance, can also explain differences between culture and in silico based drug susceptibility outcome [42].

Another observation was the emergence of the G406D mutation in embB in patient P54 (Table 3), however, a mixed call was detected in this codon and may indicate selection of the mutant population within a single strain, probably due to inappropriate treatment. Patient P35 presented emergence of the Q497R mutation in embB, and in addition a shift of RIF resistance genotype, initially presenting Q432P mutation in rpoB and posteriorly the S450L mutation. Besides that, in the first isolate the sub-lineage detected was L4.1.1 (X) and in the second was L4.3 (LAM), and also a genomic difference of 484 SNPs was observed between both, indicating a TB recurrence by reinfection [38].

We observed a statistically significant lower frequency of resistance related mutations to STR, CAP, ETH and AMK in strains belonging to sub-lineage 4.3 (LAM) in comparison with non-LAM sub-lineages. Lower frequency of mutations related to resistance in some lineages has been described, including lower frequency of S315T mutation in katG in LAM spoligofamily when compared to Haarlem [43]. Another study showed a higher mutation rate in Mtb strains of lineage 2 and as such, a faster emergence of resistance conferring mutations in strains of lineage 2, accompanied by development of DR-TB [44].

Upon comparing among patient treatment outcome, risk factors for TB development, phenotypic and genotypic resistance and Mtb lineage, we observed a significant association of EMB phenotypic (Chi-square p = 0.006) and genotypic (Fisher’s exact test p = 0.0007) resistance with TB mortality. EMB resistance has been described as a risk factor for mortality, mainly when associated with MDR pattern [43]. In our study, genotypic p-XDR-TB was significantly associated with EMB genotypic resistance, while genotypic p-XDR-TB profile was significantly associated with higher morbidity (Chi-square p = 0.03). In addition, the presence of any of the five risk factors (HIV, diabetes, smoking cigarette, alcohol and/or illicit drug use) was also significantly associated with morbidity but independent of genotypic or phenotypic EMB resistance. The increased risk of mortality in MDR-TB and p-XDR-TB patients has been described [46,47] and is related with increasing level of drug resistances, whereas accumulating of even more factors increases the probability of an unfavorable outcome [48,49]; EMB resistance could be an important stage of resistance accumulation in Mtb.

We observed a similar proportion of main drug resistance conferring mutations, frequencies and genes affected detected by all pipelines tested and this is in agreement with other studies [50,51,52,53,54].

Unprecedented was the detection of 20 novel mutations in 32 strains using TB Profiler 5.0, including a novel frameshift mutation katG_c.2070delC in an isolate that is likely associated with INH resistance, considering the impact this may have on the catalase-peroxidase structure and activity, accompanied by phenotypic resistance against INH. The same isolate also presented the mutation ahpC_c.-81C>T, classified in WHO catalogue as with uncertain significance, but may be acting in synergy with the newly described to confer resistance; this needs to be better investigated. This isolate was misclassified as genetically susceptible to INH by all pipelines but TB Profiler, and the same was observed in other isolates with the mutation ahpC_c.-81C>T, highlighting the importance of using updated pipelines.

Only four of 32 isolates were fully phenotypically characterized for the nine drugs available in this study. Unfortunately, thirteen of the isolates in which we detected a novel mutation were not phenotypically characterized for the drug of interest, making it impossible to establish the genetic and phenotypic correlation. Among the four isolates with the frameshift mutation pncA_c.193_200dupTCCTCGTC, two had no PZA DST available and two were PZA resistant (these two were from the same patient); the mutation pncA_c.305dupC was observed in one resistant and one with DST not available; finally, the mutations pncA_c.452dupT and pncA_c.75_79delCGCGC were both presents in resistant isolates, all suggesting resistance association. On the other hand, the frameshift mutation pncA_c.443_444dupGC was detected in both susceptible and resistant isolates, making it difficult to interpret. Other isolates observed with novel frameshift mutations in pncA and katG genes and phenotypically susceptible could be explained by an acquired low-level resistance, under the antibiotic concentration threshold used here, or due to some laboratorial error in phenotypic resistance determination. Phenotypic susceptible isolates harboring confident resistance mutations have been reported in other studies [55,56].

Among ten other mutations we detected that were absent in the WHO catalogue (Supplementary Table S5), eight [57,58,59,60,61,62,63,64] in twelve isolates had not yet reported in Brazil [63,65], four (katG_p.Leu634Phe, pncA_p.Asp63His, pncA_p.Gly23Val, pncA_c.521_522insT) were present in phenotypically resistant isolates in our study. Another mutation associated with resistance was detected for the first time in Brazil, rpoB_p.Ile491Phe, a mutation outside of the 81 bp hotspot of the rpoB gene, what make it undetectable by commercial assays such as GeneXpert MTB/RIF and MTBDRplus [66,67]. In Eswatini, for example, the presence of this mutation is a significant problem because of its prevalence of >60% in MDR strains [68]. In our study, this mutation was observed in one isolate genotypically and phenotypically MDR, that also carried the rpoB_p.Ser493Leu mutation, an uncertain variant according to WHO. However, their combination can be acting in synergy because the former is a borderline resistance mutation [69].

The comparison among the pipelines evaluating their performance through sensitivity, specificity, accuracy, positive and negative predictive values calculation for drug resistance prediction, using phenotypic resistance results as a reference, showed a better overall performance of resistance prediction by the in-house WGS pipeline. For all pipelines, the sensitivity of prediction of resistance to RIF and INH was higher than 80% but the specificity was lower than 95% and therefore below that of the WHO recommendations [70]. High positive predictive value-PPV (>80%) was observed in all pipelines used for detection of resistance conferring mutations for RIF, INH, OFL, AMK and similar to that described by others [71,72]. Lower PPV was observed for detection of resistance mutations whose action mechanisms are more complex and genetic bases of resistance is less understood, such as resistance to PZA, KAN and CAP. Unreliable PZA resistance prediction has been reported in various studies [25,71,72,73] and in order to enhance PZA resistance prediction, the following strategies have been described: (i) improving the mutation library of the pipelines mainly by including detection of indels in pncA, (ii) the “non-wild type sequence” approach, that consist in interpret any non-synonymous mutations or indels in pncA as genotypic PZA-resistant strains, (iii) besides manually checking of the sequence reads [74,75]. Not surprisingly, the best sensitivity performance for PZA resistance prediction was achieved by TB Profiler 5.0, the most up-to-date pipeline for detection of DR-associated mutations, accounting with the largest PZA-mutations catalogue among all pipelines used and able to detect novel mutations, including fourteen novel indels observed in this study.

Another important finding in our study is that, compared to conventional DST, more than twice as many cases were classified as p-XDR-TB using WGS. This difference in proportion could be beyond the accuracy of genomic tests to second-line resistance detection. However, one should take into account that the criteria for performing phenotypic DST tests to second line drugs in Brazil are causing the underestimates of the detection of resistance to such drugs. These results highlight the importance of using WGS for epidemiologic surveillance and control of DR/MDR/p-XDR-TB, as these discrepancies would not be detected and reported without it use.

Without any doubt, a major limitation of this study is that phenotypic tests for second line drugs were not done in all samples. Among the 161 isolates classified as phenotipically MDR-TB only 22 were tested to all second line drugs and only 10 were tested for nine drugs; one of 10 Poly-R-TB and six of the 27 p-XDR-TB were submitted to the full DST. This further emphasizes the importance of using molecular methods for TB diagnosis and comprehensive resistance analyzes. Another limitation is the small sample size from regions outside Rio de Janeiro, Distrito Federal and Sao Paulo States, therefore not representing the national scenario of DR-TB in Brazil.

In conclusion, the evaluated pipelines for prediction of MTBC lineage and drug resistance work well in the Brazilian sampling studied here and our data favor the use of WGS in cases without conventional DST. Although the in-house WGS pipeline performed slightly better in general, all tools performed well in predicting DR for RIF, INH, AMK, KAN and CAP. Importantly they allowed detection of p-XDR-TB strains that otherwise would probably have been unreported. Analysis of isolates of DR-TB patients that were sampled sometimes years apart, demonstrated that several accumulated drug resistance mutations showing that resistance evolution occurs by acquisition of mutations in the same strain, probably due to noncompliance of treatment. Phylogenetic analysis and lineage characterization contribute to better understand the MTBC genotypes spreading in the population, so WGS improves the knowledge of MTBC population structure and evolution and offers rapid and reliable assessment of resistance related mutations allowing faster access to effective treatment. A surprisingly low level of reinfection was observed in an area with high TB incidence and, even in patients with DR and prolonged treatment. We believe that these findings highlight the importance of the need for active surveillance throughout the national territory in order to avoid further aggravation of the TB scenario in Brazil.

4. Materials and Methods

4.1. Sample Collections and Phenotypic Drug Resistance

We used a convenience sample for this study of 298 Mtb clinical isolates available at the National Tuberculosis Reference Center Professor Hélio Fraga (CRPHF; Rio de Janeiro, Brazil). The sample set included isolates of patients from 16 different Brazilian States with presumed DR-TB due to a history of treatment failure, treatment abandonment, relapse TB or contact with DR-TB patients. The isolates were therefore submitted to phenotypic drug susceptibility testing (DST).

DST was performed to the first line drugs RIF, INH, PZA, EMB and to the second line drugs STR, OFL, AMK KAN and CAP, using the liquid MB/BacT system (Organon Teknika Corp., Durham, NC, USA). The following cut-off values were used: RIF (1.0 mg/L), INH (0.1 mg/L), EMB (5.0 mg/L), PZA (100.0 mg/L), STR (1.0 mg/L), OFX (2.0 mg/L), AMK (1.0 mg/L), KAN (4.0 mg/L) and CAP (2.5 mg/L).

4.2. Whole Genome Sequencing and in-House Pipeline

The Mtb genomic DNA was extracted as described previously [76]. To perform WGS, sequencing libraries were prepared as described previously [77], and sequenced on an Illumina platform at Division of Infectious Diseases and Environmental Health of the Norwegian Institute of Public Health. To describe the MTBC genetic diversity and genotypic resistance profile, a Single Nucleotide Polymorphism (SNP)-based Mtb lineage and sub-lineage, and resistance related genes (11 coding genes and four promoter regions) was assigned first using an in-house pipeline developed at the Swiss Tropical and Public Health Institute (Swiss TPH) as described before [78], hereafter called in-house WGS pipeline. In brief, FASTQ files were processed with Trimmmomatic v 0.33 (SLIDINGWINDOW:5:20) [79] to remove Illumina adaptors and trim low quality reads. Overlapping reads were then merged with SeqPrep v 1.2. Duplicated reads were marked by the MarkDuplicates module of Picard v 2.9.1. The resulting reads were mapped to a reconstructed ancestral sequence of MTBC as described previously [80] using the BWA-MEM v 0.7.13 algorithm. Pysam v 0.9.0 was used to exclude reads with alignment score lower than (0.93*read_length)-(read_length*4*0.07)): this corresponds to more than seven miss-matches per 100 bp. SNPs were called with Samtools v 1.2 mpileup and VarScan v 2.4.1. SNPs annotation used was that of Mycobacterium tuberculosis H37Rv reference strain (NC_000962.3). Repetitive regions of the genome as PE/PPE/PGRS were excluded as described before [81].

4.3. Phylogenetic Analysis

A maximum likelihood phylogenetic tree using 100 bootstraps was constructed using an alignment containing all variable positions from all quality filtered genomes with ≥15x average coverage using MEGAX v 10.2 [82], and the resulting tree was rooted using M. canettii (Genbank accession number: NC_019950.1). The tree was visualized using iTOL v 6.7 [83].

4.4. Lineage Classification Pipelines Comparisons

Mtb lineages classification by the in-house WGS pipeline was compared with three online available pipelines: TB Profiler version 5.0 [84], KvarQ version 0.12.2 [10] but using the SNP-barcode used elsewhere [11], and RD-Analyzer version 1.01 - a region of difference based analyzer of Mtb from sequence reads [85].

4.5. Genotypic Resistance Detection and Pipelines Evaluation

In order to identify the drug resistance related variants among the included Mtb genomes, we used the in-house WGS pipeline and the command line versions of the pipelines Mykrobe Predictor version 0.3.3 [12], TB Profiler version 0.3.4 [73], KvarQ version 0.12.2 [10] and TB Profiler version 5.0 (https://github.com/jodyphelan/TBProfiler), a recent pipeline version that accounts for the updated WHO catalogue of mutations conferring resistance (https://www.who.int/publications/i/item/9789240082410).

4.6. Cluster Analysis and Recent Transmission

Genomic clusters were delineated using the tool snp-dists v0.7.0 (https://github.com/tseemann/snp-dists). A SNP-based distance matrix was built using the alignment of 170 high quality genomes which not presented mixed calls. A genomic cluster was defined when two or more isolates presented genomes with ≤12 SNPs of difference, which encompasses not only linked cases such as household contacts (genetic distance from zero to five SNPs), but also related cases (five to 12 SNPs). These thresholds are considered as representing very likely recent transmissions event, especially in settings using isolates from chronically infected patients and DR-TB [86,87,88]. To estimate recent transmission, 170 genomes, corresponding to one genome per patient appropriately identified (patient ID), were used. Samples without patient identification and/or duplicated ID were excluded of this analysis. The reconstructed maximum likelihood tree using 100 bootstraps was visualized in GrapeTree [89]. Clustering rate was calculated using the following formula (nc – c)/n, where nc = the total isolates in cluster, c = the number of clusters and n = the total number of isolates.

4.7. Patient data

Patient clinical data available, included TB treatment outcome, diabetes, alcohol use, illicit drug, cigarette smoking, HIV coinfection status and sex, were obtained from the "Sistema de Informação de Agravos de Notificação" (SINAN) and "Sistema de Informação de Tratamentos Especiais de Tuberculose" (SITETB), two national databases for disease surveillance. The project was approved by Research Ethics Committee of Federal University of Rio de Janeiro (CAEE 10126919.2.0000.5257).

4.8. Statistical Analyses

For prediction of drug resistance, we calculated sensitivity, specificity, accuracy and positive and negative predictive values using the statistical software R (Version 3.6.1). First and second-line anti-TB drugs prediction were evaluated for the four pipelines, using MB/BacT DST as the reference standard.

Free-marginal Kappa Randolph’s statistics [90] with a 95% confidence interval was used to determine agreement among Mtb lineage classification by the four pipelines. The same procedure was followed to determine agreement among phenotypic drug resistance detected by MB/BacT and the genotypic drug resistance detected by the four pipelines, using four levels of agreement: <0.40 (poor), 0.40-0.59 (fair), 0.60-0.80 (moderate/good) and >0.80 (excellent).

Pearson’s chi-square test was used to determine associations among treatment outcome, risk factors, clusters, lineage, genotypic and phenotypic drug resistance using the chisq.test function in the statistical software R (Version 3.6.1). In case the expected value is less than 5, Fisher’s exact test was calculated using the fisher.test function in R.