Volume 25, Number 3—March 2019
SNP-IT Tool for Identifying Subspecies and Associated Lineages of Mycobacterium tuberculosis Complex
Mycobacterium tuberculosis complex (MTBC) encompasses a group of organisms that cause tuberculosis (TB) in humans and animals. TB in humans is caused mainly by M. tuberculosis but also by other members of MTBC, including the less well understood animal-associated subspecies M. bovis, M. caprae, M. pinnipedii, M. suricattae, M. orygis, M. microti, and M. mungi (1–6). The global burden of zoonotic TB is thought to be both underestimated and increasing (7); however, accurate assessment of prevalence is made difficult by a lack of clinical diagnostic tools and surveillance (8).
Efforts to differentiate members of MTBC and study the phylogeny of the complex have thus far included analysis of large genomic deletions (9), variable-number tandem-repeats (VNTR), spacer oligonucleotide typing (spoligotyping), multilocus sequence typing, and, more recently, single-nucleotide polymorphism (SNP)–based phylogenies (10). Numerous tools now exist that make in silico predictions of lineages within the complex from whole-genome sequencing (WGS) data using a variety of approaches, including the detection of single SNPs from both unassembled and mapped genomes, comparison of de Brujin graphs, and MinHash-based comparisons (11–14). None of these tools has yet been calibrated to reliably differentiate among all subspecies, particularly the animal-associated ones, whose incidence and clinical significance are likely to be underestimated as a result.
The host ranges of the various MTBC subspecies differ, which has serious implications for contact investigations and source case finding. For example, M. microti is found in wild cats and rodents and causes human infection, usually in association with rodent contact (15). In contrast with infections caused by M. bovis, most M. microti infections have been reported to cause pulmonary TB, which raises the possibility of onward transmission, although such transmission has not yet been reported (16). M. pinnipedii, which causes TB in seals, is sometimes transmitted to humans during outbreaks in zoos or wildlife parks (17). Although isolated mostly from gazelle species, M. orygis has also been seen in humans in recent years, although how humans contract this bacterium is still unclear (2). M. mungi causes disease in banded mongooses, the dassie bacillus causes disease in rock hyraxes, and M. suricattaecauses disease in meerkats, but none of these bacteria is currently known to cause disease in humans (3,5).
The spectrum of clinical phenotypes associated with human infection by MTBC animal lineages is largely unknown, partly because the identification of these organisms is currently difficult. Accurate identification of the causative subspecies in all cases would enable characterization of disease associated with animal lineages and of diversity in clinical phenotypes, which would contribute to better disease management. A higher level of knowledge on the spread and host range of the subspecies would also provide a better basis on which to study the history of the evolutionary development of the complex as a whole. Now that routine WGS is being performed by Public Health England (PHE) for all MTBC isolates, we sought to use these data to estimate the burden of animal-associated TB in England. Therefore, we identified a broad panel of SNPs that define each subspecies, lineage, and sublineage within the MTBC and assessed them using a new SNP-based tool, SNPs to Identify TB (SNP-IT).
We defined a set of isolates (N = 323) from which to identify SNPs associated with subspecies, lineages, and sublineages within the MTBC (Figure 1). We identified isolates from the collection of the National Institute for Public Health and the Environment (RIVM; Bilthoven, Netherlands) using a combination of spoligotyping patterns, SNPs, restriction fragment length polymorphism (RFLP) patterns, the hybridization patterns in the HAIN Genotype MTBC assay, polymorphic GC-rich repeat sequence (PGRS) profiles, and VNTR patterns in accordance with current and previous standard practice (Appendix Table 1). We identified isolates from a whole-genome sequencing archive held at the University of Oxford (Oxford, UK) using the SNP typing system of Stucki et al. (18) for Mtb lineages 1–4 and by clustering on a maximum likelihood tree with the isolates from the Netherlands for lineages 5 and 6. In addition, we used published strains (n = 5) to be able to include M. suricattae, M. mungi, and the dassie bacillus, which were not present in the Oxford or RIVM collections. The nomenclature we adopted for this study is summarized in Appendix Table 2.
We applied parallel bioinformatics approaches to assess applicability across pipelines. As such, we independently mapped reads from Illumina platforms (https://www.illumina.com) to 2 different versions of the H37Rv reference genome. We mapped reads to NC000962.3 with Breseq version 0.28.1 (http://barricklab.org/twiki/bin/view/Lab/ToolsBacterialGenomeResequencing), using a minimum allele frequency of 80% and minimum coverage of 5 times for SNP calls. Separately, we mapped reads again to NC000962.2, for which we used Snippy version 3.1 with default settings (minimum coverage 10 times, minimum allele frequency 90%) (19,20). We extracted all SNPs shared exclusively by isolates of each subspecies, lineage, and sublineage identified by both pipelines. The lineage-defining positions for lineage 4 are not variants with respect to the reference, itself lineage 4, but are uniquely conserved positions. We therefore identified these positions by mapping a core SNP alignment to a maximum-likelihood tree using Mesquite version 3.30 (21). These nucleotide loci were added to the catalog of phylogeny-determining SNPs.
All newly sequenced genomes are available from the National Center for Biotechnology Information under project accession no. PRJNA418900. The SNP-IT tool, including all relevant reference libraries used for this study, is available as an open-source package online (https://github.com/samlipworth/snpit).
To validate the algorithm, we compiled an independent collection of genomes (N = 511) using clinical isolates sequenced by PHE Birmingham, UK, identified as MTBC and not included in the calibration set (Figure 1). We augmented them with data from the European Nucleotide Archive and the National Center for Biotechnology Information Sequence Read Archive to increase the representation of animal subspecies (N = 103; Appendix Table 2). To maximize inclusion of animal isolates from the public archives, we used the new Colored Bloom Graph (CBG) software (22). Using CBG, we searched a snapshot of the Sequence Read Archive (to December 2016, N = 455,632) with our new set of reference kmers for Mykrobe predictor (see Comparison with Existing Tools).
We compared FASTA files of the whole genome (https://www.ebi.ac.uk/Tools/sss/fasta) to the catalog of phylogeny-determining SNPs to make predictions for PHE isolates, whereas for isolates downloaded from the nucleotide archive, we created only limited variant calling format files using Snippy (only SNPs with respect to the reference genome were included to increase computational efficiency). To ensure that genomic loci defining lineage 4 were included, we used a mutated reference genome to create these limited variant calling format files. We compared SNPs in the query sample with reference libraries of lineage-specific SNPs for each clade. We assigned query genomes to particular subspecies or lineages if >10% of lineage- or subspecies-specific SNPs were detected in the strain in question. We assessed all predictions against the maximum-likelihood phylogeny. For M. mungi, we could locate only 1 genome in the public sequence libraries, so we could not validate this subspecies.
To assess the caseload across the different members of the MTBC seen by the PHE laboratory in Birmingham, we applied the algorithm to 3,128 MTBC genome sequences from consecutively obtained clinical isolates. H37Rv is routinely sequenced by the laboratory on WGS plates; these isolates were not removed, and their identification served as an internal control.
Comparison with Existing Tools
We first compared strain characterization by our new SNP-IT tool with those of KvarQ (https://github.com/kvarq/kvarq), TB-Profiler (http://tbdr.lshtm.ac.uk), and Mykrobe predictor (http://www.mykrobe.com/products/predictor) on default settings and then after integrating our updated SNP library. To enable our new data to be integrated with published SNP libraries (23) and for practical reasons when modifying existing tools, we created a minimal SNP dataset. We filtered our larger SNP catalog for synonymous SNPs that occurred in coding regions (as annotated by SnpEff version 4.3 ) and selected 1 representative SNP for each subspecies, lineage, and sublineage at random. We then modified the existing software packages to include reference SNPs (or kmers for Mykrobe predictor) for the subspecies, lineages, or sublineages that they initially failed to identify.
Calibration and Validation
In total, we identified 13,893 SNPs (median of 229 SNPs per group, interquartile range 296) as predictive of taxonomic and phylogenetic groups of interest (Appendix Table 3). The greatest number of phylogenetic SNPs was seen in M. canettii (n = 6,837) and the fewest in M. bovis (n = 23). Subspecies that arise from common deep branches, such as M. microti and M. pinnipedii (Figure 2), have lower numbers of unique phylogenetic SNPs (n = 128 for M. microti and n = 301 for M. pinnipedii) than those that do not, such as M. orygis (n = 781). All predictions made by SNP-IT across all the subspecies, lineages, and sublineages were consistent with the maximum-likelihood phylogeny for all isolates in the validation set (Table 1).
Determining Prevalence of Animal Subspecies in a Collection of Clinical Isolates
We retrospectively applied SNP-IT to clinical isolates sequenced as part of the routine PHE diagnostic workflow in Birmingham to estimate the prevalence of the animal subspecies among MTBC samples. Of 3,128 samples from 2,106 patients for which there was a whole-genome sequence available, we identified 24 as M. orygis, 3 as M. microti, 34 as M. bovis, and 1 as M. caprae (Table 2). In the case of M. orygis, we further investigated whether there was any genomic signal of possible person-to-person transmission. We identified 2 such instances, 1 in which the pairwise genetic distance between 2 patients was 0 SNPs and a second in which it was 6 SNPs (Appendix Figure 1).
Phylogenetic SNPs in Drug Resistance–Associated Genes
Using a previously published list of drug resistance–associated genes for M. tuberculosis (25), we searched all subspecies for phylogenetic SNPs in drug resistance–associated genes (AppendixTable 4). All subspecies contain unique phylogenetic SNPs (N = 95 in total) in these genomic regions, but on the basis of our data, we were unable to determine whether any of these mutations are linked to lineage-specific resistance because we did not have the corresponding phenotypic drug susceptibility testing data.
Comparison with Existing Software
Compared with SNP-IT for the clinical set of isolates, Mykrobe predictor reported M. orygis as M. tuberculosis West African lineage and M. pinnipedii as M. microti. KvarQ identified all animal-associated subspecies only as “animal lineage.” TB-Profiler was unable to delineate among animal subspecies, which were all reported as M. bovis/M. tuberculosis West African lineage.
After we modified the KvarQ, TB-Profiler, and Mykrobe predictor databases with our minimal SNP catalog, all systems agreed on the identity of all the MTBC isolates in the clinical set. SNP-IT was unable to identify 10 samples because <10% of type-specific SNPs were present in these strains. This result was because our pipeline made no call at the lineage informative sites because of the presence of a minor allele, most likely the result of contamination or a mixture of 2 strains in the sample. However, TB-Profiler and Mykrobe predictor were both able to identify 2 of these isolates as mixed Beijing lineage/M. orygis using our new minimal SNP dataset.
SNP typing is a powerful method for discriminating among members of MTBC, which are often not discernible by conventional laboratory methods. The SNP databases of Stucki, Coll, and Comas are currently used as the knowledge base for KvarQ, TB-Profiler, and Mykrobe predictor (18,23,26). None of these databases, however, provides adequate resolution for the animal subspecies. In contrast, SNP-IT was able to assign subspecies, lineages, and sublineages to all samples in the validation set with 100% accuracy compared with maximum-likelihood phylogeny. Implementing this fine-resolution algorithm into a routine diagnostic workflow would be a major step toward understanding the epidemiology and pathogenicity of the less common members of MTBC. All 3 existing systems tested (KvarQ, Mykrobe predictor, and TB-Profiler) were identical in performance when given the same SNP reference database, demonstrating that the clinically meaningful differences highlighted in a recent review are easily ameliorated (27).
By applying SNP-IT to a clinical dataset, we discovered an unexpectedly high number of animal subspecies among MTBC isolates, particularly M. orygis, from humans in the United Kingdom. This recently described member of the complex has a host range that includes waterbucks, gazelles, rhesus monkeys, cows, and rhinoceri (2,28,29). Several human cases have been described in patients in the Netherlands of South/Southeast Asian origin (2), but no cases have been described in the United Kingdom. Human-to-animal transmission has been described in 1 case in New Zealand (30). Given that zoonotic TB is associated with higher rates of extrapulmonary disease and may be less likely to grow in culture (31,32), retrospective interrogation of WGS libraries, as in this study, is likely to underestimate the true burden of disease.
Given the large amount of resources aimed at controlling bovine TB, it is noteworthy that another zoonosis is seen at a similar rate in this collection of clinical isolates. This finding raises questions about the host range and transmission of M. orygis, with potential implications for TB control both in animals and humans. To recognize the particulars of the clinical phenotype, epidemiology, and optimal management strategy of M. orygis infection, it is first crucial to accurately distinguish these cases from M. tuberculosis West African lineages (5,6). This discernment is currently not possible by either the Hain Genotype MTBC molecular probe or existing SNP-based platforms. We identified 2 pairs of nearly identical M. orygis isolates that could be compatible with either person-to-person transmission or, possibly, common exposure to the same infected animal. From an international perspective, the role of M. orygis in zoonotic transmission in Africa, Asia, and other high-prevalence settings with extensive animal contact is poorly understood and may warrant further investigation.
All the animal subspecies had phylogenetic SNPs in drug resistance–associated genes. When these genes are not known to be associated with drug resistance, they can be helpfully annotated as such by diagnostic algorithms and excluded for the purpose of predicting susceptibility. An unavoidable weakness of any SNP-based approach is its vulnerability to null-calls as a result of minor alleles at informative positions or a lack of coverage. SNP-IT uses the entire library of subspecies/lineage/sublineage-defining SNPs such that this weakness is not an issue unless it occurs at most of these positions. An additional limitation is that, although we have sought to calibrate SNP-IT using the most diverse collection of samples available to us, it may not be able to correctly identify isolates that originate from deeper phylogenetic branches than those in our calibration set.
In conclusion, in this study we demonstrate a higher-than-expected burden of zoonotic TB in a large collection of clinical isolates from the United Kingdom. The SNP-IT tool we have developed will help researchers to examine the epidemiology of zoonotic TB in a global context, as well as optimizing the disease’s clinical management. As more healthcare systems begin to routinely use WGS, there is an opportunity to accurately diagnose the causative subspecies of TB in all cases, which will enable identification of previously underrecognized zoonoses and reverse zoonoses and implementation of control interventions in the interests of One Health.
Dr. Lipworth is a clinical fellow with the Modernising Medical Microbiology research group in the Nuffield Department of Medicine, University of Oxford. His research interests lie in the clinical applications of whole-genome sequencing, particularly diagnostics and molecular epidemiology. Ms. Jajou is a PhD candidate supervised by Dick van Soolingen at the National Institute for Public Health and the Environment in the Netherlands. Her primary research interest is the use of whole-genome sequencing to improve diagnostics and public health practice in tuberculosis.
We thank Phillip Fowler his work improving the SNP-IT Python package and Pretin Davda for his assistance in linking isolates to patient level.
This research was supported by the National Institute for Health Research (NIHR) Health Protection Research Unit in Healthcare Associated Infections and Antimicrobial Resistance at the University of Oxford in partnership with Public Health England (PHE) and by Oxford NIHR Biomedical Research Centre. T.P. is a NIHR senior investigator. The report presents independent research funded by NIHR. The views expressed in this publication are those of the authors and not necessarily those of the National Health Service, NIHR, the Department of Health, or PHE. Z.I. is a Sir Henry Dale Fellow jointly funded by the Wellcome Trust and the Royal Society (grant no. 102541/Z/13/Z).
- Magee J, Ward A. Mycobacterium. In: Whitman WB, editor. Bergey’s Manual of Systematics of Archaea and Bacteria. Hoboken (NJ): John Wiley & Sons; 2015. p. 2–4.
- van Ingen J, Rahim Z, Mulder A, Boeree MJ, Simeone R, Brosch R, et al. Characterization of Mycobacterium orygis as M. tuberculosis complex subspecies. Emerg Infect Dis. 2012;18:653–5.
- Parsons SDC, Drewe JA, Gey van Pittius NC, Warren RM, van Helden PD. Novel cause of tuberculosis in meerkats, South Africa. Emerg Infect Dis. 2013;19:2004–7.
- Mostowy S, Cousins D, Behr MA. Genomic interrogation of the dassie bacillus reveals it as a unique RD1 mutant within the Mycobacterium tuberculosis complex. J Bacteriol. 2004;186:104–9.
- Alexander KA, Laver PN, Michel AL, Williams M, van Helden PD, Warren RM, et al. Novel Mycobacterium tuberculosis complex pathogen, M. mungi. Emerg Infect Dis. 2010;16:1296–9.
- Fabre M, Hauck Y, Soler C, Koeck J-L, van Ingen J, van Soolingen D, et al. Molecular characteristics of “Mycobacterium canettii” the smooth Mycobacterium tuberculosis bacilli. Infect Genet Evol. 2010;10:1165–73.
- Olea-Popelka F, Muwonge A, Perera A, Dean AS, Mumford E, Erlacher-Vindel E, et al. Zoonotic tuberculosis in human beings caused by Mycobacterium bovis-a call for action. Lancet Infect Dis. 2017;17:e21–5.
- Müller B, Dürr S, Alonso S, Hattendorf J, Laisse CJM, Parsons SDC, et al. Zoonotic Mycobacterium bovis-induced tuberculosis in humans. Emerg Infect Dis. 2013;19:899–908.
- Brosch R, Gordon SV, Marmiesse M, Brodin P, Buchrieser C, Eiglmeier K, et al. A new evolutionary scenario for the Mycobacterium tuberculosis complex. Proc Natl Acad Sci U S A. 2002;99:3684–9.
- Jagielski T, van Ingen J, Rastogi N, Dziadek J, Mazur PK, Bielecki J. Current methods in the molecular typing of Mycobacterium tuberculosis and other mycobacteria. Biomed Res Int. 2014;2014:645802.
- Steiner A, Stucki D, Coscolla M, Borrell S, Gagneux S. KvarQ: targeted and direct variant calling from fastq reads of bacterial genomes. BMC Genomics. 2014;15:881.
- Coll F, McNerney R, Preston MD, Guerra-Assunção JA, Warry A, Hill-Cawthorne G, et al. Rapid determination of anti-tuberculosis drug resistance from whole-genome sequences. Genome Med. 2015;7:51.
- Ondov BD, Treangen TJ, Melsted P, Mallonee AB, Bergman NH, Koren S, et al. Mash: fast genome and metagenome distance estimation using MinHash. Genome Biol. 2016;17:132.
- Bradley P, Gordon NC, Walker TM, Dunn L, Heys S, Huang B, et al. Rapid antibiotic-resistance predictions from genome sequence data for Staphylococcus aureus and Mycobacterium tuberculosis. Nat Commun. 2015;6:10063.
- van Soolingen D, van der Zanden AGM, de Haas PEW, Noordhoek GT, Kiers A, Foudraine NA, et al. Diagnosis of Mycobacterium microti infections among humans by using novel genetic markers. J Clin Microbiol. 1998;36:1840–5.
- Xavier Emmanuel F, Seagar A-L, Doig C, Rayner A, Claxton P, Laurenson I. Human and animal infections with Mycobacterium microti, Scotland. Emerg Infect Dis. 2007;13:1924–7.
- Kiers A, Klarenbeek A, Mendelts B, Van Soolingen D, Koëter G. Transmission of Mycobacterium pinnipedii to humans in a zoo with marine mammals. Int J Tuberc Lung Dis. 2008;12:1469–73.
- Stucki D, Malla B, Hostettler S, Huna T, Feldmann J, Yeboah-Manu D, et al. Two new rapid SNP-typing methods for classifying Mycobacterium tuberculosis complex into the main phylogenetic lineages. PLoS One. 2012;7:e41253.
- Seemann T. Snippy: Rapid haploid variant calling and core SNP phylogeny; 2015 [cited 2017 May 1]. https://github.com/tseemann/snippy
- Deatherage DE, Barrick JE. Identification of mutations in laboratory-evolved microbes from next-generation sequencing data using breseq. Methods Mol Biol. 2014;1151:165–88.
- Maddison WP, Maddison DR. Mesquite: a modular system for evolutionary analysis; 2017 [cited 2017 May 1]. http://mesquiteproject.org
- Bradley P, den Bakker H, Rocha E, McVean G, Iqbal Z. Real-time search of all bacterial and viral genomic data. bioRxiv; 2017 [cited 2017 Dec 18]. https://www.biorxiv.org/content/early/2017/12/18/234955.abstract
- Coll F, McNerney R, Guerra-Assunção JA, Glynn JR, Perdigão J, Viveiros M, et al. A robust SNP barcode for typing Mycobacterium tuberculosis complex strains. Nat Commun. 2014;5:4812.
- Cingolani P, Platts A, Wang L, Coon M, Nguyen T, Wang L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin). 2012;6:80–92.
- Walker TM, Kohl TA, Omar SV, Hedge J, Del Ojo Elias C, Bradley P, et al.; Modernizing Medical Microbiology (MMM) Informatics Group. Whole-genome sequencing for prediction of Mycobacterium tuberculosis drug susceptibility and resistance: a retrospective cohort study. Lancet Infect Dis. 2015;15:1193–202.
- Comas I, Homolka S, Niemann S, Gagneux S. Genotyping of genetically monomorphic bacteria: DNA sequencing in Mycobacterium tuberculosis highlights the limitations of current methodologies. PLoS One. 2009;4:e7815.
- Schleusener V, Köser CU, Beckert P, Niemann S, Feuerriegel S. Mycobacterium tuberculosis resistance prediction and lineage classification from genome sequencing: comparison of automated analysis tools. Sci Rep. 2017;7:46327.
- Thapa J, Paudel S, Sadaula A, Shah Y, Maharjan B, Kaufman GE, et al. Mycobacterium orygis–associated tuberculosis in free-ranging rhinoceros, Nepal, 2015. Emerg Infect Dis. 2016;22:570–2.
- Rahim Z, Thapa J, Fukushima Y, van der Zanden AGM, Gordon SV, Suzuki Y, et al. Tuberculosis caused by Mycobacterium orygis in dairy cattle and captured monkeys in Bangladesh: a new scenario of tuberculosis in south Asia. Transbound Emerg Dis. 2017;64:1965–9.
- Dawson KL, Bell A, Kawakami RP, Coley K, Yates G, Collins DM. Transmission of Mycobacterium orygis (M. tuberculosis complex species) from a tuberculosis patient to a dairy cow in New Zealand. J Clin Microbiol. 2012;50:3136–8.
- Sanoussi CN, Affolabi D, Rigouts L, Anagonou S, de Jong B. Genotypic characterization directly applied to sputum improves the detection of Mycobacterium africanum West African 1, under-represented in positive cultures. PLoS Negl Trop Dis. 2017;11:e0005900.
- Dürr S, Müller B, Alonso S, Hattendorf J, Laisse CJM, van Helden PD, et al. Differences in primary sites of infection between zoonotic and human tuberculosis: results from a worldwide systematic review. PLoS Negl Trop Dis. 2013;7:e2399.
TablesCite This Article
Original Publication Date: 1/31/2019
1These authors contributed equally to this article.
2Joint senior authors.