The non-clonality of drug resistance in Beijing-genotype isolates of Mycobacterium tuberculosis from the Western Cape of South Africa

Selection of Strains for Sequencing

M. tuberculosis HN878 was obtained from the NIH TB Vaccine Testing and Research Materials Laboratory at Colorado State University. Fourteen additional drug-resistant strains of the Beijing genotype were selected for whole-genome sequencing and comparative analysis. The isolates were sampled from a database maintained at the DST/NRF Centre, Stellenbosch University, representing drug resistant TB cases from a variety of hospitals in the Western Cape, South Africa (see Table 1). The samples were all selected to be isoniazid-resistant, and were chosen to span a range of degrees of drug resistance, including cases that are mono-resistant (isoniazid), MDR (isoniazid and rifampicin), pre-XDR (MDR plus resistance to either a fluoroquinolone or an aminoglycoside), and XDR. No other biases (e.g. age, gender, HIV status) were applied in the sampling. All 14 isolates have the Beijing spoligotype, 00000000003771. Six of the isolates investigated in this study are from sublineage 1, representing the atypical Beijing genotype (ancestral, family F31). IS6110-RFLP analyses have classified these strains as being members of cluster R86. The remaining 8 isolates are from cluster R220 in sublineage 6, which represents "typical" Beijing genotype strains (more recently evolved). The IS6110 RFLP patterns used for determining cluster membership of these isolates are shown in Figure 1.

Strain ID	Origin (hospital)	Sub-lineage	Family	Cluster	DR type	Combined resistance profiles determined by DST	Combined sensitivity profiles determined by DST
R1207	George	1	F31	R86	MDR	INH, RIF	EMB
X132	Khayelitsha	1	F31	R86	pre-XDR	INH,RIF,AMI,CAP,STR,KAN	ETH,OFL,EMB
X28	Dysselsdorp	1	F31	R86	XDR	INH,RIF,AMI,CAP,OFL,STR, KAN	ETH, EMB
R1746	Mossel Bay	1	F31	R86	MDR	INH, RIF	EMB
X156	Brooklyn Chest	1	F31	R86	pre-XDR	INH,RIF, AMI,CAP,STR,KAN	ETH,OFL,EMB
X85	George	1	F31	R86	XDR	INH,RIF, AMI,OFL,KAN,ETH	CAP,EMB
R1909	Worcester	6	F29	R220	MDR	INH,RIF,EMB
R1842	George	6	F29	R220	Mono	INH	RIF, EMB
X122	Huguenot	6	F29	R220	pre-XDR	INH,RIF,OFL	ETH,AMI,EMB
R1390	Stellenbosch	6	F29	R220	Mono	INH	RIF,EMB
X189	Victoria Hos.	6	F29	R220	XDR	INH,RIF,AMI,CAP, OFL,KAN,	ETH,STR
R1505	George	6	F29	R220	MDR	INH, RIF	EMB
R1441	George	6	F29	R220	Mono	INH	RIF,EMB
X29	Retreat	6	F29	R220	pre-XDR	INH,RIF,AMI,STR,KAN	EMB,OFL,CAP,ETH

INH = isoniazid RIF = rifampicin STR = streptomycin EMB = ethambutol ETH = ethionamide AMI = amikacin KAN = kanamycin CAP = capreomycin OFL = ofloxacin

Drug-Susceptibility Testing

Drug susceptibility testing for isoniazid (0.2 μg/ml), rifampin (1 μg/ml), ethambutol (7.5 μg/ml), ethionamide (20 μg/ml), amikacin (30 μg/ml) and ofloxacin (2 μg/ml) was performed using the indirect proportion method on Middlebrook's medium [36]. Additional drug susceptibility testing was done using the proportion method on Middelbrook's 7H10 medium containing ofloxacin, amikacin, kanamycin and capreomycin at the critical concentration of 2 μg/ml, 5 μg/ml, 5 μg/ml and 10 μg/ml, respectively.

DNA Preparation and Sequencing Reaction

Sequencing of the genomes of 15 Beijing strains of M. tuberculosis was carried out on an Illumina Genome Analyzer II (Illumina, Inc.). In this study, Illumina Paired-End sequencing method (PE) was used. The cetyltrimethylammonium bromide (CTAB)-lysozyme method was used for extraction and purification of genomic DNA [37]. Suitable genomic DNA samples for the GAII were prepared as described on the sample preparation protocol (Illumina). 2-3 μg of genomic DNA was used for sample preparation. Genomic DNA was sheared by a nebulizer to generate DNA fragments for the Illumina Paried-End Sequencing method. The specific oligonucleotides (Illumina adapters) designed for PE sequencing were ligated to both ends of DNA fragments with the TA cloning method. Adapter-ligated DNA fragments of length 350-400 bp were isolated from a 2% agarose gel (Certified low-range Ultra Agarose, BIO-RAD) by using QIAquick Gel Extraction Kit (QIAGEN). Then the fragments were amplified by PCR reaction to generate the DNA library (15-30 ng/μl). The median size of the library was estimated by examining the 2% agarose gel image. The molarity of the DNA library was estimated as described in the sample preparation protocol (Illumina). The DNA libraries (5 picomoles each), including the ϕX174 control (bacteriophage DNA), were loaded on the flow cell for the cluster generation and sequencing. 79 cycles of images were collected, representing two 36-bp reads (paired ends) plus 7 bases representing id tags for multiplexing. The images were analyzed using version 1.4 of the GAPipeline software supplied by Illumina, producing files with 10-20 million pairs of 36 bp reads for each genome.

Sequence Determination and Bioinformatics

The reads were analyzed by comparative genome assembly to determine the complete sequence of each genome using custom software developed in our lab. The reads were first mapped (aligned) against the genome of H37Rv as a reference sequence. The mapping of each read was accomplished by ungapped alignment to the reference genome (including the reverse complement strand), allowing at most two mismatches out of 36 nucleotides. Initially, the reads that were pair-mates were mapped independently; subsequently, locations of reads for which the paired-end did not match within 300 bp were discarded. The mapped reads were used to assemble a list of the nucleotides observed at each position within the reference genome contributed by all the reads that overlapped it. Base calls were made by a maximum likelihood calculation, computed as the product of the probabilities for each base at each position using uncertainties estimated during image analysis. Sites where apparent differences were observed were subjected to local contig-building, in order to determine whether the difference was due to a nucleotide substitution or a small insertion or deletion.

For each genome, a list of verified differences was prepared and used to modify the reference genome to produce an intermediate ('edited') genome. Then the process was repeated by re-mapping the reads against the edited genome, and re-calling nucleotides at each position. For any sites that still had zero coverage, the base from the reference strain was used. This primarily included regions with exceptionally high GC-content (80-90%).

Large-scale deletions were identified by analyzing paired-end data for reads whose paired-end maps an unusually long distance away (i.e. >300 bp, compared to average read length of ~200). Large-scale insertions were determined by building contigs that spanned fragments, neither of whose paired ends mapped into the H37Rv genome. These ~200 bp fragments were then assembled into larger contigs using Newbler (software from Roche, inc.), localized in the genome using paired-end information, and identified by performing BLAST searches on the NCBI website.

Virtual spoligotyping was performed by aligning (without gaps) all the reads obtained for each strain against each of the 43 spacer sequences (26-bp oligos) from the direct repeats (DR) regions [38]. The number of matching reads for each spacer was counted, considering both forward and reverse-complement sequences, and accepting up to 1 nucleotide mismatch. Spacers with 0 matches were interpreted as missing.

Whole-genome alignments of the Beijing sequences to other mycobacterial strains were generated using MUMMER version 3.20 [39] to identify maximal stretches of perfectly matching regions, selecting optimal order-preserving assignments of matches between the genomes based on the longest increasing subsequence algorithm [40], and then using in-house sequence-alignment tools to determine alignments of the intervening regions with mismatches and/or gaps.

The complete genome sequences for M. tuberculosis strains HN878, R1207 (representative of R86 cluster), and X122 (representative of R220 cluster) have been deposited in GenBank with accession numbers ADNF01000000, ADNH01000000, and ADNG01000000, respectively.

Whole-Genome Sequencing

In order to establish a complete genome sequence to use as a reference for sequencing of other members of the Beijing strain family, HN878 was chosen as a representative strain and the genome was sequenced using Solexa sequencing technology. HN878 is a fully drug-susceptible member of the modern (typical) W-Beijing family which was isolated in Houston, TX as part of a TB outbreak in the 1990's, and has been used as a reference strain in multiple studies to characterize virulence and other properties of the Beijing strain [5, 7, 8]. HN878 was sequenced in paired-end mode using 36 bp reads. Mapping of reads to unique spacer sequences in the direct-repeats (DR) region confirms that HN878 has the characteristic 000000000003771 spoligotype associated with the Beijing strain family. The genome sequence of M. tuberculosis H37Rv (NCBI accession: NC_000962) was used as a reference sequence for comparative assembly of the HN878 genome. The mean depth of coverage was 70.8x, and 98.4% of the genome was covered by at least one read. The remaining uncovered regions are primarily associated with the family of PGRS genes, and are likely due to inefficient sequencing of GC-rich regions. In HN878, 91.5% of sites with zero coverage (60,454/66,065) were located in PGRS genes.

Among 1546 SNPs relative to H37Rv, 1296 occur in protein-coding regions. Only 1063 of 3989 genes contain a polymorphism of any type, and the remaining 73% of the genes in the genome are identical with the sequence in H37Rv. 75 genes contain frameshift mutations (see Additional file 1, Table S1). No mutations typically associated with drug resistance were found in the following genes, consistent with the pan-susceptible phenotype of this clinical isolate: inhA, katG, ethA, gyrA, iniABC, kasA, ndh, rpoB, rpsL, rrs, pncA, and embB. HN878 has the katG R643L and gyrA S95T alleles, confirming its membership in principle genetic group 1 [4]. Further details on SNPs and indels in HN878 (Table S2) relative to H37Rv, including IS6110 transpositions (Table S3) are described in the Additional file 1.

There was a high concordance of IS6110 insertion sites between R220 strains and HN878 (see Figure 2). Of the 21 insertion sites in HN878, four are absent in R220, although an additional four novel insertion sites are present in this cluster (see Additional file 1, Table S4). Similar to HN878, R220 has a single insertion in the NTF region, grouping it with HN878 in the modern sublineage of Beijing strains [34]. The R220 strains lack the RD150 deletion characteristic of sublineage 7 [33, 41], having Rv1671-Rv1674 intact and placing them in sublineage 6. The six R86 strains were all found to have the Beijing spoligotype and a common set of 14 IS6110 insertion sites (Additional file 1, Table S5), including one insertion in the dnaA-dnaN region. However, only 6 of these sites were shared with HN878. The R86 strains do not have any insertion in the NTF region (~3.48 Mb), classifying them as "ancestral" or "atypical" Beijing lineages [34].

Table 2 shows that the number of SNPs and indels was lower between the R220 strains and HN878 compared to the R86 strains, suggesting they are more closely related to HN878. The number of SNPs and indels compared to HN878 were 665-799 and 49-61, respectively, for the cluster R86 strains, and 267-392 and 35-54, respectively, for the cluster R220 strains.

The clinical isolates in the R220 cluster are fairly homogeneous at the genomic level. Among 1573 polymorphic sites between the R220 strains and H37Rv (SNPs not involving PPE genes, PGRS genes, or repetitive elements), 1234 differences (78.4%) were shared among HN878 and all 8 R220 strains. Some of these differences from H37Rv might be shared with other strain families like LAM. There were 113 sites (7.2%) where the R220 strains shared a difference from H37Rv but HN878 did not, and there were 91 sites (5.8%) where HN878 differed from H37Rv but the R220 strains did not. Each of the strains except one had a small number of unique SNPs (1-12).

R86 strains are found to have a set of SNPs compared to H37Rv that only partially overlaps with the SNPs in HN878. Out of a selected set of 1885 SNPs, 1126 (59.7%) were found to be common among HN878 and all six R86 strains relative to H37Rv, HN878 has 297 (15.8%) unique SNPs, and the R86 strains have 256 (13.6%) shared SNPs not in HN878. Thus, while both HN878 and the R86 share many polymorphisms presumably acquired in a common ancestor of the Beijing family, they have each diverged significantly, reflecting the difference in their lineages (sublineage 1 for R86 versus sublineage 6 for HN878).

When SNPs at synonymous sites are categorized into specific base-pair substitutions, the majority of replacements (~60%) are found to be G:C to A:T transversions (68/112 in the R86 cluster, and 75/126 in the R220 cluster) (Figure 3). However, a similar bias is observed among synonymous SNPs in CDC1551 or F11 (186/374 combined) when compared to H37Rv (χ² = 7.36 < 11.07 for the distribution of SNPs in the R220 cluster versus those in CDC1551/F11, which means the difference in substitution patterns is not significant at the p = 0.05 confidence level, df = 5), suggesting that the effect is not due to a Beijing-specific shift in the types of mutation.

Drug-Resistance Mutations

Rifampicin resistance is explained in three R86 strains (R1207, X132, and X28) by the D435V mutation in rpoB, while the other three R86 strains (R1746, X156, and X85) have S450L (equivalent to amino acids 516 and 531 in the conventional rpoB numbering based on E. coli). Both mutations are known to confer resistance to rifampicin [46]. S450L (TCG->TTG) is the most frequently observed mutation in cultures selected on rifampicin in vitro, implying lowest fitness cost [32]. D435V (GAC->GTC) is much less common in vitro (~1%) [47], but is frequently observed clinically [46, 47]. Four of the five R220 isolates that are rifampicin-resistant (R1909, X122, X189, and R1505) have a mutation in rpoB (either S450L or H445Y). However, isolate X29 (rifampicin-resistant) does not have any mutations in rpoB, and conversely, isolate R1842 (rifampicin-sensitive) has S450L in rpoB, which in turn may reflect laboratory error. The H445Y mutation is less frequent among clinical isolates, but also occurs in the rifampicin binding site in RNA polymerase.

Resistance to aminoglycosides can result from a variety of polymorphisms, generally related to the small subunit of the ribosome. Mutations that confer resistance to streptomycin and similar compounds are usually found in the 500-bp or 900-bp region of rrs, the 16S rRNA, or in rpsL, one of the ribosomal proteins [48]. All six R86 strains have the a514c mutation in rrs, conferring streptomycin resistance. X28 and X156 showed streptomycin resistance in drug-susceptibility testing; data on the other four strains was not available. Kanamycin binds a different site on the ribosome and shows cross-resistance with amikacin and capreomycin [49]. The most common mutation associated with kanamycin resistance is a1401g in rrs[50]. All R86 strains except isolate R1746 have the a1401g mutation in rrs. Four of the six strains were kanamycin-resistant (the other two were not tested). The capreomycin sensitivity of isolate X85 is inconsistent with the kanamycin/amikacin resistance, and might be an anomaly in drug-susceptibility testing. In the R220 cluster, isolates X189 and X29 are resistant to kanamycin and amikacin, and both strains show the expected a1401g mutation in rrs, whereas X122 (amikacin-sensitive) does not. R1842 has rrs c517t and X122 has K43R in ribosomal protein rpsL, each of which should confer streptomycin resistance, though this was not tested.

Resistance to fluoroquinolones is typically caused by mutations in gyrA[51]. Isolates X28 and X85 were reported to be ofloxacin-resistant, and X132 and X156 were reported to be ofloxacin-sensitive (data for R1207 and R1746 was not available). Three different amino acid mutations were observed in gyrA among these six strains: D94H in isolate R1207, D94N in X28, and D94G in X85. The SNPs in isolates X28 and X85 explain their ofloxacin resistance, and the fact that X132 and X156 have wild-type gyrA sequences is consistent with their ofloxacin sensitivity. The ofloxacin resistance of isolates X122 and X189 can be explained by mutations in D94 in gyrA. Isolate X29, which is sensitive to oflaxacin, has a wild-type gyrA sequence.

All six R86 strains have the mutation M306I in embB, commonly associated with resistance to ethambutol [52], although all six strains were reported to be sensitive to ethambutol by drug-susceptibility testing. It is well known that drug-susceptibility testing under-reports ethambutol resistance [53]. Note that, while all 6 strains have a mutation from Met to Ile, they use different codons; 'ATA' (for R1207, X132, and X28) and 'ATC' (for R1746, X156, and X85). Many of the R220 strains have a mutation in Met306 in embB, either to Val or Ile, although most of the strains tested sensitive to ethambutol as well.

Although resistance to pyrazinamide was not tested, 10 of the 14 strains showed mutations in pncA (pyrazinamidase). One R86 isolate (X29) has a frameshift mutation of 8 bp in pncA, while another (X189) had a large deletion of 881 bp knocking out the entire coding region, along with the adjacent gene, Rv2044c. Strain R1909 has mutation D8N, and strain X122 has Y103* (truncation mutation). All six R220 strains show mutations in pncA: three strains have a frameshift mutation (+c in codon 172), and three have amino acid substitution (C14R). Multiple mutations, including frameshifts, throughout pncA are associated with resistance to pyrazinamide [54–56].

Non-clonal Acquisition of Drug-Resistance Mutations

In order to better understand the origins of drug resistance among the Beijing strains, a phylogenetic tree was constructed to identify their evolutionary relationships. Subsequently, drug resistance mutations were mapped onto the tree to determine whether the patterns of drug resistance (associations between strains) could be explained by the same topology (i.e. inheritance). A master set of 727 polymorphic sites (SNPs only) with good depth of coverage (≥10x, i.e. sites covered by at least 10 reads) across all sequenced strains was selected. A subset of 704 sites was produced by removing those involved in drug resistance (e.g. inhA, katG, gyrA, rpoB, pncA, rrs, embB). The 704 sites were used to construct a maximum parsimony tree using dnapars in PHYLIP 3.66 [57]. Figure 4 shows the phylogenetic relationship with representative branch lengths, clearly showing the tight clustering of the two groups, cluster R86 (X85, X156, R1746, X132, X28, R1207) and the R220 cluster (R1441, R1505, R1390, X29, R1842, X189, R1909, X122), indicating that both groups have evolved independently from a common ancestor.

The number of SNPs differentiating the strains within each cluster is small but non-zero. Since the relationships are difficult to see in Figure 4a, the phylogeny was re-drawn as a cladogram (without meaningful branch lengths) in Figure 4b, and each branch is labeled with the number of SNPs associated with it. The observation that cluster R86 isolates share 236 SNPs while the R220 isolates share 303 unique SNPs suggests that these two clusters are largely clonal, and that isolates within the clusters are differentiated from one another by a limited number of SNPs. For example, the similarity between X29 and R1842 is supported by two unique SNPs: Q219P in Rv2571c and V16V in proW. The SNPs associated with each branch point are listed in Additional file 1, Table S7. It is also interesting to note that some individual isolates have continued to diverge, for example isolate R1746 has accumulated 12 unique polymorphisms. These differentiating SNPs likely represent recent evolutionary events. According to this phylogeny, HN878 is more closely related to the R220 cluster (sharing 195 SNPs, compared to 0 shared with the R86 cluster, using H37Rv as the outgroup), but is distinguished from the R220 members by a further 195 unique differences (108 + 87). It is important to note that these polymorphisms do not necessarily imply anything about the overall population structure, due to the non-random nature of the sample. Nonetheless, they represent the diversity and inter-relationships among the 14 clinical isolates sequenced.

When mutations related to drug resistance were projected onto this phylogeny, the patterns were frequently found to disagree with the phylogeny above. The most parsimonious explanation for this observation is that mutations conferring resistance to a given drug have arisen multiple times independently within both the R220 and R86 clusters. Figure 5 shows the pattern of mutations in gyrA conferring fluoroquinolone resistance. In these 14 Beijing strains, mutations in Asp94 appear to have arisen 6 times independently, each associated with a unique isolate (i.e. they were not clustered in the cladogram). In fact, Asp94 is mutated to 4 different amino acids (Gly: X85, X122, Asn: X28, His: R1207, Ala: X29, R1842, X189), again supporting the uniqueness of these events.

Mutations in rpoB related to rifampicin resistance are shown in Figure 6. In the R86 cluster, the six isolates divide into two groups: three isolates with the S450L mutation, and three with the D435V mutation. Similarly, the S450L mutation also explains the rifampicin resistance in a group of five isolates of the R220 cluster, clearly demonstrating that each cluster has acquired rifampicin resistance independently. It is important to note that three of the eight R220 isolates did not bear the S450L mutation. Although many R220 isolates have the S450L polymorphism [58], as of 2006 only 73% of R220 isolates were RIF-resistant [27], indicating that this polymorphism had not yet achieved fixation in the population, and thus the population structure has not yet been taken over by a distinct MDR clone. Of the three R220 isolates without the S450L mutation, two had wild-type rpoB sequences and were RIF-susceptible, and the third had a distinct RIF-resistance mutation: H445Y.

With respect to isoniazid resistance, Figure 7 shows that the c-15t inhA promoter mutation arose independently in both the R220 and R86 clusters. Whereas the c-15t mutation is found in all 8 strains of R220, it is only found in 2 of the 6 strains in the R85 cluster. However, the g-17t promoter mutation occurs in three other strains in the R86 cluster. The c-15t and g-17t inhA promoter mutations in the R86 cluster both occurred in the context of the katG S315T mutation, which occurs on an earlier branch and is found in all the strains of the R86 cluster. Simultaneous mutations in katG and the inhA promoter are frequently observed in clinical isolates [42], although it is unexpected that the inhA promoter mutations, which confer lower-level resistance to isoniazid, apparently occurred second. In contrast, the Ile194Thr mutation in inhA in isolates X189 and R1909 appears to have been acquired subsequent to the c-15t inhA promoter mutation common to the R220 strains, probably conferring higher-level resistance.

When the mutations in rrs and rpsL are mapped onto the phylogeny of the Beijing isolates (Figure 8), it is found that the R86 strains all shared mutation a516g in rrs, suggesting that streptomycin resistance arose in a common ancestor. The rrs a1401g mutation, responsible for kanamycin resistance, was observed in 5 of the 6 cluster R86 isolates, but not in isolate R1746. Therefore, it most likely has arisen twice within this cluster (unless there was a reversion to wild-type in R1746). In the R220 cluster, mutations (and most likely aminoglycoside resistance) were less prevalent. Isolate R1842 has a unique c517t rrs mutation that presumably also causes streptomycin resistance, and X122 is the only strain that has the K43R mutation in rpsL.

Finally, among the 14 Beijing isolates (both clusters), there appear to be six distinct mutations of Met306 in embB (Figure 9).