Large-scale whole genome sequencing of M. tuberculosis provides insights into transmission in a high prevalence area

1Despite the huge global burden of tuberculosis, the factors influencing transmission remain poorly understood. Compared to other bacteria, the genome of Mycobacterium tuberculosis is stable and genetic variation was thought to be limited, but with increased sequencing, greater diversity has been recognized. Based on the genotype, M. tuberculosis has seven lineages: three ‘ancient’ (lineage-1 and two Mycobacterium africanum lineages), and three ‘modern’ (lineages-2, 3, 4), and one intermediate (lineage-7), recently described in Ethiopia. The lineages may vary in propensity to transmit and cause disease, but results are inconsistent and there is considerable strain-to-strain variation within lineages.

2Lineage-2 (Beijing) strains are associated with increasing spread and drug resistance in some areas but not others (European Concerted Action on New Generation Genetic Markers, 2006), and with a lower or higher proportion of extrapulmonary tuberculosis. M. africanum has been associated with lower virulence, and lineage-1 with faster sputum smear conversion. In low incidence settings, lineage is often associated with immigrant sub-groups, and while host–pathogen co-evolution has been suggested, it is difficult to disentangle the effects of lineage and host susceptibility on pathogenesis.

3Since the 1990s, methods such as RFLP based on the insertion element IS6110 have been used to distinguish clusters of patients with shared DNA-fingerprint patterns, suggesting recent transmission, but within the clusters, these methods cannot distinguish who transmitted to whom. Whole genome sequencing provides far greater resolution, and if data are collected in a whole population over several years, single nucleotide polymorphisms (SNPs) can be used to construct transmission networks. In low-incidence settings small numbers of SNPs have been found between epidemiologically linked patients, although the maximum SNP difference to ‘confirm’ a link is not yet established. No population-based study to-date has applied long-term large-scale whole genome sequencing in a high prevalence area, it is much more challenging to interpret transmission networks when there are many possible sources of infection. Yet understanding transmission in high prevalence areas would have the greatest public health benefit.

4As part of the Karonga Prevention Study in Malawi, we assess transmission using whole genome sequencing in the whole district over 15 years. We show decreasing transmission over time and marked variation between M. tuberculosis lineages 1–4 which are unconfounded by host differences.

5Between September 1995 and September 2010, there were 2332 person-episodes of culture-confirmed tuberculosis in Karonga District. Whole genome sequences that passed quality control were available for 1687 (72%). The distribution of patients with and without sequences available was very similar by age, sex, and HIV status. The proportion with sequences available was the highest in 2002–2006 (82%) and was higher in those with smear-negative pulmonary disease (78%) than in those with smear-positive disease (71%) and extra-pulmonary disease (67%).

6The phylogenetic tree is shown in Figure 1. Most M. tuberculosis strains (68%) were lineage-4, with 16% lineage-1, 4% lineage-2, and 12% lineage-3 (Table 1). Lineage-4 strains were more common in the earlier years. Lineage-1 strains were more common in HIV-positive and older patients and less common in recurrent tuberculosis. Lineage-2 strains were more common in younger patients and were all drug sensitive. Lineage-3 strains were associated with recurrent tuberculosis and with isoniazid resistance. There was no association between lineage and having been born or recently resident outside the district. The associations of lineage with HIV status and recurrent tuberculosis persisted after adjusting for age, sex, and year. The association between lineage and recurrent tuberculosis was also present when restricted to those with drug-sensitive strains, and the association between lineage and isoniazid resistance was also present when restricted to those with first episode tuberculosis.

7Figure 2 shows the SNP distances between all possible pairs of samples in the data set (including more than one per individual in some cases). Peaks corresponding to large numbers of SNPs represent comparisons between lineages. On the basis of the distribution, we chose cut-offs at 5 and 10 SNPs for distinguishing links. Similar figures were drawn for the mutation rate (Figure 2—figure supplement 1). We have previously shown that patients with relapse had up to 8 SNPs difference, and these cut-offs are similar to those used in other studies.

8To construct the transmission network, we included links of up to 10 SNPs difference. We included one sample per person-episode of disease and excluded extra-pulmonary cases as they cannot transmit. Example clusters are shown in Figure 3, and the full transmission network in Figure 3—figure supplement 1. Overall, after excluding relapses (recurrences with ≤10 SNPs difference from the initial episode), 66% of patients were in clusters with at least one other patient. Clusters ranged in size from 2 to 36 (Figure 4A), with 23% of patients in clusters of 10 or more. The size of the clusters varied by lineage (Figure 4B): compared to lineage-4 (the commonest lineage), lineage-2 and lineage-3 strains were more likely to be clustered and in larger clusters and lineage-1 strains were less likely to be clustered and were in smaller clusters. The median cluster size and interquartile range (IQR) for lineages 1–4 were 3 (1, 6), 13 (7, 24), 7 (2, 22), and 3 (1, 8), respectively. The p-values for differences between lineages were similar if non-clustered strains were excluded.

9Overall, of 824 links with 0–10 SNPs identified in the networks, 255 (31%) had 0 SNPs different, 182 (22%) had 1 SNP, 127 (15%) had 2 SNPs, 77 (9%) had 3 SNPs, 52 (6%) had 4 SNPs, 32 (4%) had 5 SNPs, and 99 (12%) had 6–10 SNPs different. The number of SNPs correlated with the time between disease onset in the pairs of individuals linked in the network (Figure 4C): linear regression r2 = 10%, p < 0.001. The regression coefficient suggests a mutation rate of 0.26 SNPs/year (95% CI 0.21–0.31). The regression results were the same if sputum collection dates were used instead of disease onset dates.

10The within-patient mutation rate was calculated in 74 individuals with multiple specimens, including 51 relapses, allowing ≤10 SNPs, and using the first and last specimens if there were more than two. The estimated mutation rate was 0.45 SNPs/year (95% CI 0.15–0.75), r2 = 11%, p = 0.004 (Figure 4—figure supplement 1).

11Figure 4D shows the number of SNPs in the likely transmissions identified from the network, by lineage. Lineage-2 had the lowest number of SNPs per transmission, and lineage 1 the highest. The median mutation rates per year for the different lineages were lineage-1, 0.58 (IQR 0.11–1.9); lineage-2, 0.11 (0–0.66); lineage-3, 0.35 (0–1.1); lineage-4, 0.40 (0–1.2) (p = 0.004, equality-of-medians test). The regression of number of SNPs by number of days showed no clear differences between lineages (Figure 4—figure supplement 2).

12We investigated the number of SNPs in the likely transmissions by smear status, HIV status, and isoniazid resistance of the initial and subsequent cases. There were no differences by the characteristics of the first case, but transmissions to smear-positive subsequent cases had slightly more SNPs than those to smear-negative subsequent cases (p = 0.05); and those to HIV-negative subsequent cases had slightly more SNPs than those to HIV positive subsequent cases (p = 0.02). Using mutation rates, the results were similar, but with smaller differences by smear status and HIV status of the subsequent case (p = 0.06 and 0.08, respectively).

13For further analysis of transmission, we excluded 77 uncertain links (i.e., with 6–10 SNPs and mutation rate ≥0.003 SNPS/day, Figure 2—figure supplement 1).

14A case of tuberculosis was defined as being due to recent infection if a source case was identified in the network within the previous 5 years, and not being due to recent infection if no source was identified or if the closest source (in terms of number of SNPs) was more than 5 years earlier. Overall, 38% of patients had evidence of recent infection (Table 2). This was the highest for lineage-2 (65%) and the lowest for lineage-1 (31%). Linkage with a recent source case was less common in older age groups, in those who had been living outside the district, and in more recent years, with the proportion linked decreasing from 45% in 1999–2001 to 30% in 2008–2010. These trends persisted after adjusting for each other (Table 2). There was no association of linkage with sex, HIV status, sputum smear status, or isoniazid resistance and adjusting for these did not affect the results. The effect of village of birth was lost after adjusting for recent residence.

15From the network, 32% of individuals were linked as likely sources of infection to at least one other individual. Individuals were sources for up to 12 others, with 293 (22%) linked to one, 76 (6%) linked to two, 22 (2%) linked to three, 14 (1%) linked to four, and 26 (2%) linked to five or more.

16Table 3 shows the association of characteristics of the index episode with the likelihood of transmission, using ordered logistic regression. There were more transmissions from those with positive smears, and with tuberculosis in the earlier years. Lineage-2 and lineage-3 strains were more likely to transmit than lineage-4, and these differences were more marked after adjustment for year, age, sex, and smear status. Place of birth and recent residence were weakly associated with onward transmission, and further adjusting for these or the other factors in the table did not affect the results.

17Comparing those with any transmissions vs those with none in a logistic regression model gave very similar results (not shown). Restricting the links to those within 3 years of the index episode, there was still a strong trend with year: the odds ratios from the ordered logistic regression analysis, adjusted for lineage, age, sex, and smear status, for the year groups 1999–2001, 2002–2004 and 2005–2007, compared to 1995–1998, were 0.47 (95% CI 0.33–0.66), 0.35 (0.25–0.50), and 0.37 (0.25–0.54), respectively.

18This is the largest whole genome sequencing study of M. tuberculosis transmission to-date, and the first to use a network approach. We show that this approach is feasible and that with long-term, population-wide data, important inferences can be made about transmission. In this population, although lineage-4 has been present for longer, lineages are not now associated with area of birth or recent residence, so differences by lineage are unlikely to be confounded by associations with host sub-populations.

19The mutation rates in this study are consistent with those from other settings and in vitro. This is the largest study to measure between-patient mutation rates. Although the confidence intervals on the estimate are narrow, there is considerable variation as others have found. The measure assumes the correct source has been identified and uses the time interval between dates of when the disease was first diagnosed or specimen collection as necessarily crude approximations of the time since divergence of the samples from their common ancestor. Furthermore, there is a bottleneck on transmission: most infections probably arise from one or very few organisms, which may be minority strains in the first case. The within-patient estimate of mutation rate does not have the same measurement problems and gave a consistent result.

20We showed striking differences between the lineages in cluster size, the proportion of disease due to recent transmission and transmissibility. Lineage-1 formed the smallest clusters, with the largest SNP differences. Patients with lineage-1 strains were the least likely to have disease due to recent transmission and were less likely to transmit and cause new cases than those with lineages 2 or 3. These observations suggest a lower propensity to cause disease, which may explain lineage-1's association with HIV infection, if it is less likely to cause active disease in those who are not immunosuppressed. Lineage-1 strains have been associated with lower virulence in animal models.

21Lineage-2 formed large clusters with small SNP differences. It had the highest proportion of disease due to recent transmission and the highest proportion of transmissions. Increased virulence in lineage-2 has been suggested previously, often in association with drug resistance, but in this population all lineage-2 strains were drug sensitive. Despite these associations that suggest higher virulence and transmissibility, the proportion of cases due to lineage-2 did not increase over the period. We have previously reported that lineage-2 was first detected in this area in 1991, initially increased, and then plateaued from around 2000. This may explain the lower proportion of lineage-2 strains in the oldest age group. The high proportion of linked cases with lineage-2 could reflect few imported (and therefore unlinked) cases, although there was no association between lineage and immigration.

22In contrast, lineage-3 increased as a proportion of tuberculosis cases over time. It was associated with an intermediate proportion of disease due to recent infection and high transmission. In this population, it is also associated with relapse. Lineage-4 had smaller cluster sizes than lineages 2 and 3. It remains the most common lineage in this population, although the proportion has fallen over time.

23Over the period of the study, the proportion of cases due to recent transmission decreased from 46% to 30%, and the proportion of cases transmitting and giving rise to new cases of tuberculosis also fell markedly. This correlates with a reduction in tuberculosis incidence over this period. It suggests a considerable success of the tuberculosis and HIV control programmes, despite the potential for M. tuberculosis transmission in antiretroviral clinics.

24We found no association with HIV infection in the proportion of disease due to recent infection (in contrast to our findings with RFLP in the earlier period) or in transmissibility. Social clustering of HIV-infected individuals may increase the opportunities for transmission to susceptible individuals who manifest disease, balancing out any decreased transmissibility. The change from our earlier findings could be due to the reduced transmission in the population and to the increasing use of isoniazid prophylaxis and antiretroviral therapy in HIV-positive individuals.

25In this study, we had high quality whole genome sequence data on 72% of culture-positive patients over 15 years. While this is a high proportion, links will be missed, and the best link found may not be the correct one (especially when there are multiple patients with identical strains). Missing links will lead to underestimation of the proportion of disease due to recent transmission and of transmissions. The missing and wrongly attributed links are likely to be randomly distributed, leading to non-differential misclassification of linkage, and underestimation of associations with lineage and other factors.

26This large, long-term study provides strong evidence for differences in transmission patterns and virulence between the M. tuberculosis lineages, particularly high transmissibility and virulence for lineages 2 and 3 and low transmissibility and virulence for lineage-1, which are unrelated to drug resistance, HIV infection, or host sub-population.

27In Karonga District, northern Malawi (population approximately 300,000), project staff at the hospital and peripheral health centres identify individuals with suspected tuberculosis, and sputum and other specimens are taken. All diagnosed tuberculosis patients are interviewed, and HIV-tested, after counselling and if consent is given. The incidence of new smear-positive tuberculosis in adults in the district has fallen from 124/100000/year to 87/100000/year over the period of this study, with about 6% isoniazid resistance and <1% multidrug resistance. Adult HIV prevalence in the area is around 10%.

28Approval for the study was given by the ethics committee of the London School of Hygiene & Tropical Medicine (#5067) and the Malawian National Health Sciences Research Committee (#424). Informed consent was obtained from all participants.

29Culture is performed in the project laboratories in Malawi, with species identification and drug susceptibility testing in the UK Mycobacterium Reference Laboratory. RFLP was performed on cultures from all patients from late 1995–2008. We processed all available stored DNA samples or cultures from 1995 to 2010 for whole genome sequencing at the Sanger Institute, using Illumina HiSeq 2000, paired-end reads of length 100 base-pairs.

30We used trimmomatic software (http://www.usadellab.org/cms/?page=trimmomatic) to remove low-quality reads and low-quality 3′ ends of reads, keeping only reads ≥50 base-pairs long, with nucleotides >Q27 (equivalent to a risk of error of <0.2% per read per base-pair).

31We mapped reads for each sample against the H37Rv reference genome (Genbank assession: AL123456.3), using the BWA-mem algorithm (http://bio-bwa.sourceforge.net/). We excluded samples with average genomic coverage less than 10-fold.

32We identified SNP positions using SAMtools (http://samtools.sourceforge.net/). Sample genotypes were called using the majority allele (minimum frequency 75%) in positions supported by at least 20-fold coverage; otherwise we classified them as missing (thus ignoring heterozygous calls). We excluded samples with >15% missing genotype calls, to remove possible contaminated or mixed samples or technical errors. (The proportion of mixed strains is low in this setting). We excluded genome positions with >15% missing genotypes, and those in highly repetitive and variable regions (e.g., PE/PPE genes).

33In the final analysis, 94% of the M. tuberculosis genome was analysed for variants. Median coverage was 88-fold, mean 127. Spoligotyping was performed in silico using SpolPred. Lineages were defined from spoligotype families.

34We calculated SNP distances between sequences using the ape library in the R statistical package (http://cran.r-project.org/). We computed a maximum-likelihood phylogenetic tree including all samples, using RAxML, using the GTRCAT model.

35For the transmission network, we used the SeqTrack package in R, using one sample per person-episode of disease, excluding episodes of extrapulmonary tuberculosis (as these cannot transmit). This builds a minimum-spanning tree, minimizing the genomic distance between links and keeping the disease onset dates coherent. Based on our data, we allowed up to 10 SNPs difference for inclusion in the networks. The suitability of the cut-off was assessed by examination of the SNP differences within and between patients. We have previously shown in 92 patients in this data set with repeat samples from the same or different episodes of disease that, using the same pipeline, there is a clear bimodal distribution, with pairs of samples either having up to 8 SNPs between them or more than 100 SNPs. Furthermore, among 187 pairs of individuals with epidemiological links, 62 had ≤10 SNPs, 9 had 10–99 SNPs, and 116 had ≥100 SNPs.

36Statistical analysis used STATA 13 (http://www.stata.com/). We estimated the between-patient mutation rate using linear regression of the number of SNPs by time between disease onset dates (taken as the date of first evidence of tuberculosis—the earliest of date of collection of the first positive sample, or registration or treatment) in the patients connected in the network. This analysis was repeated using dates of specimen collection. For comparison, we calculated the within-patient mutation rate, in individuals with more than one specimen from the same episode of disease or from a relapse, also using a cut-off of ≤10 SNPs.

37We compared the size of the clusters by lineage. Among the likely transmissions (≤10 SNPs), we examined the number of SNP differences and mutation rates by lineage, and characteristics of the index and subsequent case, using non-parametric tests (Wilcoxon rank-sum and the equality-of-medians test).

38For the analyses of risk factors for disease due to recently acquired infection and for transmissibility, we classified links with 6–10 SNPs different as uncertain unless the mutation rate was <0.003 SNPs/day, to allow for larger changes over long time periods. Those patients with uncertain links were excluded from the risk factor analyses.

39The SeqTrack network shows the most likely source of infection for each case. For groups of cases with zero SNPs between them the one closest in date was chosen. A case was defined as due to recently acquired infection if the most likely source was within 5 years, and not due to recent infection if there was no source identified or if the source was earlier than this (even if there were other closely related strains within 5 years).

40In this analysis of ‘backwards’ links, we used the first 3 years of data only for identifying previous links. We examined risk factors for disease due to recent infection among individuals with their first episode of tuberculosis, using logistic regression. The multivariable analysis included lineage, age, sex, and year a priori, and other factors if they were associated with recent infection after adjustment for these, or if they confounded other variables.

41The SeqTrack network links can also be used to examine forward transmission that results in disease. In this analysis, we used the last 3 years of data only to identify transmissions that had taken place, to allow time for transmissions causing new cases. We used ordered logistic regression to assess risk factors for transmission and the number of transmissions. In the multivariable analysis, we adjusted for lineage, age, sex, year, and sputum smear status of the index case a priori, and assessed confounding by other factors. We repeated the analysis using logistic regression, comparing any transmissions vs none. Since those in later years had less time for transmission to be detected, we examined the effect of calendar period using transmission within 3 years of the index case.

42Software sources for in-house programs will be made available on sourceforge.net (http://sourceforge.net/projects/patogenico/). Raw data can be obtained from the European Nucleotide Archive at EMBL-EBI (project accessions: ERP000436 and ERP001072).