Inferring patient to patient transmission of Mycobacterium tuberculosis from whole genome sequencing data

Abstract

Background: Mycobacterium tuberculosis is characterised by limited genomic diversity, which makes the application of whole genome sequencing particularly attractive for clinical and epidemiological investigation. However, in order to confidently infer transmission events, an accurate knowledge of the rate of change in the genome over relevant timescales is required.

Methods: We attempted to estimate a molecular clock by sequencing 199 isolates from epidemiologically linked tuberculosis cases, collected in the Netherlands spanning almost 16 years.

Results: Multiple analyses support an average mutation rate of ~0.3 SNPs per genome per year. However, all analyses revealed a very high degree of variation around this mean, making the confirmation of links proposed by epidemiology, and inference of novel links, difficult. Despite this, in some cases, the phylogenetic context of other strains provided evidence supporting the confident exclusion of previously inferred epidemiological links.

Conclusions: This in-depth analysis of the molecular clock revealed that it is slow and variable over short time scales, which limits its usefulness in transmission studies. However, the superior resolution of whole genome sequencing can provide the phylogenetic context to allow the confident exclusion of possible transmission events previously inferred via traditional DNA fingerprinting techniques and epidemiological cluster investigation. Despite the slow generation of variation even at the whole genome level we conclude that the investigation of tuberculosis transmission will benefit greatly from routine whole genome sequencing.

Figure 1

Maximum likelihood phylogeny of 199 M.tuberculosis strains from 11,879 single nucleotide polymorphisms. Lineages are indicated as described by [22].

Figure 2

Pairwise genetic distances between M. tuberculosis isolates from patients linked by contact tracing. The genetic distances consist of the number of identified single nucleotide polymorphisms that differed between the genomes of two linked isolates.

Figure 3

Poor correlation between time and number of SNPs accumulated in the secondary case isolate for drug resistant and sensitive isolates. Three pairs were excluded (see main text). SNPs conferring drug resistance were also removed. Resistant isolates are classed as isolates phenotypically resistant to atleast isoniazid, streptomycin, ethambutol or rifampicin.

Figure 4

Per lineage root to tip plot. Lineages were rooted using their topology in the entire maximum likelihood tree, and the number of SNPs accumulated from the root was plotted against date of isolation. Correlation is poor for all lineages, with r squared values of 0.002,0.03,0.006 and 0.06 for Euro American, East African.

Figure 5

Date of collection vs. root to tip SNP distance of the 5 largest clusters. A. Clusters were rooted using their topology in the entire maximum likelihood tree. Linear regression was fitted using Path-O-Gen [21]. B. Data combined from A. Time represents days since first isolation in the cluster. Shaded area indicates 95% confidence interval.

Figure 6

Exclusion of epidemiologically linked pairs based on phylogenetic position. Red indicates primary case isolate, blue is the secondary case isolate. A: excluded pair 1. B: excluded pair 2. C: Excluded pair with SNP difference of 149. N09501026 shares a deletion with an isolate in a different study [36]. D: Example of expected phylogenetic positioning of direct transmission pairs, brackets indicate paired isolates.

Figure 7

Pairwise SNP differences between isolates. A: Pairwise SNP difference between all 199 isolates, the peaks represent the pairwise differences between the major lineages. Box indicates subsection shown in B: pairwise SNP differences under 150 SNPs for linked and unlinked pairs. There are many unlinked pairs that have SNP distances which overlap with the distribution of SNP distances for linked pairs.