Genetic signatures of Mycobacterium tuberculosis Nonthaburi genotype revealed by whole genome analysis of isolates from tuberculous meningitis patients in Thailand

Selection of strains

Three isolates, CSF3053, 46-5069 and 43-13838, identified to belong to Nonthaburi genotype by IS6110-RFLP, were selected from the stock of samples collected from the CSF of TBM patients at the Drug Resistant Tuberculosis Research Fund Laboratory, Department of Microbiology, Faculty of Medicine Siriraj Hospital, Mahidol University, Thailand.

Genomic DNA extraction

Stock culture of selected strains, stored at −70 °C in MH79 broth containing 15% glycerol, were subcultured on Loewenstein-Jensen medium and incubated for 4 weeks at 37 °C. DNA extraction was carried out using cetyltrimethylammonium bromide (CTAB)-lysozyme enzymatic method as earlier described (Larsen et al., 2007).

Spoligotyping

Spacer oligonucleotide typing, a polymerase chain reaction (PCR) based method used in typing M. tuberculosis was performed following the methods earlier described (Gori et al., 2005).

Whole genome sequencing and analysis

Genomic DNA samples isolated from the three isolates were sequenced at Macrogen Inc., Seoul, South Korea on the HiSeq 2000 platform with insert size of 300 bp (Illumina, San Diego, CA, USA) yielding 100 bp paired end reads. The qualities of the sequences were assessed with FastQC software (www.bioinformatics.babraham.ac.uk/projects/fastqc) to determine the parameters used for trimming. Bases with quality of less than 5, reads with average of quality less than 20 for every four bases, and reads with lengths that are less than 45 bases were discarded using Trimmomatic software (Bolger, Lohse & Usadel, 2014) (version 0.33). The trimmed sequences were aligned to the reference strain M. tuberculosis H37Rv (NC_000962.3) using the short reads aligner, Bowtie2 (version 2.2.0) (Langmead & Salzberg, 2012). The genomic coverage was estimated using Bedtools (version 2.18) (Quinlan & Hall, 2010). The fold coverage is estimated as the number of reads supporting a particular nucleotide position on the genome. Variant calling was performed on the aligned sequences using the Genome Analysis Tool Kit (GATK) (version 3.3) haplotype caller (McKenna et al., 2010) with minimum calling confidence threshold set at phred score 30. Point allelic variation at any position within the genome when compared with the reference H37Rv genome (NC_000962.3) is considered a single nucleotide polymorphism (SNP).

SnpEFF (Cingolani et al., 2012) (version 4.0) software was used to annotate the SNPs. The SNPs were filtered using standard hard filtering parameters according to GATK Best Practices Recommendations (DePristo et al., 2011; Van der Auwera et al., 2013). Variants with QualByDepth <2.0, FisherStrand >60, RMSMapping quality <40, MappingQualityRankSumTest <−12.5 and ReadPosRankSumTest <−8 were filtered. All SNPs were confirmed using Integrated Genomic Viewer (IGV) (James et al., 2011) (version 2.0). The SNPs were further grouped according to the functions of the genes in which they were found in the genome when compared to the reference genome H37Rv (NC_000962.3). We evaluated SNPs in groups of genes considered to be essential, drug resistance related, virulence related, contain known epitopes and associated with efflux pumps.

The Whole Genome Shotgun project has been deposited at DDBJ/EMBL/GenBank under the accession numbers LGCH00000000, LGCG00000000and LGCF00000000. The versions described in this paper are LGCH01000000, LGCG01000000and LGCF01000000for CSF3053, 46-5069 and 43-13838 respectively. The raw sequences have been deposited to the short read archive (SRA) of NCBI under accession numbers SRX1094547, SRX1094546and SRX1094545for isolates CSF3053, 46-5069 and 43-13838 respectively.

Determination of principle genetic group, lineage and sequence type

Nucleotide alleles at positions 7585 and 2154724 were investigated to determine the principal genetic group of the isolates as earlier defined (Sreevatsan et al., 1997). To determine the lineage of the isolates, SNPs specific to different lineages as earlier reported (Coll et al., 2014b) were investigated.

Draft genome assembly

The paired-end raw reads of the isolates were assembled into draft genomes by using the de novo assembly algorithm of CLC Genomics Workbench (version 7.5) which works by using a de Bruijn graph (http://www.clcbio.com/). The minimum contig output was set at 200 bp long. Annotation of the draft genome was performed by Rapid Annotation using Subsystem Technology (RAST) (http://www.nmpdr.org/) and by NCBI Prokaryotic Genome Annotation Pipeline (PGAP) (http://www.ncbi.nlm.nih.gov/genome/annotation_prok/).

Comparison of Nonthaburi isolates with isolates from other lineages

The SNPs that are common to the three isolates were compared with 92,000 SNPs from 1,601 genomes of MTBC previously reported (Coll et al., 2014a) (http://pathogenseq.lshtm.ac.uk/phytblive/index.php). These include 121, 390, 189, 856, 17, 11, and 6 genomes from lineages 1, 2, 3, 4, 5, 6, and 7 respectively. Eleven samples from M. bovis were also included.

Large sequence polymorphism determination

Regions of differences when compared with reference strain H37Rv (NC_000926.3) were determined by using the indel and structural variants determination tool of CLC Genomics Workbench (version 7.5) (http://www.clcbio.com) and Bedtools (version 2.18) (Quinlan & Hall, 2010). The regions of deletions were confirmed with PCR using primers CF (CATCCGCACCGAACCTGTAA) and CR (AACCGTTCACGACAAGCAAC), AF (GCCCAACCTGATTGGTTTCG) and AR (CAAACGCTCGCCATGATCTC), BF (TCGACTGCCATACAACCTGC) and BR (ACTTCCGGTGGTAACAGTGC) respectively for RD239, RD147c and newly identified deletion of 500 bp between 3501224-3501724 (M. tuberculosis H37Rv (NC_000962.3 genome numbering). The reactions were performed with initial denaturation at 94 °C and 30 cycles of denaturation for 1 min, annealing of primers at 60 °C for 1 min and extension with platinum Taq DNA polymerase for 1 min at 68 °C. Final extensions were performed at 68 °C for 10 min. The reactions were performed as recommended by the manufacturer of the DNA polymerase.

Drug susceptibility testing

The susceptibility of the isolates to first line drugs and other second-line anti-tuberculosis drugs was investigated using the standard agar proportion method (Larsen et al., 2007). The drug concentrations used in the test comprise 0.2 mg/l isoniazid, 1.0 mg/l rifampicin, 2.0 mg/l streptomycin, 5.0 mg/l ethambutol, 1.0 mg/l linezolid, 6.0 mg/l amikacin, 5.0 mg/l ethionamide, 2.0 mg/l paraaminosalycic acid, 2.0 mg/ml ofloxacin, 2.0 mg/l moxifloxacin, 2.0 mg/l gatifloxacin, 1.0 mg/ml sitafloxacin, 6.0 mg/l kanamycin, 2.0 mg/l ciprofloxacin, 2.0 mg/l levofloxacin, and 3.0 mg/l clarithromycin. Growth equal to or more than 1% on drug containing media compared to drug free media was recorded as drug resistance. The phenotypic drug testing was performed on the initial isolates from the patients and repeated on the stock cultures.

Ethical approval

The study was approved by the Institutional review board (IRB) of Faculty of Medicine Siriraj Hospital, Mahidol University SiEC No. 152/2549.

Genome assembly

The sequences of the isolates were assembled and annotated as described in Methods. 159 contigs with N₅₀ of 69,028, 173 contigs with N₅₀ of 63,852, and 177 contigs with N₅₀ of 63,019 contigs were obtained for CSF3053, 43-5069 and 46-13838 respectively. All isolates have 65.5 % guanine/cytosine (GC) content, typical of mycobacteria. The draft genomes have an average size of 4,364,461 bp, 4,154 genes, 48 rRNAs and 64 pseudogenes. Details of the assembly and annotation are shown in Table 1.

Single nucleotide polymorphisms

Point allelic variations at any position within the genome when compared with the reference H37Rv genome (NC_000962.3) were investigated.

In total, 2,202 positions were found to have similar allelic changes (SNPs) in all isolates as shown in Fig. 1. 1,963 are in coding regions (754 synonymous, 1209 (61.6%) non synonymous) and 239 are intergenic. In this study, CSF3053, 46-5069 and 43-13838 have 10, 7 and 49 unique SNPs respectively. 43-13838 and CSF3053 have 23 SNPs in common, CSF3053 and 46-5069 have 99 SNPs in common, while 43-13838 and 46-5069 have 7 SNPs in common. Using the SNPs, the isolates were found to belong to lineage 1 with the presence of allele C/A and G/C at positions 2154724 and 7585 resulting in katG R463L and gyrA S95T respectively (Sreevatsan et al., 1997). Using a recently developed SNP barcode (Coll et al., 2014a), the isolates were found to be specific to Indo Oceanic lineage 1.2.1, with nucleotide changes G/A at position 615938, C/A at position 3479545, G/C at position 4244420 and G/C at position 9260.

The 2,202 SNPs that were found to be common to the isolates in this study were compared with 92,000 SNPs from 1,601 genomes of MTBC that were analyzed previously. These include 121, 390, 189, 856, 17, 11, and 6 genomes from lineages 1, 2, 3, 4, 5, 6, and 7, respectively. Eleven samples from M. bovis were also included (Coll et al., 2014a). The common SNPs were used to position the strains on a phylogenetic tree compared to other strains and lineages of MTBC as shown in Fig. S1. Nucleotide change G/C at position 2342203 was found only in the isolates in this study when compared with the 1,601 MTBC genomes. There is evidence from macrophage systems that strain-to-strain variability affects phenotypic outcomes (McEvoy et al., 2012). Phylogeographic strain variation may therefore have considerable effect on the development of new diagnostic tools, vaccines and drugs.

SNP C/T at position 3378828 was reported to be unique to members of lineage 1 (Coll et al., 2014a). Although this SNP was found in many genomes belonging to lineage 1, we found out that it was absent in the three isolates in this study and in 6 other Nonthaburi isolates from Thailand and the Netherlands used in previous studies which are grouped under lineage 1. This indicates that the allele change at this position may be specific only to a sub-branch of lineage 1. Synonymous SNP T/C at position 28910, non-synonymous SNP C/T at position 152178 resulting in Thr344Ile in pepA gene and intergenic SNP C/T at position 1180580 were found only in Nonthaburi isolates from Thailand. They were not found in any genome belonging to lineages 2, 3, 4, 5, 6 and 7. Within lineage 1, these SNPs were found only in Thai Nonthaburi isolates, from previous study (Coll et al., 2014a), and the isolates in this study. They were however absent in the Nonthaburi genotype isolates from the Netherlands. pepA gene is a probable serine protease with the exact function unknown. It is in the intermediary metabolism and respiration functional category. Its mRNA was found to be upregulated after 96 h of starvation (Betts et al., 2002), suggesting its role in the adaptation of mycobacteria to extreme conditions. The association of the SNPs at these positions with Thailand warrants further investigation.

Large sequence polymorphism

Region of difference RD239 that is specific to lineage 1 of MTBC and previously reported RD147c, not specific to lineage 1, were found in all the three isolates. In addition, a region of deletion of 500 bp between 3501224-3501724 (M. tuberculosis H37Rv (NC_000962.3 genome numbering) comprising Rv3135 (ppe50), was observed in all isolates. The details of the deletions as well as the affected open reading frames are shown in Table 2. The deletions were confirmed with PCR (see Figs. S2, S3 and S4). The PE-PPE protein class, while not well characterized, represents the third most abundant category of mycobacterial proteins and showed the most consistent expression during infection (Kruh et al., 2010). Although PPE50 has a yet unknown function, it was listed among promising therapeutic target in tuberculosis treatment based on its expression, and homology to human and other microbial proteins (Raman, Yeturu & Chandra, 2008). The deletion of this gene may be a means of evading recognition by the host immune system.

Deletions have been shown to have a wide range of effects on M. tuberculosis including association with an increased probability of transmission (Tsolaki et al., 2004).

Polymorphisms in drug resistance associated genes

Despite being isolated from patients with severe form of tuberculosis, drug susceptibility tests results show that the three isolates are susceptible to first line drugs; isoniazid, rifampicin, ethambutol and streptomycin, and to quinolones: ciprofloxacin, ofloxacin, gatifloxacin, moxifloxacin, levofloxacin, and sitafloxacin. They were also found to be susceptible to linezolid, amikacin, ethionamide, paraaminosalicylic acid, kanamycin and clarithromycin.

However, 37 SNPs were found in drug-resistant related genes reported in TBdream database and other earlier published reports (Sandgren et al., 2009). Nineteen are synonymous while 18 are non synonymous. Non synonymous mutations Gly312Ser of kasA gene and Ile73Thr in efpA were previously reported to be associated with isoniazid resistance (Mdluli et al., 1998; Ramaswamy et al., 2003), but were found in our isolates. Association between these mutations and resistance to isoniazid needs to be confirmed. iniA gene and Rv1592c were reported to be associated with tolerance to isoniazid (Colangeli et al., 2005; Ramaswamy et al., 2003). In our analysis, mutations His481Gln in iniA gene and Ile322Val in Rv1592c were found. These positions may not be associated with the supposed roles of these genes in isoniazid resistance.

Polymorphism exists at position 237 of nudC in M. tuberculosis isolates (Wang et al., 2011). In particular, the amino acid change Gln237Pro in nudC is found in the Indo Oceanic and West African lineages. It was demonstrated to prevent dimer formation and results in the loss of activity of the enzyme. It was also shown to degrade the active forms of isoniazid and ethionamide (Wang et al., 2011). We however found this codon change in all isolates in this study. This suggests the non-involvement of the amino acid change at this position in resistance to both drugs.

Mutations Cys110Tyr in embR, Thr270Ile and Asn394Asp in embC, Pro913Ser in embA and Glu378Ala in embB, were previously reported to be involved in ethambutol resistance (Ramaswamy et al., 2000; Srivastava et al., 2009). However, these mutations were found in this study. Mutation Ser257Pro in rmlD was suspected to be involved in isoniazid and ethambutol resistance (Ramaswamy et al., 2000). This was however found in all isolates considered in this study. Mutations Glu21Gln in gyrA, Ile322Val in Rv1592c, Arg463Leu in katG, and Arg93Leu in cycA were found to be common to the isolates in this study. They have also been reported to be common to pan-susceptible and drug-resistant M. tuberculosis sequence type 10 Beijing isolates (Regmi et al., 2015). Our results confirm that these mutations are polymorphic rather than being involved in drug resistance. The details of the synonymous and non-synonymous SNPs found in drug-resistant related genes and the predicted protein variation effects are shown in Table 3.

Polymorphisms in virulence genes, efflux pump related genes, and essential genes

Oftentimes, mutations provide selective advantage to an organism in a particular environment. Some non-synonymous mutations in rpoC gene have been shown to result in higher competitiveness in vitro and have higher fitness in vivo evidenced by their prevalence across patient populations (Comas et al., 2012). In this study, we found Ala172Val mutation in rpoC gene in all isolates.

We also sought to determine polymorphisms in genes that play important roles in the survival and pathogenesis of M. tuberculosis. Of particular interest are the genes that are involved in the evasion of the host immune system. SNPs in 37 mycobacteria virulence related genes were found to be common to the isolates. Twenty nine of the SNPs are non-synonymous. Polyketide synthases (PKs) are group of genes involved in the synthesis of polyketides which are structurally complex compounds produced by organisms for survival advantage. Some mycobacteria PKs genes such as pks15, pks1, pks10, pks12, pks5, and pks7 are known to be involved in virulence (Reed et al., 2004; Rousseau et al., 2003; Sirakova et al., 2003; Tsenova et al., 2005). Insertion of 7 base pairs was found in pks15/1 junction in all isolates. The presence of the 7 base pair insertion leads to a frame shift that results in the loss of stop codon of pks15. This results in a continuous transcription of pks15 and pks1. This was previously associated with the more virulent phenotype of the modern Beijing family, but such claim has since been refuted as it can be found across the seven lineages. The implication of the insertion needs further experiments to understand. Two mutations Ile474Met and Thr604Ala were found in nuoG gene. nuoG is a probable NADH dehydrogenase, reported to be involved in apoptosis inhibition (Velmurugan et al., 2007). Mutation Arg463Leu was found in katG, a gene previously implicated in inhibiting antimicrobial effectors of the macrophage (Ng et al., 2004). Protein kinases such as pknD and pknG are important virulent factors of M. tuberculosis. pknD has been reported to play a role in the infection of the host’s central nervous system by M. tuberculosis (Be, Bishai & Jain, 2012; Cowley et al., 2004). Gln472Pro mutation in pknD was found in all isolates. virS is a transcription regulator that belongs to AraC family. Its attenuation in a mouse model resulted in an increased animal survival (Gupta, Jain & Tyagi, 1999; Singh et al., 2003). We found mutation Leu316Arg in this gene in all isolates.

Stop codon was gained after Arg305 in PStA1, an inorganic phosphate ABC transporter. Stop codon was however lost in Rv1504. The stop codon was replaced with glutamine as codon 200. Rv1504 and PstA1 were reported to be involved in the adaptation and survival of mycobacteria in macrophages (Brodin et al., 2010; Rengarajan, Bloom & Rubin, 2005).

Non-synonymous and synonymous SNPs were found in other genes involved in various other functions related to virulence such as synthesis of complex and simple lipids, cell wall proteins, lipoproteins, cholesterol metabolism, secretion systems, protein kinases, metal transporter proteins, two component systems and other proteins of unknown functions (Table S1).

Efflux pumps play roles in drug resistance, cell physiology, detoxification and virulence of M. tuberculosis (Nikaido, 2009). Ten synonymous SNPs and 15 non-synonymous SNPs were found in efflux pump related genes. One stop codon was gained by Rv2994, a predicted transmembrane protein involved in efflux system (Table S2).

Twenty eight SNPs were observed in genes with known epitopes, 11 are synonymous while 17 are non-synonymous (Table S3).

In addition, 316 SNPs were found in essential genes, 135 are synonymous, 181 are non-synonymous. A start codon was lost in pabB gene. pabB is a cell membrane associated gene that encodes para-aminobenzoate synthetase component-I involved in the biosynthesis of p-aminobenzoate, a precursor of folate biosynthesis (Sassetti, Boyd & Rubin, 2003; Zheng et al., 2008). The details of the position, nucleotide change, amino acid change and the genes involved are presented in Table S4.

The association of the SNPs or deletions reported in this study to TBM needs further investigations. This can be done by comparing them with variations from PTB cases, to determine exclusive associations with TBM. Furthermore, The involvement of the reported allelic changes in the functions of the various genes from which they were found can be verified by site directed mutagenesis in laboratory strains of M. tuberculosis, and subsequent animal experiments.