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

To improve understanding of the factors influencing tuberculosis transmission and the role of pathogen variation, we sequenced all available specimens from patients diagnosed over 15 years in a whole district in Malawi. Mycobacterium tuberculosis lineages were assigned and transmission networks constructed, allowing ≤10 single nucleotide polymorphisms (SNPs) difference. We defined disease as due to recent infection if the network-determined source was within 5 years, and assessed transmissibility from forward transmissions resulting in disease. High-quality sequences were available for 1687 disease episodes (72% of all culture-positive episodes): 66% of patients linked to at least one other patient. The between-patient mutation rate was 0.26 SNPs/year (95% CI 0.21–0.31). We showed striking differences by lineage in the proportion of disease due to recent transmission and in transmissibility (highest for lineage-2 and lowest for lineage-1) that were not confounded by immigration, HIV status or drug resistance. Transmissions resulting in disease decreased markedly over time. DOI: http://dx.doi.org/10.7554/eLife.05166.001


Introduction
Despite 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 (Homolka et al., 2010). 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) (Comas et al., 2009), and one intermediate (lineage-7), recently described in Ethiopia (Firdessa et al., 2013). The lineages may vary in propensity to transmit and cause disease (Thwaites et al., 2008;Homolka et al., 2010;Parwati et al., 2010;Gagneux, 2012), but results are inconsistent and there is considerable strain-to-strain variation within lineages (Portevin et al., 2011;Mathema et al., 2012).
Lineage-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 (Click et al., 2012) or higher (Kong et al., 2007) proportion of extrapulmonary tuberculosis. M. africanum has been associated with lower virulence (de Jong et al., 2008), and lineage-1 with faster sputum smear conversion (Click et al., 2013). 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 (Reed et al., 2009;Gagneux, 2012;Pareek et al., 2013).
Since the 1990s, methods such as RFLP based on the insertion element IS6110 (van Embden et al., 1993) have been used to distinguish clusters of patients with shared DNA-fingerprint patterns, suggesting recent transmission (Small et al., 1994), 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 (Bryant et al., 2013;Walker et al., 2013Walker et al., , 2014. In low-incidence settings small numbers of SNPs have been found between epidemiologically linked patients (Kato-Maeda et al., 2013), although the maximum SNP difference to 'confirm' a link is not yet established (Perez-Lago et al., 2014). No population-based study to-date has applied long-term large-scale whole genome sequencing in a high prevalence area (Luo et al., 2014;Walker et al., 2014), 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.
As 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.

Results
Between September 1995 and September 2010, there were 2332 person-episodes of cultureconfirmed 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%). eLife digest Tuberculosis is an important public health threat around the globe and is particularly common in developing countries. It is difficult to control the spread of the disease because the bacteria that cause it can spread when an infected individual coughs or sneezes. It may take years for an infected individual to develop symptoms of tuberculosis so it can be hard to trace the source of an outbreak, and people infected with HIV are particularly susceptible to the disease.
The bacterium that causes the majority of cases of tuberculosis is called Mycobacterium tuberculosis. There are several different varieties or 'lineages' of M. tuberculosis, and it is thought that they may vary in their ability to spread and cause disease. However, the results of previous studies have been inconsistent and there also seems to be a lot of variation between strains within the same lineage.
In this study, Guerra-Assunç ã o et al. used an approach called whole genome sequencing alongside more traditional methods to study the spread of tuberculosis in Malawi. They sequenced the genomes of every available sample of M. tuberculosis collected from patients in the Karonga district of Malawi over a 15-year period. This produced high-quality DNA sequence data about the bacteria responsible for almost 1700 cases of disease.
Using this massive amount of data, Guerra-Assunç ã o et al. constructed networks that showed how the bacteria had spread in the community. This revealed that there were differences between the ability of the various M. tuberculosis lineages to cause disease and to spread in communities. For example, lineage 1 was less likely than the other lineages to cause disease soon after infecting an individual and was less able to spread.
The data also show that the proportion of cases of disease due to recent infection declined substantially during the 15-year period. This indicates that the tuberculosis and HIV control programmes in the area have been successful.
Guerra-Assunç ã o et al.'s findings show that it is possible to understand how tuberculosis is transmitted on a large scale. The next challenge is to understand why the lineages differ in their ability to cause disease and spread between individuals.
The 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.   Figure 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 (Guerra-Assuncao et al., 2014), and these cut-offs are similar to those used in other studies (Bryant et al., 2013;Walker et al., 2013).

Transmission network
To 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 nonclustered strains were excluded.
The regression results were the same if sputum collection dates were used instead of disease onset dates. The 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), r 2 = 11%, p = 0.004 (Figure 4-figure supplement 1). Figure 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).
We 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).
For further analysis of transmission, we excluded 77 uncertain links (i.e., with 6-10 SNPs and mutation rate ≥0.003 SNPS/day,

Recent infection
A 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,   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.

Transmissibility
From 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. Table 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.
Comparing 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

Discussion
This 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.
The mutation rates in this study are consistent with those from other settings (Bryant et al., 2013;Walker et al., 2013) and in vitro (Ford et al., 2013). This is the largest study to measure betweenpatient 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 (Lurie, 1964), 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.
We 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 (Narayanan et al., 2008;Reiling et al., 2013).  Lineage-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 (Parwati et al., 2010), 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 (Glynn et al., 2005a. 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.
In 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.
Over 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 (Mboma et al., 2013). It suggests a considerable success of the tuberculosis and HIV control programmes, despite the potential for M. tuberculosis transmission in antiretroviral clinics.
We 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 [Houben et al., 2009) 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.
In 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.
This 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.

Patients
In Karonga District, northern Malawi (population approximately 300,000), project staff at the hospital and peripheral health centres identify individuals with suspected tuberculosis ), In this analysis individuals are defined as linked ('backwards links') using the cut-offs described in the text and if the closest link was with a patient within the previous 5 years. Extrapulmonary, recurrent cases, and cases before 1999 were excluded. Odds ratios (OR) calculated using logistic regression. *In this model a dummy variable was used for the 32 individuals with missing data on recent residence. †Test for trend. DOI: 10.7554/eLife.05166.012 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 (Mboma et al., 2013). Adult HIV prevalence in the area is around 10%. Approval 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.

Cultures and sequencing
Culture is performed in the project laboratories in Malawi, with species identification and drug susceptibility testing in the UK Mycobacterium Reference Laboratory (Mboma et al., 2013). RFLP was performed on cultures from all patients from late 1995-2008 (Glynn et al., 2005b). 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.

Read quality filtering
We used trimmomatic software (http://www.usadellab.org/cms/?page=trimmomatic) to remove lowquality 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).
We identified SNP positions using SAMtools (http://samtools.sourceforge.net/) (Li et al., 2009). 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).
We 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.

Transmission mapping
For the transmission network, we used the SeqTrack package in R (Jombart et al., 2011), 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 (Guerra-Assuncao et al., 2014). Furthermore, among 187 pairs of individuals with epidemiological links, 62 had ≤10 SNPs, 9 had 10-99 SNPs, and 116 had ≥100 SNPs. Statistical 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.
We 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-ofmedians test).
For 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.

Disease due to recently acquired infection
The 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).
In 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.

Transmissibility
The 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.

Repositories for data and software
Software 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).