Unraveling Mycobacterium tuberculosis genomic diversity and evolution in Lisbon, Portugal, a highly drug resistant setting

Genotypic analysis

The 24-loci MIRU-VNTR genotyping technique grouped the MDR-TB isolates into three major clusters: Lisboa3-A, Lisboa3-B and Q1 (Figure 1). Use of the 12-loci set groups Lisboa3-A and -B in a single cluster (Lisboa3, data not shown). Only the Lisboa3-B and Q1 clusters were found to be associated with XDR-TB isolates. Eight of the ten XDR-TB isolates belonged to either Lisboa3-B or Q1 cluster, and one of remaining two strains was found to be Q1-related, raising the possibility of ancestral Q1 XDR followed by posterior divergence from this cluster. No XDR-TB isolate was found to belong to Lisboa3-A cluster.

Genomic analysis

All 56 clinical isolates underwent whole genome sequencing (WGS) The total number of identified SNPs (point mutations differing from H37Rv) ranged between 488–1465 (mean: 928.0, 26.7% in non-coding regions) (Table 1). Of the SNPs on coding regions, 58.5% were considered non-synonymous substitutions yielding a mean non-synonymous/synonymous ratio (N_s/S) of 1.41 (Table 1). AG, CT, GA and TC transitions were found to be the most frequent substitution types (see Additional file 1), which is reflected by a mean transversion/transition ratio (T_v/T_s) of 0.62. Overall, across the 56 clinical isolates and 19 publicly available reference strains (F11, CDC1551, KZN1435, KZN4207, KZN605, KZN_R506, KZN_V2475, UT205, RGTB327, RGTB423, CCDC5180, CCDC5079, CTRI-2, BTB05_552, BTB05_559, S96_129, HN878, R1207, and X122), 9419 genome-wide SNPs were identified by mapping to the reference genome of M. tuberculosis H37Rv. The number of small insertions and deletions (indels) detected upon read mapping ranged between 15–175 indels per isolate with a size between 1–59 bp (Table 1).

Global phylogenetic analysis using WGS

Using WGS data, the 56 clinical isolates and 19 publicly available strains were assigned into established six SNP Clusters Groups (SCG) and three Principal Genetic Groups (PCG) [14, 20]. Overall, at least one isolate belonging to each SCG and subgroups was included in the subsequent analysis. Forty-four (78.6%) of the 56 clinical isolates belonged to SCG 5, reflecting the high prevalence of these strains in Lisbon Health Region (Table 1).

A phylogenetic tree was inferred from a set of 9419 genome-wide SNPs (Figure 2). It reveals that the two main genetic clusters associated with XDR-TB in the region, Q1 and Lisboa3, constitute two genetically close but distinct clades within the SCG 5. The MIRU-VNTR Lisboa3-A cluster was found to form a monophyletic group within the Lisboa3 clade. The MIRU-VNTR Lisboa3-B clade designation was therefore considered as paraphyletic in the light of a genome-wide SNP phylogeny. The sequenced strain closest to the Lisboa3-Q1 clade is M. tuberculosis UT205, a virulent Colombian isolate that according with the present phylogeny shares a more recent common ancestor with Q1 strains than these do with Lisboa3 strains.

Global evolution through large sequence polymorphisms

Genomes in the M. tuberculosis complex can downsize, through Large Sequence Polymorphisms (LSP) or Regions of Difference (RD), and 89 have been previously identified [21–23]. Across the 75 isolates, 29 (of 89) were detected as absent in at least one isolate (see Additional file 2). The most prevalent RDs detected were RD149 (64 isolates), RD152 (45), RD174 (43), RD3 (64), RD6 (54) and RD^RIO (43). As expected, all 43 strains bearing the RD174 (LAM) deletion also had the RD^RIO deletion [24]. Both deletions constitute a distinct sub-lineage within the Euro-American lineage [23] and were detected only among SCG 5 strains. UT205, like the Q1 and Lisboa 3 samples, had both deletions, confirming its phylogenetic proximity with these M/XDR associated strains (see Additional file 2).

All nine isolates from the SCG 2 had the RD105 deletion characteristic of the East-Asian clade. Of these, other RD deletions were present (RD207 9 isolates, RD181 8, RD142 2). Moreover, other RD deletions associated with specific lineages were detected: RD750 (East-African-Indian lineage, SCG 3a, 1 isolate), RD115 (Euro-American lineage, Americas-Europe sublineage, SCG 5, 8), RD183 (Euro-American lineage, Americas-Europe sublineage, SCG 3c, 1) RD193 (Euro-American lineage, Americas-Europe sublineage, SCG 4, 3), RD219 (Euro-American lineage, Americas-Europe sublineage, SCG 6a, 3), RD761 (Euro-American lineage, South Africa sublineage, SCG 5, 1 (F11 strain)) and RD724 (Euro-American lineage, Central Africa sublineage, SCG 5, 3).

No RD region was found to be absent in RGTB423 and only RD^RIO deletion was detected in RGTB327. Strain RGTB423 has been found to belong to SCG 1 and PGG 1 [25], but in-silico PCR analysis showed that the strain had the pks15/1 7 bp frameshift deletion and the TbD1 deletion indicative of a modern Euro-American strain [23]. Nevertheless, this classification is incongruent with the SCG and PGG classification [11]. On the other hand, RGTB327 was found to have the RD^RIO deletion only and in silico PCR of the pks15/1 and TbD1 loci also pointed towards a modern Euro-American strain, despite the fact that deletion RD174 was not detected. Further sequencing of these two assembled strains may be required to resolve incongruences.

Other structural variants were also searched employing different methodologies (see Additional files 3, 4, 5, 6, 7 and 8). From the deletions described in the supplementary material we highlight characteristic deletions for the Lisboa3 clade and the ARS6559 isolate (complete Lisboa3 subtree in Figure 3) (112 bp, position 2727803, PPE41 gene), as well as the Q1 clade strains (297 bp, position 3929891, ORF PE_PGRS53). Both deletions were also validated by mapping coverage, but further laboratory confirmation is required.

Microevolution towards multidrug and extensively drug resistance

Given the relative high number of sequenced strains present in both Lisboa3 and Q1 clades it was possible to trace the microevolutionary path reflecting the genomic changes accompanying the resistance acquisition process. We considered the subtrees containing the Lisboa3 and Q1 clades plus one or two strains for the Lisboa3 and Q1 subtrees, respectively, included as outgroups for the ensuing analysis (Figure 3). In particular, we inferred the changes in candidate resistant gene mutations at the nodes of the trees.

The Lisboa3 subtree, including the outgroup strain ARS6559, was found to be characterized by a 5 bp deletion on the iniA gene. There is a common acquisition of high-level INH resistance through a inhA double mutation (in node B). The data also reflect the acquisition of RIF resistance in three separate occasions, twice in the Lisboa3-B strains by a rpoB S450L (equivalent to E. coli S531L) and in Lisboa3-A lineage by a rpoB D435V (equivalent in E. coli to D516V). Acquisition of XDR can be seen in the two branches: the first by acquisition of an eis G-10A, gyrA S91P and tlyA Ins755GT mutations (node B1); and, by an eis G-10A and gyrA D94G mutations (node E1). EMB resistance is likely to have been acquired twice by embB M306V and P397T mutations. The latter mutation has been previously reported in one EMB resistant isolate [26]. PZA resistance was found to be acquired on multiple independent occasions through pncA mutations.

The Q1 subtree included two other Q1-related strains as outgroups. Here, it is possible to distinguish the acquisition of INH low-level resistance by an inhA C-15 T mutation (node B) from the acquisition of a higher INH resistance level by an inhA missense mutation (I194A, node C) [27]. Some of the isolates present in the subtree were found outside the Q1 MIRU-VNTR cluster, but share more recent common ancestors with other strains in the clade, potentially indicating subsequent MIRU-VNTR divergence. The Q1 clade has, therefore, been defined as all isolates bearing the gidB A80P mutation characteristic of this cluster and associated with STP intermediate-level resistance previously described by some of us [10]. A more linear resistance acquisition dynamic was found for this clade. EMB resistance was acquired on two possible occasions, through an embB M423T (node C) and M306V (node D) mutations. RIF resistance development, leading to MDR-TB, was found to be acquired by a rpoB S450L mutation (node D), although a second mutation on rpoB (L731P) was later developed (node E). Resistance to PZA, injectable second-line drugs and FQs occurred once by mutations on pncA (V125G, node D), rrs (A1401G, node F) and gyrA (D94A, node G), respectively. Interestingly, isolates IHMT308_08 and IHMT361_08 did not show the two latter mutations in rrs and gyrA genes, and therefore inconsistent with both strains positioning in the Q1 subtree.

A further observation is that M/XDR development in the Lisboa3 subtree appeared to be accompanied by a higher genomic diversification, translated in the number of SNPs and small indels (Additional files 9 and 10). This observation is probably in line with an earlier emergence of the Lisboa3 clade and prolonged circulation in the community leading to a higher intra-clade diversity when compared to Q1 strains. Moreover, isolates from the Lisboa3 and Q1 clades were found to bear a mean proportion of 0.73% (range: 0.2-1.8%) and 0.85% (range: 0.2-1.6%) unique SNPs, respectively, in comparison with the total SNP count of each strain. Both clades were found to share a pool of 654 (67.7-90.3%) and 626 (68.2-85.2%) common SNPs, respectively (Additional file 11). This intra-cluster degree of genomic uniqueness is comparable with the data reported by Niemann et al. for the comparison of two Beijing isolates from the same outbreak clone [13].

Mutational compensation for RIF-resistance

We also screened the remaining RNA polymerase subunits, RpoB and RpoZ, but only eight non-synonymous mutations were identified in RpoB, concomitantly with other RIF resistance associated mutations in RpoB (Table 2). Five RpoB mutations (P45A, T328N, L452P, L731P and I1106T) were predicted to affect protein function after SIFT score analysis (SIFT score <0.05) (Table 2).

Insertion sequence (IS) mapping and functional consequences for genomic stability

Transposition events from ISs can have a profound effect on strain physiology given the possibility of interference with gene expression by ORF knock-out or gene upregulation resulting from upstream transposition [32, 33]. For all strains included in the phylogenetic analysis, we attempted to map the site of all ISs annotated as mobile elements in the genome of M. tuberculosis H37Rv, namely IS6110. Some complex inversions were found to be predominantly transpositional events from multi-copy mobile-elements, such as IS6110. The analysis revealed the presence of IS6110, IS1081, IS1547, IS1557 and IS1558 in multiple copies, but differing in size or annotated sequence at both extremities. For this reason these ISs have been excluded from the mapping analysis.

Variability was only observed for IS1561 and IS1532 (Additional file 15). As expected, IS1561 was not detected in all isolates bearing the RD^RIO deletion, whereas IS1532 is absent in isolates bearing the RD6 deletion found on different SCGs. For IS6110, a total of 251 candidate insertion sites have been obtained (Additional file 16), classified as of high (160), medium (18) or lesser (73) confidence. Almost half (125 (49.8%)) of the 251 ISs were observed on the positive strand. A total of 105 (41.8%) insertion sites were found to be intergenic, from which 64 (25.5%) were in the same orientation with an upstream ORF, known to exert a putative upregulatory effect. For these latter insertion sites the distance from the 3′ end to the upstream ORF ranged between 0–939 bp (47 (18.7%) less than 300 bp). Thirty-three sites were found to be within PE/PPE genes, while three other insertion sites were located 18–38 bp upstream of a PPE gene.

Lisboa3 and Q1 clades were found to share 7 IS6110 sites but were differentiated by IS6110 insertions on positions 889015 (intergenic) and 4183431 (Rv3732 knock-out) for Lisboa3 and, on 2582457 (intergenic) for Q1 isolates (Figure 4). Moreover, we have found that strains belonging to Lisboa3-A MIRU-VNTR cluster (rpoB D435V clade on Figure 3-A) share three distinct IS6110 insertion sites on Rv1682 (position 1906425), Rv2818c (position 3125900) and Rv3096 (position 3465467). Strains from the XDR-TB Lisboa3 B1 clade (Figure 3-A) share a distinctive IS6110 site on the plcC gene (position 2628462). Although no common IS6110 site was found for the SCG 5 strains, SCG 2 strains were found to share three IS6110 sites: an intergenic site on position 888786; on Rv1754c (position 1986639); and, on Rv2820c (position 3127931). SCG 4 strains were found to also share three IS6110 sites on mmpS1 (position 483580), PPE46 (position 3377326) and PPE47 (position 3379768). One hundred and fifty-three (60.0%) sites were found to be specific to a single isolate.

Interestingly, an IS6110 insertion in the NTF locus (position 3493907) was detected in six out of the eight Beijing strains included in the analysis, which is a characteristic of the Beijing/W family [34, 35] (Additional file 16). Hence, two of the three Beijing isolates recovered in Lisbon Health Region were found to belong to the Beijing/W family. No relation with the New York City Beijing/W MDR clade was found as a second insertion in the NTF locus was not detected in any strains [34, 35]. Curiously, a SCG 6a strain (HPV113_08) shared the latter insertion site with the Beijing/W strains, although only one end was detected which can be indicative of another genomic rearrangement. A SCG5 strain (HPV157_06) was found to have an IS6110 67 bp upstream of the characteristic IS6110 insertion site of the Beijing/W family, however in a different orientation. Both insertion sites are found within Rv3128c. This latter gene has an in-frame amber nonsense mutation in H37Rv and for this reason any functional consequence of IS6110-mediated ORF disruption is highly questionable.

Strains belonging to PGG2 were found to have a significantly lower number of IS6110 copies when compared with PGG1 strains (Kruskal-Wallis test, p <0.001). Given the reduced number of PGG3 strains no statistical comparison was possible to perform.

Differential substitution ratios highlight different genomic adaptation strategies

A statistically significant difference in the N_s/S ratio was observed between Lisboa3/Q1 and Beijing strains and others, but only the Lisboa3 and Q1 result met a multiple comparison threshold (Additional file 17). The only significant T_v/T_s ratio difference occurred for differences between Lisboa3 and Q1 clusters (Q1 greater, mean difference: 0.045, p =0.033) (Additional file 17).

These ratios were also found to vary across the genome and across the different Clusters of Orthologous gene Groups (COGs). For each strain, we have computed the N_s/S and T_v/T_s ratio for the different genomic quadrants and for each COG. Overall quadrant N_s/S and T_v/T_s comparison, showed that N_s/S ratio varied along the chromosome such that the second quadrant had a lower N_s/S ratio when compared with the other three quadrants and that the first quadrant had the highest N_s/S mean ratio (Kruskal-Wallis, p <0.001) (Additional file 18). No statistical difference was observed between the third and fourth quadrant. Regarding the T_v/T_s ratio, an approximately inverse situation was found as no statistical difference was observed between the first, third and fourth quadrants. The second quadrant showed however, a significantly higher T_v/T_s ratio than the three other quadrants (Kruskal-Wallis, p <0.001) (Additional file 18).

When these results were stratified by genetic clade, it was found that in the first quadrant the Beijing strains showed a statistically lower N_s/S ratio upon comparison with Q1 and other non-clustered (NC) strains, but not Lisboa3 (Additional file 19). No statistical difference was found in this quadrant for T_v/T_s ratio. In the second quadrant, Lisboa3 strains showed a statistically significant reduced N_s/S ratio compared with the other three groups of strains, while Beijing strains presented a higher N_s/S ratio than the remaining groups (Additional file 19). Inversely, the T_v/T_s ratio on the second quadrant was significantly higher for Beijing strains when compared to Q1 and other NC strains, but not to Lisboa3 strains (Additional file 19). The analysis of the third quadrant showed no statistical difference for Ns/S ratio while Beijing strains showed a higher T_v/T_s ratio on comparison with Lisboa3 and other NC strains, but not Q1 strains. Lisboa3 strains showed a reduced T_v/T_s ratio on this latter quadrant when compared to all other groups. In the fourth quadrant, only a statistical difference was observed for a Q1 reduced N_s/S ratio when comparing with the other strain groups and no significant difference was observed for the T_v/T_s ratio (Additional file 19).

These results show that the N_s/S and T_v/T_s ratio measures appear to vary on a strain and chromosome region dependent mode. Data stratification by isolate and quadrant showed that the T_v/T_s ratio was found to correlate negatively with the N_s/S ratio (Pearson, p <0.001). Correlation between overall isolate N_s/S and T_v/T_s ratio was also attempted but no correlation was found (Pearson, p =0.433).

The comparison of the N_s/S and T_v/T_s ratios across the different COGs also yielded strain dependent results. On comparison with the other three strain groups: Lisboa3 strains showed higher Ns/S ratios on COG groups D (Cell Cycle Control, Mitosis and Meiosis) and P (Inorganic Ion Transport); Q1 strains showed higher Ns/S ratios on COG group V (Defense Mechanisms); and, Beijing strains showed higher N_s/S ratios on COG groups F (Nucleotide Transport and Metabolism), K (Transcription), N (Cell Motility), O (Post translation Modification, Protein turnover and Chaperones) and Q (Secondary Metabolites Biosynthesis, Transport and Catabolism) (Additional file 20). Regarding the T_v/T_s ratio no significant difference was observed for Lisboa3 strains, but higher ratios were observed for Q1 strains in COG groups J (Translation), L (Replication, Recombination and Repair), M (Cell Wall, Membrane Biogenesis) and, for Beijing strains in COG group C (Energy Production and Conversion) (Additional file 21).

These results support the notion of a differential mode of evolution and adaptation to the human host by accumulation/selection of a higher degree of non-synonymous mutations at genes belonging to specific functional categories.

According to recent work by Namouchi et al. [36], the N_s/S ratio varied along the phylogenetic tree, such that terminal branches had a higher N_s/S ratio than inner branches. We have computed the N_s/S and T_v/T_s ratio for the inner nodes assigned in the subtrees in Figure 4 and compared with the respective ratios calculated for the tips of the subtrees. Contrary to the data of Namouchi et al. [36] we have verified that both subtrees had ≈ 6% and ≈ 12% lower Ns/S ratios at the tips of Lisboa3 and Q1 subtrees, respectively, when compared with the inner nodes of the tree (Independent t-test, p <0.001). For the T_v/T_s ratio, the opposite was found: higher T_v/T_s ratios were observed at the tips in comparison with the inner nodes (Mann–Whitney test, p <0.001).

M. tuberculosisgenomic distinctiveness in Lisbon

For at least two decades the Lisbon Health Region in Portugal has been characterized by a high-level of drug resistance, at first MDR-TB, and later XDR-TB, mainly caused by a particular group of strains: the Lisboa family. Presently, this drug resistance is due almost in its entirety to an endemic circulation of the Q1 and Lisboa3 phylogenetic clades. Present data from 24-loci (not 12-loci) MIRU-VNTR allowed the subdivision of the Lisbon3 cluster in two other clusters herein designated as Lisboa3-A and –B. This data suggests two independent outbreaks, over the years, dated back to the 90s when the discrimination of Lisboa strains was identified by distinct rpoB mutations [8]. The Q1 spoligotyping data has revealed that this cluster is in fact intimately related with the B cluster identified in the 90s outbreak (unpublished data). Phylogenetic analysis based on previously published sets of SNPs [14, 37] revealed that Lisboa3 and Q1 strains formed distinct monophyletic evolutionary clades within the SCG 5 and PGG 2. Interestingly, M. tuberculosis F11 and the XDR-TB associated KZN strains, both originating from South Africa, also belong to SCG 5. Nevertheless a clear distinction is highlighted in the proposed phylogeny. This distinctiveness is also reflected by the RD comparison, but Lisboa3, Q1 and KZN strains appear to have an incongruent phylogeographic association using the RD typing. All these strains belong to the Euro-American lineage according to the RD classification proposed by Gagneux et al. [23]. However, the KZN strains included in the analysis showed to be positive for RD115, associated with an Americas/Europe sublineage, despite the fact that these strains are a major public health concern in South Africa, namely, the XDR-TB outbreak in KwaZulu Natal [38, 39]. The Lisboa3 and Q1 strains were on the other hand positive for RD174, associated with a West-African sublineage, but constitute a major public health concern in Europe. Present knowledge recognizes that RD174 is also associated with RD^RIO, an LSP that has initially been discovered in Rio de Janeiro, Brazil but was later found to be widespread. Historic ties connect Portugal, Brazil and West African Countries and a possible ancestor for these two clades might lie in Africa, more specifically on Portuguese Speaking African Countries. These phylogeographic incongruences are consistent with human migratory events out from, and back into, the African continent [12]. Moreover, these results also highlight that more is still needed to fully grasp the genetic diversity present within the SCG5 and LAM family as it encloses a high genetic diversity allied with a broad geographical distribution [40].

Another question still seems pertinent as to which selective advantages do these two clades possess allowing such high prevalence in this setting especially since other strains, e.g. pre-XDR-TB Beijing strains, also do circulate but at an apparent lesser prevalence? TB caused by RD^RIO strains has shown to be associated with weight loss, hemoptysis, higher bacillary loads and progression to cavitary disease [21, 41]. This deletion encompasses several PPE genes that have shown to be a potential source of immune variation (reviewed in [42, 43]) and hence, may constitute a pathogenic adaptation strategy to immune evasion. Higher bacillary loads are associated with a higher secondary case rate [44–46] and if in fact the absence of these genes truly plays an important role towards an increased virulence, or even transmissibility, it may be a factor that has contributed to the high prevalence of RD^RIO strains in this setting simultaneously contributing to the emergence and spread of M/XDR-TB strains.

Besides previously described RDs, the additional structural variants that were identified and that may be clade-specific could carry functional consequences that reflect host adaptation and selection. The finding that a 112 bp deletion is present among Lisboa3 clade strains, with a more restricted distribution than RD^RIO, affecting gene PPE41 might also provide additional clues and contribute to a higher virulence or transmissibility. PPE41 has been previously described has having an immunodominant nature and shown to activate a CD4⁺ and CD8⁺ mediated T cell response leading to an enhanced IFN-γ response as well as induce a strong humoral response [47, 48]. The deletion found might constitute a means of immune evasion and constitute a selective advantage over other circulating strains. More specifically, a stronger humoral response to PPE41 was found among extra-pulmonary TB patients [48]. The selective advantage provided by this deletion might therefore also be related with the fact that Lisboa strains were first identified among HIV-infected patients, which is associated with an increase in extra-pulmonary TB.

Phylogenetic context and microevolutionary trajectory of Lisboa3 and Q1 clades

The use of SNPs as molecular markers has contributed to an improved understanding of the evolutionary history of the M. tuberculosis complex. In the present study, given the availability of genomewide SNP data, a SNP-based phylogeny was deduced from the genomic data and, overall, the proposed phylogeny appears to be consistent with other SNP-based phylogenies although as already pointed out: SCG 3 does not exist as a monophyletic lineage but instead as a paraphyletic one. The original report by Filliol et al.[20] proposed a minimum number of sixteen SNPs that allowed assignment of any strain to an SCG but not to its subgroupings. A later erratum showed that SCG 3a belonged in fact to PGG1 while SCG 3b and 3c belonged to PGG2 as confirmed by our results. Alland et al. [37] proposed instead a set of nine SNPs that allowed strain assignment to any SCG and each subgroup [37].

The phylogeny constructed in the present study contributes nevertheless to demonstrate the uniqueness of Lisboa3 and Q1 strains in a global context and will comprise a future framework for genome-wide association studies (GWAS).

The phylogeny proposed also enabled a microevolutionary perspective on the path towards MDR and XDR. As expected, in the Lisboa3 and Q1 clades, INH resistance was found to be mediated by double inhA promoter/structural mutations, recently described by some of us to contribute to INH high-level resistance [27]. The acquisition of inhA C-15 T mutation was found to have occurred independently in both lineages, and in Q1 cluster it was possible to determine that C-15 T mutation was acquired at a first stage of INH high-level resistance development. In Lisboa3 it was not possible to determine which mutation appeared in the first place since no Lisboa3 isolate with a single inhA mutation was found. Recent work by Fenner et al. [49] suggested that inhA promoter mutations, more specifically C-15 T mutation, might be associated with Lineage 1 (Indo-Oceanic/SCG 1) [11, 49]. Nevertheless, another earlier study from Brimacombe et al. [50] showed that SCG 1 and 5 had all the mutations of interest towards INH resistance [50]. In our view, the fact that INH resistance in both Lisboa3 and Q1 clades is associated with inhA mutations, instead of the of the more usual KatG mutations, is possibly related with selective pressures exerted by the drug regimen itself.

The analysis of Lisboa3 subtree has further highlighted the M/XDR evolutive process in this clade. We have recently proposed an evolutionary path regarding drug resistance acquisition dynamics based on the acquisition of an eis promoter mutation as the first-step from MDR to XDR [6]. However the SNP phylogeny proposed is consistent with a twice and independent acquisition of an eis promoter mutation. Given this phylogeny it is not possible to establish any order of mutation acquisition. Nonetheless, instead of a single event, our analysis supports multiple development of XDR-TB in the same phylogenetic clade. Two different transmission chains involving strains with the RpoB S450L, instead of one, are also more likely since it is proposed that this mutation has also been acquired twice and independently [8].

Also important, he Lisboa3 XDR lineage characterized by gyrA D94G and eis G-10A mutations (node E1) will most likely present resistance to KAN, but not to CAP and AMK. If drug susceptibility testing to KAN is not included in the standard second-line drug panel of tested drugs, the strains belonging to this lineage will have an undetected XDR phenotype. An exception to this is the strain FF359_98 that bears a rrs A1401G mutation that leads to high-level KAN, AMK and CAP resistance [51].

One striking phylogenetic incongruence was found in two Q1 strains that lacked both second-line injectable drug and FQ genetic resistance determinants and at the same time sharing a recent common ancestor resistant to these two classes of drugs. These two strains were genotypically and phenotypically susceptible to amikacin (AMK), capreomycin (CAP) and any of the FQs tested. Two explanations may be considered: a phylogenetic misplacement, although the branches had a good statistical support or, these strains may descend from a reverter ancestor. Although theoretically possible, events such as these may be extremely rare. Only one report has documented an in-patient reversion of an isogenic strain from INH resistant to susceptible [52].

Compensatory evolution and RIF-resistance

The acquisition of further mutations in rpoA, rpoB or rpoC genes following RIF resistance development was recently demonstrated, using Salmonella as a model organism, to have an important role in fitness compensation, leading to a reduction in the doubling-time to values closer to the wild-type [53].It as also been demonstrated that rpoC gene has been target of convergent evolution [54]. In our microevolutionary analysis we have detected a RpoC mutation (K1152Q) occurring in the same branch as a RpoB S450L (equivalent to S531L in RpoB E. coli numbering). It is the first description of a putative compensatory mutation within the Lisboa3 clade, contributing to the success of one of its sub-lineages through the amelioration of the resistance fitness cost [28]. RIF compensatory evolution has been the subject of two recent studies that showed a high prevalence of rpoA and rpoC mutations mapped to the RpoA-RpoC interaction region [29, 30, 55]. The rpoC mutation described in a Lisboa3 sub-lineage does not fall in this region, nor was it described in these studies [29, 30]. Nevertheless, two other putative compensatory mutations mapping to the RpoA-RpoC interaction region were found in other isolates not belonging to the Lisboa3 or Q1 clades (Additional file 14). The putative role of these two latter mutations is only substantiated by the bioinformatic analysis of residue conservation. However, the putative compensatory role of the Lisboa3 K1152Q RpoC is further substantiated by their co-occurrence in the same branch as the RIF resistance determining mutation in rpoB. Furthermore, none of these putative compensatory mutations was previously described and may constitute novel polymorphisms associated with molecular RIF resistance compensation [18, 29, 30].

RpoB mutational analysis also allowed the identification of five putative compensatory mutations, of which one (L731P) was found to be acquired in the Q1 clade following RIF resistance acquisition through another rpoB mutation. This latter mutation was found to be homoplasic as it was also detected in different SCGs, which also points towards the usefulness of this mutation to counteract fitness costs imposed by the acquisition of other RIF resistance associated mutations. Mutations outside the RIF resistance determining region on rpoB gene have been described previously on RIF-resistant isolates with no mutations on this region [56, 57]. The mutations herein described as putatively compensatory were only considered as such if a mutation in the RRDR was already present providing further support for the compensatory role of the former.

The role of compensatory mutations in other loci and associated with compensation to resistance to other drugs than RIF have been identified and studied, namely, mutations on ahpC for INH or on rrs for second-line drug aminoglycosides [58–60]. Nevertheless, no compensatory mutations were identified in these genes (data not shown).

Still, it is yet possible that other mechanisms underlying resistance or compensation might lie elsewhere in the genome as even the role that synonymous SNPs play in gene expression must be reckoned with. Such an example in M. tuberculosis comes from a recent and elegant study by Safi et al.[61] in which a synonymous SNP on Rv3792 was found to act as an hypermorphic mutation on a downstream gene (embC), leading to an increase in EMB resistance [61]. In the same study, another type of mutation was found to be a key player at the multistep process of EMB resistance development – a neomorphic mutation on gene Rv3806c that increased the turnover of the decaprenylphophoryl-β-D-arabinose pathway, which also led to an increase in EMB resistance [61]. Two other recent studies from Zhang et al.[62] and Farhat et al.[54] also point to other genes that may be at play and under positive selection concerning drug resistance in M. tuberculosis[54, 62]. It becomes clear that functional characterization of the significant portion of the M. tuberculosis genes of unknown function must catch up the pace of high-throughput sequencing if a broader understanding of the genomic adaptation process is to be obtained.

IS6110transposition role in gene integrity and regulation

Insertion site mapping revealed a high genomic stability of insertion sequences other than IS6110. In fact, we have verified that only deletion events were responsible for variability regarding presence/absence of an insertion sequence other than IS6110. On the other hand, as demonstrated in this study IS6110 is a highly polymorphic marker, probably due to its rapid transposition rate [63, 64].

The finding that 65.5% of the IS6110 insertion sites mapped were located intragenically is in line with previous reports [65]. Considering that ≈ 91% of M. tuberculosis genome is composed by coding regions [66], it highlights the deleterious effects of transposition into certain genes essential to viability or to the successful completion of the pathogen’s infectious cycle. PGG1 strains, including the Beijing strains, were found to bear a higher number of IS6110 copies than PGG2 strains. IS6110 copy number is presumed to be under negative selection [67], however, in certain circumstances, it is the insertion site per se that might provide a selective advantage and not the copy number.

Considering the data obtained in this study, IS6110 is unarguably the species’ most important mobile element when considering transposition impact on genomic integrity. IS6110 appears to have an important role in genomic re-shaping towards adaptation either through localized disruption of putative antigenic targets (e.g. PPE/PE genes) or through its mobile promoter activity located in the IS6110 3′ end, capable of inducing transcription or upregulating the expression of downstream genes under stressful conditions [68]. We have found a considerable number of insertion sites to be within PPE genes and a more reduced number of sites to be upstream of PPE genes. PPE genes appear to have been positively selected in pathogenic mycobacteria, have important immune and antigenic potential, and some can induce a shift towards a Th2-type response [42, 43]. Not only the IS6110-mediated disruption of PPE genes might constitute a mean of immune evasion but it is also conceivable that upregulation of specific PPE genes might affect the Th1/Th2 response balance.

Remarkably, Lisboa3 and Q1 did not show any IS6110-mediated disruption of a PPE ORF, nor did we found any IS6110 upstream of a PPE gene. This fact perhaps demonstrates a different mode of evolution and host adaptation that does not require PPE gene modulation through IS6110 transposition.

Nevertheless, the maximum distance between an IS6110 and a downstream gene so that this 3′ promoter can exert its influence on gene expression is unknown. The results reported by Safi et al. [68] demonstrate that an IS6110 in M. tuberculosis 210 located 297 bp upstream of Rv1468c was associated with a threefold increase in transcription upon macrophage infection. Our results show that 15% of the mapped sites are located upstream of an ORF in proper orientation and at a distance of less than 300 bp which, at the light of present knowledge, fulfills the necessary assumptions to exert a putative upregulatory effect on those ORFs. Also considering the diversity of ORFs interrupted by IS6110 copies, gene knock out studies and assessment of downstream gene expression are necessary if a functional role for specific transposition events is to be established.

Spoligotyping lineage association with specific IS6110 sites has already been demonstrated, highlighting the phylogenetic informativeness of this marker [69]. Our results also support an association of specific IS6110 sites with strain lineage at both global and local levels.

Genome-wide SNP dynamics

Notably, the comparison of the distribution of SNPs by COG showed that N_s/S ratios vary through COG in a lineage-dependent manner. Although, Lisboa3 and Q1 isolates might be overrepresented in the analysis due to the high prevalence in the community, we have shown that Lisboa3 and Q1 present statistically different N_s/S ratios from the remaining isolates. We propose that differences in N_s/S COG might highlight different evolution strategies selected during host-pathogen interaction and adaptation.

Moreover, an overall higher N_s/S ratio was observed for the first quadrant and an overall lower N_s/S ratio was observed for the second quadrant revealing heterogeneous N_s/S ratios negatively correlated with the T_v/T_s ratio. The biochemical nature behind this T_v/T_s ratio heterogeneity requires further studies as it may be driving localized higher non-synonymous mutation rates with functional impact on strain evolution. The precise genes affected by non-synonymous mutations within these COG categories merit further studies as each COG includes a considerable number of genes, that mutated might enhance the transmissibility or drug resistance, and should be analyzed in a systems biology perspective using in silico models [70–72].

The finding that terminal branches of the subtrees analyzed had lower N_s/S ratios than the inner branches was contrary to the findings of Namouchi et al. [36]. Namouchi et al. suggested that non-synonymous changes might be purged by natural selection yielding higher N_s/S ratios. However, an opposite view is also possible: non-synonymous mutations are favored by natural selection yielding the same higher N_s/S ratio, especially as a mean of adaptation in an organism devoid, or with low, horizontal gene transfer such as M. tuberculosis[73, 74]. From our data, we can say that non-synonymous mutations may be favored in the inner branches of both subtrees as a mean to develop and adapt to drug resistance, yielding a higher N_s/S ratio that is consistent with a reduced purifying selection [12]. Nevertheless, the differences between the results from both studies might lie in the fact that the N_s/S ratio analysis herein presented was performed for two sets of strains that are in a much more closer time frame in order to understand microevolution within two clades.

WGS and molecular epidemiology

In the present study we have shown that in Lisbon, Portugal, where the MDR-TB situation had already escalated to a XDR-TB situation, it is mainly caused by transmission of two unique phylogenetic clades. We have previously shown that XDR-TB was already a reality in Portugal during the 1990s [6] but noteworthy, the data from the present study clearly shows that these strains belong to the same phylogenetic Lisboa3 M/XDR sublineages that are still presently in circulation. The uniqueness of these strains was revealed by a distinct phylogenetic placement within SCG 5.

The Lisboa3 clade belongs to a much broader group of strains that usually share at least 95% of MIRU-VNTR pattern similarity: the Lisboa family [7]. One of the future goals is to better understand the populational structure of this family of strains, from which the Lisboa3 clade has differentiated, as strains belonging to this family have previously shown the potential to evolve to MDR-TB [4].

Another important point coming from the present student is the utility of WGS for epidemiological surveillance and strain typing. WGS presents an advantage over the classical typing methods (RFLP-IS6110, Spoligotyping or MIRU-VNTR) as it enables picturing transmission at a much higher resolution and ascertain isolate relatedness using well described models of molecular evolution. In the present study, WGS allowed strain discrimination within MIRU-VNTR clusters and distinguish between three independent Lisboa3 MDR sub-clades. Despite the technical difficulties in data analysis, as WGS costs converge towards the cost of MIRU-VNTR, the former is likely to replace MIRU-VNTR as the gold-standard for molecular epidemiological surveillance and strain typing. WGS can also enable more focused contact tracing by reducing the number of plausible genomically linked cases to investigate, leading to an improved case detection. WGS-assisted routine surveillance is still far away for many settings, but as this technology becomes gradually available to the mycobacteriology laboratory it will also be expected a greater understanding of TB transmission. In a recent study, Walker et al. have defined a threshold of 12 SNPs of difference, above which recent transmission can be excluded [75]. In our study, the number of unique SNPs to each isolate determined for the Lisboa3 and Q1 clades are consistent with ongoing recent transmission. This finding allied with the genomic uniqueness of these strains are of special importance not only locally but in a macro-epidemiological context. It is likely that these strains may spread to other parts of the world, due to increasing global travel and migratory waves, and be the cause of additional public health concern [76–78]. A recent report of an RD^RIO strain recovered from a remote location in Tibet alerts to this possibility [79].

It is also worth having in consideration that the host residing bacilli population has a certain degree of heterogeneity that can be overlooked through WGS but nonetheless lead to a higher than expected genomic difference after transmission. Similarly, selection during drug treatment might also artificially extend genomic distances. In this regard, classical genotyping using more stable markers might prove useful. The present study also stresses the need of further genomic studies in order to contribute to a M. tuberculosis genome-wide evolutive scenario, representative of different settings.

This, together with clinical data, will ultimately enable GWAS with a positive impact in TB management.

Isolates and genetic data

The study consists of 56 M. tuberculosis clinical isolates (source: 55 Lisbon, 1 Oporto) recovered from laboratories and hospital units across Lisbon Health region. This set of isolates comprises a convenience sample of M. tuberculosis clinical isolates received for genotypic analysis at the Mycobacteria Laboratory from the Faculty of Pharmacy of the University of Lisbon. The sample is composed by drug resistant isolates plus additional susceptible isolates found to be genetically close (MIRU-VNTR) to the drug resistant isolates. All isolates underwent drug susceptibility testing for INH, RIF, STP, EMB and PZA and second-line drugs using standard methods (see [4]). DNA extraction was performed from culture growth on Lowenstein-Jensen medium slants using the cetyl trimethylammonium bromide methodology [80]. The DNA was used in genotyping by the 24-loci MIRU-VNTR method (see previous work, [81]). Extracted DNA was also subjected to whole-genome (101 bp paired end) sequencing at the KAUST genomics facility using the Illumina HiSeq 2000 platform (500 bp insert size). We also complemented this data using sequences in the public domain (F11, CDC1551, KZN1435, KZN4207, KZN605, KZN_R506, KZN_V2475, UT205 RGTB327, RGTB423, CCDC5180, CCDC5079, CTRI-2, BTB05_552, BTB05_559, S96_129, HN878, R1207, and X122 (all from the NCBI database).

Genomic variant detection

The raw Illumina sequencing data was aligned to the H37Rv reference genome using the Burrows-Wheeler Alignment Tool v.0.6.1, yielding high coverage data for all isolates (mean read depth per position, mean 249.9, range 44–1411 fold; mean 99.1% genome covered, range 98.6 - 99.9%) (Table 1) [82]. Single nucleotide polymorphisms (SNPs) and small indels (<30 bp) were called using SAMtools software (v0.1.18) [83]. Other small indels (<100 bp) were detected using the software Pindel [84]. Only variants supported by at least ten sequence reads were considered. Detection of larger structural variants was performed using the SVMerge v1.2 pipeline combining Pindel v0.2.4 t, BreakDancer v1.1 Cpp package and, SECluster analysis outputs [85, 86]. Structural variant detection was done for each isolate alone and validation was achieved using comparison with local de novo assembly using Velvet [87]. Loci reported to be associated with regional differences [22, 23, 88] were identified using the alignments and coverage.

For insertion sequence (IS) mapping, reads containing specific oligonucleotide sequence of both 5′ and 3′ extremities (listed in Additional file 22) were extracted and, flanking genomic regions of interest concatenated in FASTA format producing one file for each extremity for each strain. Local BLAST analysis (standalone NCBI BLAST v.2.2.27+) was carried out for each file against M. tuberculosis H37Rv reference genome, minimum supporting read depth used to as a quality filter (10 for isolates with >500 fold coverage, 2 for the remaining). For IS6110 BLAST hits, a mapping quality classification scheme was established consisting in high confidence, medium confidence and lesser confidence sites. Paired sites corresponding to mapping of both 5′ and 3′ ends in all isolates on which it occurred were classified as high confidence sites. Paired insertion sites for which both ends were mapped in at least 50% of the isolates on which they were found to occur were considered of medium confidence. Insertion sites in which only one end of the IS6110 was mapped were considered of lesser confidence. Furthermore, insertion site hits mapped to M. tuberculosis H37Rv were excluded to avoid repetitive mapping.

Other bioinformatics

The genomic data of publicly available M. tuberculosis strains (format FASTA) were included in the analysis through conversion to FASTQ format reads using the program dwgsim v.0.1.10, and mapped and analyzed as described above. When necessary, DNA sequence alignment was performed using the CLC Sequence Viewer v7.6.1 (CLC bio®, Aarhus N, Denmark) and visualized in BioEdit v7.1.3.0 (T. Hall).

A MIRU-VNTR-based dendrogram was constructed in the public MIRU-VNTRplus database using the D_sw measure of genetic distance for tandem repeat loci [89] and the Unweighted Pair Group Method with Arithmetic Averages (UPGMA) clustering method. Spoligotyping profile was inferred from raw read data using the SpolPred software followed by comparison to the SITVIT WEB database [40, 90]. A phylogenetic tree based on SNPs was constructed using Seaview 4.3.5 [91] using the Maximum Likelihood method. The analysis involved 76 nucleotide sequences with a total of 11271 sites in the final dataset. Tree topology was tested using the most recent approximate Likelihood Ratio Test (aLRT) as an alternative to bootstrap.

Putative impact of selected compensatory mutations on protein function was assessed through the use of SIFT scores (available at http://sift.jcvi.org/) [31] computed from the query alignment against UniRef90 database hits (with less than 90% identity, with a median sequence conservation equal to 3.00).

Any statistical analysis was conducted using the SPSS software.

Data access

All sequencing data have been submitted to the European Nucleotide Archive (http://www.ebi.ac.uk/ena/) under study accession number ERP002611. Phylogenetic data (tree and alignment matrix) have been submitted to TreeBase under Study ID no. 16158 (URL: http://treebase.org/treebase-web/home.html).