Nimmo et al.

BMC Genomics (2019) 20:389

https://doi.org/10.1186/s12864-019-5782-2

RESEARCH ARTICLE Open Access

Whole genome sequencing Mycobacterium

tuberculosis directly from sputum identifies
more genetic diversity than sequencing
from culture
Camus Nimmo1,2* , Liam P. Shaw3,4, Ronan Doyle1,5, Rachel Williams1, Kayleen Brien2, Carrie Burgess1,
Judith Breuer1, Francois Balloux3 and Alexander S. Pym2

Abstract
Background: Repeated culture reduces within-sample Mycobacterium tuberculosis genetic diversity due to selection
of clones suited to growth in culture and/or random loss of lineages, but it is not known to what extent omitting
the culture step altogether alters genetic diversity. We compared M. tuberculosis whole genome sequences
generated from 33 paired clinical samples using two methods. In one method DNA was extracted directly
from sputum then enriched with custom-designed SureSelect (Agilent) oligonucleotide baits and in the other
it was extracted from mycobacterial growth indicator tube (MGIT) culture.
Results: DNA directly sequenced from sputum showed significantly more within-sample diversity than that
from MGIT culture (median 5.0 vs 4.5 heterozygous alleles per sample, p = 0.04). Resistance associated variants
present as HAs occurred in four patients, and in two cases may provide a genotypic explanation for phenotypic
resistance.
Conclusions: Culture-free M. tuberculosis whole genome sequencing detects more within-sample diversity than a
leading culture-based method and may allow detection of mycobacteria that are not actively replicating.
Keywords: Mycobacterium tuberculosis, Drug-resistant tuberculosis, Whole genome sequencing, Sputum, Within-patient
diversity, Heteroresistance

Background second line drug regimen for at least 9 months, which is

International efforts to reduce tuberculosis (TB) infec- also used for rifampicin monoresistance [3, 4]. With the
tions and mortality over the last two decades have only emergence of resistance to fluoroquinolones and amino-
been partially successful. In 2017, 10 million people de- glycosides (extensively drug-resistant [XDR] TB) there is
veloped TB and it has overtaken HIV as the infectious an increasing need for individualised therapy based on
disease responsible for the most deaths worldwide [1, 2]. drug susceptibility testing (DST). Individualised therapy
Drug resistance is a major concern with a steady rise in ensures patients are treated with sufficient active drugs
the number of reported cases globally and rapid in- which can prevent selection of additional resistance, im-
creases in some areas [1]. Patients with Mycobacterium prove treatment outcomes and reduce duration of infec-
tuberculosis resistant to the first line drugs rifampicin tiousness [5–8].
and isoniazid are classed as having multidrug-resistant Traditionally, phenotypic culture-based DST was used
(MDR) TB and usually treated with a standardised to identify drug resistance but this is being replaced by
rapid genetic tests that detect specific drug resistance-
* Correspondence: c.nimmo.04@cantab.net conferring mutations. Next generation whole genome se-
1
Division of Infection and Immunity, University College London, London
WC1E 6BT, UK quencing (WGS) of M. tuberculosis is being increasingly
2
Africa Health Research Institute, Durban, South Africa used in research and clinical settings to comprehensively
Full list of author information is available at the end of the article

© The Author(s). 2019, corrected publication 2019. Open Access This article is distributed under the terms of the Creative
Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use,
distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source,
provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain
Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article,
unless otherwise stated.
Nimmo et al. BMC Genomics (2019) 20:389 Page 2 of 9

identify all drug resistance associated mutations [9]. M. diversity in M. tuberculosis DNA sequenced directly
tuberculosis has a conserved genome with little genetic from sputum compared to DNA from culture [31],
diversity between strains and no evidence of horizontal whereas another study of mostly drug susceptible pa-
gene transfer [10], but more detailed analysis of individ- tients showed sequencing directly from sputum identi-
ual patient samples with WGS has identified genetically fied a slight excess of HAs relative to culture [33]. Here
separate bacterial subpopulations in sequential sputum we reanalyse heterozygous alleles (HAs) for the 12 avail-
samples [11–16] and across different anatomical sites able paired sequences with > 60-fold mean genome
[17]. This within-patient diversity can occur as a result coverage from that study [33] in addition to 21 newly
of mixed infection with genetically distinct strains or collected samples from patients with MDR-TB and
within-host evolution of a single infecting strain [18]. further explore the genomic location of the additional
Bacterial subpopulations can be detected in clinical diversity identified.
samples after sequencing reads are mapped to a refer-
ence genome where multiple base calls are detected at a Results
single genomic site. These heterozygous alleles (HAs) at Patient characteristics and drug susceptibility testing
sites associated with drug resistance (resistance associ- Whole genome sequences were obtained for 33 patients
ated variants, RAVs) may reflect heteroresistance, where from both mycobacterial growth indicator tube (MGIT)
a fraction of the total bacterial population is drug sus- culture and direct sputum sequencing. The patients were
ceptible while the remainder is resistant [19]. Identifica- predominantly of black African ethnicity (83%) and 50%
tion of genetic diversity within clinical samples may were HIV positive. First line phenotypic drug suscepti-
improve detection of RAVs over currently available rapid bility testing (DST) results identified 20 patients with
genetic tests [19] and can be achieved with freely avail- MDR-TB and one with rifampicin monoresistance. In
able WGS analysis toolkits [20–22]. Identifying RAVs addition there were two isoniazid monoresistant patients
could improve individualised therapy, prevent acquired and ethambutol resistance was detected in 7 patients.
resistance [12], and give insight into bacterial adaptation Second-line phenotypic DST was performed for patients
to the host. with rifampicin-resistant or MDR-TB and identified one
M. tuberculosis WGS is usually performed on fresh or case of kanamycin resistance (Table 1).
stored frozen cultured isolates to obtain sufficient puri- All samples had mean genome coverage of 60x or
fied mycobacterial DNA [23, 24]. However, the culture above with at least 85% of the genome covered at 20x
process can change the population structure from that (Additional file 1: Table S1). We observed greater mean
of the original sample due to genetic drift (random loss coverage depth in sputum-derived sequences than MGIT
of lineages) and/or the selection of subpopulations more sequences (median 173.7 vs 142.4, p = 0.03, Additional
suited to growth in culture [25–27]. Repeated subculture file 1: Table S1), and so mapped reads were randomly
leads to loss of genetic diversity and heteroresistance downsampled to give equal mean coverage depth in
[28]. Additionally, in the normal course of M. tubercu- each pair. A genotypic susceptibility profile was deter-
losis infection, some bacteria exist as viable non-cultur- mined by evaluating MGIT WGS for consensus-level
able persister organisms that are hypothesised to cause RAVs using a modified version of publicly available
the high relapse rate seen following treatment of insuffi- lists [22, 35]. Genotypic RAVs predicted all rifampicin
cient duration. Although these organisms may be identi- phenotypic resistance and > 95% of isoniazid phenotypic
fied in sputum by techniques such as reporter phages or resistance. Ethambutol genotypic RAVs were poorly pre-
culture with resuscitation promoting factors [29, 30] dictive of phenotypic resistance in line with findings from
they are likely to be missed by any sequencing method other studies [36] (Table 1). The patient with kanamycin
reliant on standard culture. phenotypic resistance was correctly identified by an rrs
WGS directly from sputum without enrichment is a1401g RAV. No full phenotypic fluoroquinolone pheno-
challenging [23]. It has recently been improved by de- typic resistance was identified, but several colonies from
pleting human DNA during DNA extraction [31]. We patient F1013 did grow in the presence of ofloxacin
have previously reported the use of oligonucleotide en- (although not enough to be classified as resistant). The
richment technology SureSelect (Agilent, CA, USA) to consensus sequences from this patient harboured a gyrB
sequence M. tuberculosis DNA directly from sputum E501D mutation which is believed to confer resistance to
[32] and demonstrated its utility in determining a rapid moxifloxacin but not other fluoroquinolones, which may
genetic drug resistance profile [33, 34]. explain the borderline phenotypic DST result [37].
It remains unclear to what extent WGS of cultured M.
tuberculosis samples underestimates the genetic diversity Genetic diversity
of the population in sputum samples. One previous To compare consensus sequences from sputum and
study of 16 patients did not identify increased genetic MGIT, a WGS consensus sequence-level maximum likelihood
Nimmo et al. BMC Genomics (2019) 20:389 Page 3 of 9

Table 1 Phenotypic and genotypic drug susceptibility testing (DST) results and sensitivity and specificity of genotypic DST relative
to phenotypic DST
Drug Resistance by phenotypic DST Resistance by genotypic DST Genotypic DST sensitivity Genotypic DST specificity
First line drugs
Rifampicin 21/32 (65.6%) 21/33 (63.6%) 21/21 (100%)a 21/21 (100%)
Isoniazid 22/32 (68.8%) 24/36 (66.7%) 21/22 (95.5%) 23/24 (95.8%)
Ethambutol 7/31 (22.6%) 15/34 (44.1%) 7/7 (100%) 7/15 (46.7%)
Second line drugs
Ofloxacin 0/22 (0.0%) 1/22 (4.5%) N/A 0/1 (0%)b
Kanamycin 1/22 (4.5%) 1/22 (4.5%) 1/1 (100%) 1/1 (100%)
Phenotypic DST available for first line drugs for 32 of the 33 patients, and for second line drugs for 22 patients who demonstrated rifampicin drug resistance
a
In one directly-sequenced sputum samples rifampicin RAVs were missed due to low coverage, although they were identified in the corresponding MGIT sample
b
This sample had < 1% of colonies grow in the presence of ofloxacin, so is categorised as susceptible but may have low-level or heteroresistance to
fluoroquinolones (see main text)

phylogenetic tree was constructed (Additional file 1: than one lineage or sublineage. We tested for mixed in-
Figure S1). As expected, all paired sequences were fection with the same sublineage by screening samples
closely related, with a median difference of 0.0 (range by HA frequency and then using Bayseian model based
0–1) single nucleotide polymorphisms (SNPs). Sam- clustering in samples with >10 HAs as described previ-
ples from patients F1066 and F1067 were closely re- ously [39]. This identified mixed infection in the sputum
lated with only one consensus-level SNP separating sample from patient F1096, which had 261 heterozygous
all four consensus sequences. There was no obvious alleles, greater than 10 times that in any other sample.
epidemiological link between these patients (although This patient was therefore excluded from further
this study was not designed to collect comprehensive analyses.
epidemiological information) and they lived 20 km As a first step to comparing diversity between sputum
apart in Durban. However, both patients were admit- and MGIT sequenced samples we looked at the location
ted contemporaneously to an MDR treatment facility of genetic diversity within the M. tuberculosis genome.
and sampled on the same day. DNA extraction and HAs were widely dispersed across the genome at similar
sequencing occurred on different runs. Therefore the sites in both sputum and MGIT samples. The genes with
close genetic linkage may represent direct transmis- the greatest density of HAs are shown in Table 2.
sion within a hospital setting, a community transmis- Notably, genetic diversity was found in the ribosomal
sion chain or an unlikely cross-contamination during RNA (rRNA) genes (rrs and rrl) uniquely in sputum
sample collection. samples, compared to other genes where distribution of
Having established congruence between sputum and diversity between MGIT and sputum was more balanced.
MGIT sequences at the consensus level we then com- As rRNA contains regions that are highly conserved
pared genetic diversity by DNA source. We first defined across bacteria [40], we considered the possibility that
a threshold for calling variants present as heterozygous
alleles (HAs) in our entire dataset by using a range of
minimum read count frequencies as described in the
methods (Fig. 1). Below a minimum of three supporting
reads there was an exponential increase in the number
of HAs identified, which may be indicative of the inclu-
sion of sequencing errors. To reduce this risk, we used a
threshold of a minimum of four supporting reads.
Genetic diversity may occur because of within-host
evolution or mixed infection. To identify mixed infection
we used a SNP-based barcode [38] to scan all HAs for a
panel of 413 robust phylogenetically informative SNPs
that can resolve M. tuberculosis into one of seven line-
ages and 55 sub-lineages. We found three phylogenetic
SNPs among the HAs. In all cases the heterozygous
phylogenetic SNP originated from the same sublineage Fig. 1 Variation in total number of heterozygous alleles (HAs) identified
across all 33 patients in sequences generated from sputum and MGIT
as other SNPs present at 100% frequency, and there
depending on minimum supporting read count threshold
were no cases of HAs indicating the presence of more
Nimmo et al. BMC Genomics (2019) 20:389 Page 4 of 9

Table 2 Genes with ≥2 heterozygous alleles (HAs) across all sputum samples, ordered by greatest number of HAs per base
Gene Heterozygous alleles per base Total number of heterozygous alleles Functional category
Sputum MGIT Sputum MGIT
rv1319c 0.021 0.021 33 33 Metabolism and respiration
rrs 0.016 0.000 25 0 16S ribosomal RNA
rrl 0.006 0.000 19 0 23S ribosomal RNA
ppsA 0.003 0.001 15 4 Lipid metabolism
rv2082 0.006 0.006 13 14 Unknown function
accE5 0.006 0.000 3 0 Lipid metabolism
lppB 0.005 0.005 3 3 Probable surface lipoprotein
pks12 0.000 0.001 3 10 Lipid metabolism
rv2319c 0.003 0.005 3 4 Stress protein
lppA 0.003 0.002 2 1 Probable surface lipoprotein
rpoC 0.001 0.001 2 3 RNA polymerase beta’ subunit
rv3888c 0.002 0.001 2 1 Probable membrane protein
vapC25 0.005 0.000 2 0 Possible toxin
vapC31 0.005 0.002 2 1 Possible toxin

SureSelect baits targeting rRNA genes were capturing explained by the selective nature of MGIT media which
both M. tuberculosis and other bacterial species. To will enrich M. tuberculosis sequences and the decontam-
evaluate this, metagenomic taxonomic assignment was ination step used to kill non-mycobacteria prior to cul-
performed on all reads by sampling reads that were not ture inoculation. Importantly the frequency of HAs in
assigned to M. tuberculosis (i.e. presumed contaminants other highly diverse genes between sequencing strategies
from other bacteria). We then performed a BLAST was more balanced (Table 2) in addition to the lack of
search against the most diverse genes listed in Table 2 BLAST hits of contaminating reads to these genes.
which indicated that a sizeable proportion of non-M. tu- After excluding the sample with mixed infection and
berculosis reads from directly sequenced sputum had a removing rRNA gene sequences we compared the fre-
BLAST hit of at least 30 bases to M. tuberculosis rrs and quency of HAs in sputum and MGIT. There were 265
rrl genes that encode rRNA (330 BLAST hits from spu- HAs identified across all sputum samples compared to
tum sequences vs 4 BLAST hits from MGIT sequences, 200 in MGIT samples (median 5.0 vs 4.5, p = 0.04,
median 8.5% vs 0.0%, p < 0.01, Additional file 1: Fig- Additional file 1: Table S1). In both sputum and MGIT
ure S2). There were no BLAST hits against any of samples, the majority of HAs were indels, and non-syn-
the other genes with ≥2 sputum HAs apart from onymous mutations were more commonly frameshift than
rpoC, for which there were 3 BLAST hits from spu- missense mutations (Table 3). The distribution of HAs by
tum sequences but none from MGIT sequences (me- patient is shown in Fig. 2.
dian 0.0% for both sputum and MGIT sequences),
indicating that this issue appears largely specific to Genetic diversity in drug resistance genes
rRNA. To determine if contaminating reads were HAs in drug resistance associated regions, including
contributing to HAs identified in intergenic regions, promoters and intergenic regions, were individually
we repeated this analysis for all intergenic regions assessed. Four of the 32 patients with single strain infec-
with ≥2 sputum HAs (Additional file 1: Table S2). tion had RAVs present as HAs in at least one gene,
There were no BLAST hits to any of these regions, which are shown in Table 4. Patient F1002 had three
suggesting that this is not the case. The taxonomic compensatory mutations in rpoC present at HAs in both
assignment of these contaminating reads were typical sequences. As described above, the strains from patients
of genera composing the oral flora, with a high representa- F1066 and F1067 were highly related with only one con-
tion of Actinomyces, Fusobacterium, Prevotella, and Strepto- sensus SNP difference between all four sequences. Both
coccus (Additional file 1: Figure S3). had phenotypic high level isoniazid resistance with no
This supported the hypothesis that the baits may en- consensus-level katG or inhA mutation, but had frame-
rich rRNA from other organisms so rRNA genes were shift katG mutations present as HAs which have the po-
excluded from further analysis. The difference in diver- tential to cause resistance [43]. F1066 and RF021 had
sity between sputum and MGIT sequences can be Rv1979c and pncA mutations respectively at low
Nimmo et al. BMC Genomics (2019) 20:389 Page 5 of 9

Table 3 Variants identified in MGIT and sputum derived genetic diversity between methods. All paired sequences
sequences from paired samples were closely related at the consensus level, and WGS
Sputum variants MGIT variants predicted phenotypic drug susceptibility with over 95%
Total variants 24,480 25,465 sensitivity and specificity for rifampicin and isoniazid in
Total variants present as HAs 265 (1.1%) 200 (0.8%) line with published data [44].
(% of total variants) We find that the rRNA genes have high levels of diver-
Median HAs per sample 5.0 4.5 sity in sputum samples, but believe this is due to
Variant type (% all HAs)
non-mycobacterial DNA hybridising to the capture baits.
This conclusion is borne out by the taxonomic assign-
SNP 217 (81.9%) 174 (87.0%)
ment of reads aligning to these genes in common oral
MNP 2 (0.8%) 0 (0.0%) bacteria. We therefore excluded these from further ana-
Insertion 4 (1.5%) 1 (0.5%) lysis, and recommend others using enrichment from
Deletion 24 (9.1%) 15 (7.5%) sputum do similarly. We find more diversity when se-
Complex 18 (6.8%) 10 (5.0%) quencing directly from sputum with significantly more
Coding change (% all HAs)
unique heterozygous alleles (HAs) than sequencing from
MGIT culture (p = 0.04).
Non-synonymous (missense) 93 (35.1%) 77 (38.5%)
The understanding of within-patient M. tuberculosis
Non-synonymous (frameshift) 6 (2.3%) 7 (3.5%) genetic diversity is becoming increasingly important as
Synonymous 57 (21.5%) 57 (28.5%) the detection of rare variants has been shown to im-
Intergenic 109 (41.1%) 59 (29.5%) prove the correlation between phenotypic and genotypic
Values given represent totals for 32 paired samples. SNP single nucleotide drug resistance profiles [19] and can identify emerging
polymorphism, MNP multi-nucleotide polymorphism drug resistance [11, 12]. Not including a culture step
avoids the introduction of bias towards culture-adapted
subpopulations and the impact of random chance and
frequency in sputum only which have the potential to is also likely to incorporate DNA from viable non-
confer phenotypic resistance to clofazimine (Rv1979c) culturable mycobacteria. A reduction in genetic diver-
and pyrazinamide (pncA), although no phenotypic test- sity has previously been shown with sequential M.
ing was performed for these drugs. tuberculosis subculture [25, 28], but was not con-
firmed by a study performing WGS directly from spu-
Discussion tum [31]. However, the 16 paired sputum and MGIT
In this study we performed whole genome sequencing samples compared by Votintseva [31] had a minimum
using DNA from sputum and MGIT culture in paired of 5x coverage compared to a minimum 60x coverage
samples from 33 patients and compared within-patient in this study, and were likely to contain less genetic

Fig. 2 Number of heterozygous alleles (HAs) found in directly sequenced sputum only (sputum), MGIT (MGIT) only or in both samples (shared) by patient
Nimmo et al. BMC Genomics (2019) 20:389 Page 6 of 9

Table 4 Resistance associated variants present as heterozygous alleles (HAs)

Patient ID Phenotypic resistance Mutation Frequency (MGIT/sputum) Description
F1002 Rifampicin rpoB S450 L 100%/100% High confidence resistance mutation
F1002 Rifampicin rpoC G332R [41] 82.6%/21.7% Putative compensatory mutations
F1002 Rifampicin rpoC L516P [41] 12.7%/7.7%
F1002 Rifampicin rpoC P1040S [42] 21.7%/12.3%
F1066 Isoniazid (high) katG N218 fs 0.0%/6.9% Possible resistance mutations, not
previously described
F1066 Clofazimine – not tested Rv1979c G376D 0.0%/0.5%
F1067 Isoniazid (high) katG N218 fs 10.7%/7.6%
RF021 Pyrazinamide – testing failed pncA Q122H 0%/2.5%

material as they were surplus clinical rather than ded- Conclusions

icated research samples. Directly sequencing M. tuberculosis from sputum is able
Two-thirds of the patients with MDR-TB had already to identify more genetic diversity than sequencing from
been treated for drug susceptible-TB (DS-TB), and add- culture. Characterising within-patient genetic diversity is
itional diversity in sputum samples may represent early important to understand bacterial adaptation to drug
adaptation to drug pressure. As direct sputum sequen- treatment and the acquisition of drug resistance. It also
cing does not rely on live mycobacteria, DNA from has potential to identify low frequency RAVs that may
recently killed M. tuberculosis is likely to also be se- further enhance the prediction of drug resistance pheno-
quenced, meaning that recent genomic mutations are type from genotype.
likely to be represented as HAs.
In two patients, RAVs present as HAs provided a likely Methods
genotypic basis for otherwise unexplained phenotypic re- Patient enrolment
sistance. Given the small total number of resistance muta- Adult patients presenting with a new diagnosis of spu-
tions in this study, it is not possible to draw conclusions tum culture positive TB were included in the study.
about the frequency of heterozygous RAVs in directly se- Patients were recruited in London, UK (n = 12) and
quenced sputum. However the presence of heterozygous Durban, South Africa (n = 21). All patients recruited in
RAVs in both MGIT and sputum sequences reinforces the Durban were Xpert MTB/RIF (Cepheid, CA, USA) posi-
biological importance of these mutations. tive for rifampicin resistance. Two sputum samples were
To reduce the risk of sample cross contamination, collected prior to starting the current treatment regi-
paired samples were extracted on different days, pre- men, with one inoculated into mycobacterial growth in-
pared in different sequencing libraries and sequenced dicator tube (MGIT) culture (BD, NJ, USA) and the
on different runs. However it is not possible to com- other used for direct DNA extraction. Therefore for pa-
pletely exclude the possibility of contamination during tients with drug susceptible-TB (DS-TB), sputum was
sample collection and between different samples proc- collected prior to taking any TB therapy, while patients
essed in batches. A further limitation of this study is starting MDR-TB treatment may have already taken
that it can be difficult to distinguish low frequency treatment for DS-TB if this was intiated prior to resist-
variants from sequencing error. The SureSelect library ance results being available.
preparation protocol for sputum sequencing incorpo-
rates more PCR cycles than that used for MGIT se- Ethics, consent and permissions
quencing, which may increase the risk of error. All patients gave written informed consent to participate
Where possible this could be evaluated further by in the study. Ethical approval for the London study was
performing technical sequencing replicates on ex- granted by NHS National Research Ethics Service East
tracted DNA samples, although this was not possible Midlands–Nottingham 1 (reference 15/EM/0091). Ethical
due to insufficient surplus material and financial con- approval for the Durban study was granted by University
straints. To reduce the risk of sequencing errors we of KwaZulu-Natal Biomedical Research Ethics Committee
used high read and mapping quality thresholds, and (reference BE022/13).
required a stringent 98% identity between sequenced
reads and the reference genome. Low frequency vari- Microbiology
ants of particular clinical importance could be con- MGIT samples were incubated in a BACTEC MGIT 960
firmed by resequencing the same DNA samples. (BD, NJ, USA) until flagging positive. Phenotypic DST
Nimmo et al. BMC Genomics (2019) 20:389 Page 7 of 9

data for London samples were those provided to treating Variant calling
hospitals by Public Health England. Phenotypic DST Variant calling for comparison for HA counts was per-
were performed using equivalent standardised methods. formed with FreeBayes v1.2 [51]. Variants falling in or
For Durban samples this was the solid agar proportion within 50 bases of PE/PPE family genes and repeat ele-
method (Additional file 1: Methods) and for London ments were excluded using vcfinteresect in vcflib [52].
samples the resistance ratio method [45]. For the initial analysis of genetic diversity, variants were
included if supported by ≥2 reads, with ≥1 forward and
reverse read, no read position bias, a minimum mapping
DNA extraction and sequencing
quality of 30 and base quality of 30. The minimum
Positive MGIT tubes were centrifuged at 16,000 g for 15
supporting read threshold was then increased in a
min and the supernatant removed. Cells were resus-
stepwise fashion from 2 to 15. Variant calling files
pended in phosphate-buffered saline before undergoing
where variants were supported ≥4 supporting reads
heat killing at 95 °C for 1 h followed by centrifugation at
including ≥1 forward and reverse read were used to
16,000 g for 15 min. The supernatant was removed and
compare HA frequency and location and to screen
the sample resuspended in 1 mL sterile saline (0.9% w/
for mixed infection.
v). The wash step was repeated. DNA was extracted with
The phylogenetic tree was constructed by calling vari-
mechanical ribolysis before purification with DiaSorin
ants with VarScan v2.4.0 [53] mpileup2cns as this is able
Liaison Ixt (DiaSorin, Italy) or CTAB [46]. NEBNext
to generate consensus-level calls at each reference se-
Ultra II DNA (New England Biolabs, MA, USA) was
quence base. SNPs were then used to generate a se-
used for DNA library preparation.
quence of equal length to the reference using a custom
Sputum samples for direct sequencing were heat
perl script and these sequences were combined in a
killed, centrifuged at 16,000 g for 15 min and the super-
multi-alignment fasta file. SNP sites were extracted from
natant was removed. DNA extraction was performed
this alignment using snp-sites v2.4.1 [54], and pairwise
with mechanical ribolysis followed by purification using
SNP differences calculated using snp-dists v0.6.3 [55].
DiaSorin Liaison Ixt (DiaSorin, Italy) or DNeasy blood &
Extracted SNP sites were used to generate a maximum
tissue kit (Qiagen, Germany) [46]. Target enrichment
likelihood phylogenetic tree using RaxML v8.2.12 [56]
was performed using SureSelect with a custom-designed
which was visualised using FigTree v1.4.3.
bait set covering the entire positive strand of the M. tu-
berculosis genome as described previously [33]. Batches
of 48 multiplexed samples were sequenced on NextSeq Identification of mixed infection
500 (Illumina, CA, USA) 300-cycle paired end runs with All samples were screened for evidence of mixed infec-
a mid-output kit. Sequencing was performed by the tion using described methods [39]. In brief, any sample
Pathogen Genomics Unit at University College London with 10 or fewer heterozygous SNPs, or between 11 and
in a dedicated laboratory where one sequencing run was 20 heterozygous SNPs where heterozygous SNPs were ≤
processed at a time. All paired samples were extracted, 1.5% of all SNPs was classified as not mixed. For other
prepared and sequenced on different days. The National samples, the Baysian mixture model analysis [39] was
Center for Biotechnology Information Sequence Read used where samples with a Bayesian information criter-
Archive (NCBI SRA) accession number for each sample ion value > 20 for presence of more than one strain were
is shown in Additional file 1: Table S3. assumed to be mixed.

Read mapping Metagenomic assignment

DNA sequence reads were adapter and quality trimmed Sequencing reads were classified using Kraken v0.10.6
then aligned to the H37Rv reference genome (GenBank [57] against a custom Kraken database previously con-
accession NC_000962.3) with Trim Galore v0.4.4 [47] structed from all available RefSeq genomes for bacteria,
and BBMap v38.32 [48], with mapped reads stored in an archaea, viruses, protozoa, and fungi, as well as all
output bam file. Duplicate reads were removed with RefSeq plasmids (as of September 19th 2017) and three
Picard tools v1.130 [49] MarkDuplicates and coverage human genome reference sequences [58]. The size of the
statistics generated with Qualimap v2.2.1 [50]. For each final database after shrinking was 193 Gb, covering
sample pair, the bam file with greater mean genome 38,190 distinct NCBI taxonomic IDs.
coverage was randomly downsampled to that of the To assess the proportion of contaminating reads
paired sample with Picard tools v1.130 [49] Downsam- that could generate spurious diversity when mapped
pleSam. All further analyses were performed using these to M. tuberculosis ribosomal genes, we randomly sub-
downsampled bam files. Command line parameters used sampled 100 reads taxonomically assigned as non-M.
are specified in the Additional file 1: Methods. tuberculosis and performed a BLAST search with
Nimmo et al. BMC Genomics (2019) 20:389 Page 8 of 9

blastn v2.2.28 [59] against each gene as described Genetics Institute, University College London, London WC1E 6BT, UK.
4
from the H37Rv reference genome. We only analysed Nuffield Department of Clinical Medicine, Oxford University, Oxford OX3
7BN, UK. 5Clinical Research Department, London School of Hygiene and
hits of at least 30 bases. Tropical Medicine, London WC1E 7HT, UK.

Statistics Received: 3 October 2018 Accepted: 7 May 2019

Statistical analyses were performed with Prism v8.0

(Graphpad, CA, USA). Mean coverage depth statistics,
References
number of HAs and BLAST hits of contaminating reads 1. Global Tuberculosis Report 2018. Geneva: World Health Organization; 2018
in paired samples were compared using a two-tailed 2018.
Wilcoxon matched-pairs signed rank test. 2. Murray CJ, Ortblad KF, Guinovart C, Lim SS, Wolock TM, Roberts DA, et al.
Global, regional, and national incidence and mortality for HIV, tuberculosis,
and malaria during 1990-2013: a systematic analysis for the global burden
Additional file of disease study 2013. Lancet. 2014;384(9947):1005–70.
3. Dheda K, Gumbo T, Maartens G, Dooley KE, McNerney R, Murray M, et al.
The epidemiology, pathogenesis, transmission, diagnosis, and management
Additional file 1: Supplemetary Methods, Figures and Tables. (PDF 490
of multidrug-resistant, extensively drug-resistant, and incurable tuberculosis.
kb)
Lancet Respir Med. 2017;5(4):291–360.
4. WHO treatment guidelines for drug-resistant tuberculosis. World Health
Abbreviations Organization; 2016.
DST: Drug susceptibility testing; DS-TB: Drug susceptible-tuberculosis; 5. Trauner A, Liu Q, Via LE, Liu X, Ruan X, Liang L, et al. The within-host
HA: Heterozygous allele; MDR-TB: Multidrug resistant-tuberculosis; population dynamics of mycobacterium tuberculosis vary with treatment
MGIT: Mycobacterial growth indicator tube; RAV: Resistance associated efficacy. Genome Biol. 2017;18(1):71.
variant; rRNA: Ribosomal RNA; SNP: Single nucleotide polymorphism; 6. Olaru ID, Lange C, Heyckendorf J. Personalized medicine for patients with
TB: Tuberculosis; WGS: Whole genome sequencing MDR-TB. J Antimicrob Chemother. 2016;71(4):852–5.
7. Pasipanodya JG, McIlleron H, Burger A, Wash PA, Smith P, Gumbo T. Serum
Acknowledgements drug concentrations predictive of pulmonary tuberculosis outcomes.
The authors would like to thank Sashen Moodley, Ashentha Govender and J Infect Dis. 2013;208(9):1464–73.
Colin Chetty in the Microbiology Core, Africa Health Research Institute. 8. Cegielski JP, Kurbatova E, van der Walt M, Brand J, Ershova J, Tupasi T, et al.
Multidrug-resistant tuberculosis treatment outcomes in relation to
Funding treatment and initial versus acquired second-line drug resistance. Clin
Camus Nimmo is funded by a Wellcome Trust Research Training Fellowship Infect Dis. 2016;62(4):418–30.
reference 203583/Z/16/Z. Judith Breuer receives funding from the UCL/UCLH 9. Satta G, Lipman M, Smith GP, Arnold C, Kon OM, McHugh TD.
NIHR funded Biomedical Research Centre. This work was additionally funded Mycobacterium tuberculosis and whole-genome sequencing: how close are
by National Institute for Health Research via the UCLH/UCL Biomedical we to unleashing its full potential? Clin Microbiol Infect. 2018;24(6):604–9.
Research Centre (grant number BRC/176/III/JB/101350) and the PATHSEEK 10. Sreevatsan S, Pan X, Stockbauer KE, Connell ND, Kreiswirth BN, Whittam TS,
European Union’s Seventh Programme for research and technological et al. Restricted structural gene polymorphism in the mycobacterium
development (grant number 304875). The funding bodies had no input on tuberculosis complex indicates evolutionarily recent global dissemination.
study design, analysis, data interpretation or manuscript writing. Proc Natl Acad Sci U S A. 1997;94(18):9869–74.
11. Sun G, Luo T, Yang C, Dong X, Li J, Zhu Y, et al. Dynamic population
Availability of data and materials changes in mycobacterium tuberculosis during acquisition and fixation
Original fastq files are available at NCBI Sequence Read Archive with BioProject of drug resistance in patients. J Infect Dis. 2012;206(11):1724–33.
reference PRJNA486713: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA486713/ 12. Merker M, Kohl TA, Roetzer A, Truebe L, Richter E, Rüsch-Gerdes S, et al.
Whole genome sequencing reveals complex evolution patterns of
Authors’ contributions multidrug-resistant mycobacterium tuberculosis Beijing strains in patients.
Study conception: JB, ASP, Data collection: CB, KB, Analysis and interpretation: PLoS One. 2013;8(12):e82551.
CN, LPS, RD, RW, Drafting of manuscript: CN, LPS, Revision of manuscript: FB, JB, 13. Operario DJ, Koeppel AF, Turner SD, Bao Y, Pholwat S, Banu S, et al.
ASP, Final approval of manuscript: CN, LPS, RD, RW, KB, CB, JB, FB, ASP. Prevalence and extent of heteroresistance by next generation sequencing of
multidrug-resistant tuberculosis. PLoS One. 2017;12(5):e0176522.
Ethics approval and consent to participate 14. Black PA, de Vos M, Louw GE, van der Merwe RG, Dippenaar A, Streicher
All patients gave written informed consent to participate in the study. Ethical EM, et al. Whole genome sequencing reveals genomic heterogeneity and
approval for the London study was granted by NHS National Research Ethics antibiotic purification in mycobacterium tuberculosis isolates. BMC
Service East Midlands–Nottingham 1 (reference 15/EM/0091). Ethical approval Genomics. 2015;16(1):857.
for the Durban study was granted by University of KwaZulu-Natal Biomedical 15. Eldholm V, Norheim G, von der Lippe B, Kinander W, Dahle UR, Caugant DA,
Research Ethics Committee (reference BE022/13). et al. Evolution of extensively drug-resistant mycobacterium tuberculosis from
a susceptible ancestor in a single patient. Genome Biol. 2014;15(11):490.
Consent for publication 16. Bloemberg GV, Keller PM, Stucki D, Stuckia D, Trauner A, Borrell S, et al.
Not applicable Acquired resistance to Bedaquiline and Delamanid in therapy for tuberculosis.
N Engl J Med. 2015;373(20):1986–8.
Competing interests 17. Lieberman TD, Wilson D, Misra R, Xiong LL, Moodley P, Cohen T, et al.
The authors declare that they have no competing interests. Genomic diversity in autopsy samples reveals within-host dissemination of
HIV-associated mycobacterium tuberculosis. Nat Med. 2016;22(12):1470–4.
18. Ford C, Yusim K, Ioerger T, Feng S, Chase M, Greene M, et al. Mycobacterium
Publisher’s Note tuberculosis--heterogeneity revealed through whole genome sequencing.
Springer Nature remains neutral with regard to jurisdictional claims in Tuberculosis (Edinb). 2012;92(3):194–201.
published maps and institutional affiliations. 19. Metcalfe JZ, Streicher E, Theron G, Colman RE, Allender C, Lemmer D, et al.
Cryptic micro-heteroresistance explains M. tuberculosis phenotypic
Author details resistance. Am J Respir Crit Care Med. 2017;196(9):1191–201.
1
Division of Infection and Immunity, University College London, London 20. Bradley P, Gordon NC, Walker TM, Dunn L, Heys S, Huang B, et al. Rapid
WC1E 6BT, UK. 2Africa Health Research Institute, Durban, South Africa. 3UCL antibiotic-resistance predictions from genome sequence data for
Nimmo et al. BMC Genomics (2019) 20:389 Page 9 of 9

Staphylococcus aureus and mycobacterium tuberculosis. Nat Commun. 41. Yang C, Luo T, Shen X, Wu J, Gan M, Xu P, et al. Transmission of multidrug-
2015;6:10063. resistant mycobacterium tuberculosis in Shanghai, China: a retrospective
21. Feuerriegel S, Schleusener V, Beckert P, Kohl TA, Miotto P, Cirillo DM, et al. observational study using whole-genome sequencing and epidemiological
PhyResSE: a web tool delineating mycobacterium tuberculosis antibiotic investigation. Lancet Infect Dis. 2017;17(3):275–84.
resistance and lineage from whole-genome sequencing data. J Clin 42. Eldholm V, Monteserin J, Rieux A, Lopez B, Sobkowiak B, Ritacco V, et al.
Microbiol. 2015;53(6):1908–14. Four decades of transmission of a multidrug-resistant mycobacterium
22. Coll F, McNerney R, Preston MD, Guerra-Assunção JA, Warry A, Hill-Cawthorne tuberculosis outbreak strain. Nat Commun. 2015;6:7119.
G, et al. Rapid determination of anti-tuberculosis drug resistance from whole- 43. Heym B, Alzari PM, Honore N, Cole ST. Missense mutations in the catalase-
genome sequences. Genome Med. 2015;7(1):51. peroxidase gene, katG, are associated with isoniazid resistance in
23. Doughty EL, Sergeant MJ, Adetifa I, Antonio M, Pallen MJ. Culture-independent mycobacterium tuberculosis. Mol Microbiol. 1995;15(2):235–45.
detection and characterisation of mycobacterium tuberculosis and M. 44. Coll F, Phelan J, Hill-Cawthorne GA, Nair MB, Mallard K, Ali S, et al. Genome-
africanum in sputum samples using shotgun metagenomics on a wide analysis of multi- and extensively drug-resistant mycobacterium
benchtop sequencer. PeerJ. 2014;2:e585. tuberculosis. Nat Genet. 2018;50(2):307–16.
24. Bjorn-Mortensen K, Zallet J, Lillebaek T, Andersen AB, Niemann S, Rasmussen 45. Sam IC, Drobniewski F, More P, Kemp M, Brown T. Mycobacterium
EM, et al. Direct DNA extraction from mycobacterium tuberculosis frozen tuberculosis and rifampin resistance, United Kingdom. Emerg Infect Dis.
stocks as a Reculture-independent approach to whole-genome sequencing. 2006;12(5):752–9.
J Clin Microbiol. 2015;53(8):2716–9. 46. Larsen MH, Biermann K, Tandberg S, Hsu T, Jacobs WR. Genetic manipulation
25. Depledge DP, Palser AL, Watson SJ, Lai IY, Gray ER, Grant P, et al. Specific of mycobacterium tuberculosis. Curr Protoc Microbiol 2007;Chapter 10:Unit
capture and whole-genome sequencing of viruses from clinical samples. 10A.2.
PLoS One. 2011;6(11):e27805. 47. Krueger F. TrimGalore. Available from: https://github.com/FelixKrueger/
26. Hanekom M, Streicher EM, Van de Berg D, Cox H, McDermid C, Bosman M, TrimGalore . [cited 2019 Mar 19].
et al. Population structure of mixed mycobacterium tuberculosis infection is 48. Bushnell B. BBMap. Available from: https://jgi.doe.gov/data-and-tools/
strain genotype and culture medium dependent. PLoS One. 2013;8(7): bbtools/. [cited 2019 Mar 19].
e70178. 49. Picard Tools: Broad Institute. Available from: http://broadinstitute.github.io/
27. Martin A, Herranz M, Ruiz Serrano MJ, Bouza E, Garcia de Viedma D. The picard/. [cited 2019 Mar 19].
clonal composition of mycobacterium tuberculosis in clinical specimens 50. Okonechnikov K, Conesa A, Garcia-Alcalde F. Qualimap 2: advanced multi-
could be modified by culture. Tuberculosis (Edinb). 2010;90(3):201–7. sample quality control for high-throughput sequencing data. Bioinformatics.
28. Metcalfe JZ, Streicher E, Theron G, Colman RE, Penaloza R, Allender C, et al. 2016;32(2):292–4.
Mycobacterium tuberculosis subculture results in loss of potentially clinically 51. Garrison E, Marth G. Haplotype-based variant detection from short-read
relevant heteroresistance. Antimicrob Agents Chemother. 2017;61(11): sequencing. arXiv preprint arXiv:12073907. 2012.
e00888–17. 52. Garrison E. Vcflib. Available from: https://github.com/vcflib/vcflib. [cited 2019
29. Jain P, Weinrick BC, Kalivoda EJ, Yang H, Munsamy V, Vilcheze C, et al. Dual- Mar 19].
reporter Mycobacteriophages (Phi2DRMs) reveal preexisting mycobacterium 53. Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, et al. VarScan
tuberculosis persistent cells in human sputum. MBio. 2016;7(5):e01023–16. 2: somatic mutation and copy number alteration discovery in cancer by
30. Mukamolova GV, Turapov O, Malkin J, Woltmann G, Barer MR. Resuscitation- exome sequencing. Genome Res. 2012;22(3):568–76.
promoting factors reveal an occult population of tubercle bacilli in sputum. 54. Page AJ, Taylor B, Delaney AJ, Soares J, Seemann T, Keane JA, et al. SNP-
Am J Respir Crit Care Med. 2010;181(2):174–80. sites: rapid efficient extraction of SNPs from multi-FASTA alignments. Microb
31. Votintseva AA, Bradley P, Pankhurst L, Del Ojo Elias C, Loose M, Nilgiriwala K, Genom. 2016;2(4):e000056.
et al. Same-day diagnostic and surveillance data for tuberculosis via whole 55. Seemann T. snp-dists. Available from: https://github.com/tseemann/snp-
genome sequencing of direct respiratory samples. J Clin Microbiol. 2017; dists. [cited 2019 Mar 19].
55(5):1285–98. 56. Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-
32. Brown AC, Bryant JM, Einer-Jensen K, Holdstock J, Houniet DT, Chan JZ, et analysis of large phylogenies. Bioinformatics. 2014;30(9):1312–3.
al. Rapid whole-genome sequencing of mycobacterium tuberculosis isolates 57. Wood DE, Salzberg SL. Kraken: ultrafast metagenomic sequence classification
directly from clinical samples. J Clin Microbiol. 2015;53(7):2230–7. using exact alignments. Genome Biol. 2014;15(3):R46.
58. Lassalle F, Spagnoletti M, Fumagalli M, Shaw L, Dyble M, Walker C, et al. Oral
33. Doyle RM, Burgess C, Williams R, Gorton R, Booth H, Brown J, et al. Direct
microbiomes from hunter-gatherers and traditional farmers reveal shifts in
whole genome sequencing of sputum accurately identifies drug resistant
commensal balance and pathogen load linked to diet. Mol Ecol. 2018;27(1):
mycobacterium tuberculosis faster than MGIT culture sequencing. J Clin
182–95.
Microbiol. 2018;56(8):e00666–18.
59. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al.
34. Nimmo C, Doyle R, Burgess C, Williams R, Gorton R, McHugh TD, et al. Rapid
BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.
identification of a mycobacterium tuberculosis full genetic drug resistance
profile through whole genome sequencing directly from sputum. Int J
Infect Dis. 2017;62:44–6.
35. Consortium CR. The GP, Allix-Beguec C, Arandjelovic I, bi L, Beckert P, et al.
prediction of susceptibility to first-line tuberculosis drugs by DNA
sequencing. N Engl J Med. 2018;379(15):1403–15.
36. Walker TM, Kohl TA, Omar SV, Hedge J, Del Ojo Elias C, Bradley P, et al.
Whole-genome sequencing for prediction of mycobacterium tuberculosis
drug susceptibility and resistance: a retrospective cohort study. Lancet
Infect Dis. 2015;15(10):1193–202.
37. Malik S, Willby M, Sikes D, Tsodikov OV, Posey JE. New insights into
fluoroquinolone resistance in mycobacterium tuberculosis: functional
genetic analysis of gyrA and gyrB mutations. PLoS One. 2012;7(6):e39754.
38. Coll F, McNerney R, Guerra-Assuncao JA, Glynn JR, Perdigao J, Viveiros M, et
al. A robust SNP barcode for typing mycobacterium tuberculosis complex
strains. Nat Commun. 2014;5:4812.
39. Sobkowiak B, Glynn JR, Houben R, Mallard K, Phelan JE, Guerra-Assuncao JA,
et al. Identifying mixed mycobacterium tuberculosis infections from whole
genome sequence data. BMC Genomics. 2018;19(1):613.
40. Konstantinidis KT, Tiedje JM. Prokaryotic taxonomy and phylogeny in the
genomic era: advancements and challenges ahead. Curr Opin Microbiol.
2007;10(5):504–9.