Tuberculosis (TB) caused by Mycobacterium tuberculosis (MTB) has been a global public health challenge for decades. It was estimated that there were 10.4 million new cases and 1.8 million TB-related deaths worldwide annually (World Health Organization, 2016). Despite the enormous global efforts in TB control, multidrug resistant TB (MDR-TB) with resistance to at least rifampicin (RIF) and isoniazid (INH) is starting to threaten the treatment regimens currently available. In 2016, World Health Organization estimated a total of 480,000 TB cases caused by MDR-TB. Among these MDR-TB cases, ~9.5% were extensively drug resistant TB (XDR-TB) with additional resistance toward fluoroquinolones and one of the other second-line injectable anti-TB drugs.

Understanding the molecular basis of drug resistance and the transmission of MDR-TB is critical in developing new therapeutic and preventive strategies against the disease. However, while over 95% of phenotypic RIF resistance can be explained by the mutations at 81 bp rifampicin resistance determining region (RRDR) of rpoB (McCammon et al., 2005), katG S315T missense mutation and mutations at the promoter region of inhA can only explain 50–95 and 8–43% of phenotypic INH resistance respectively (Witney et al., 2016). The discordance between phenotypic and genotypic drug resistance profile suggested that genetic variations among MDR-TB strains are still far from fully elucidated. Whole genome sequencing using short read sequencing platform allows identification of additional resistance-related mutations, but there are still limitations to the current method. First of all, the repetitive nature and high GC content of MTB genome often result in amplification bias during library preparation, which in turn causes fragmented genome assembly (Pop and Salzberg, 2008; Treangen and Salzberg, 2011). Secondly, although there are numerous pipelines for accurate Single Nucleotide Polymorphisms (SNPs) identification, bioinformatics tools in determining other genetic variations such as indels and copy number variations are still very much limited (Ratan et al., 2015; Steglich and Nubel, 2017). More importantly, comparative genomic analyses between drug-sensitive TB and MDR-TB strains often result in thousands of potential mutations. With the vast majority of these mutations being lineage-specific SNPs or neutral mutations accumulated during bacterial replication, identification of novel resistant-related mutations from other irrelevant gene candidates would require massive amount of computation analysis and experimental validations (Fan et al., 2014).

Single Molecule Real Time (SMRT) Sequencing is a recently developed sequencing platform that has greatly simplified the de novo assembly step of small microbial genomes (Roberts et al., 2013). With an average of 10–20 kbp per read, SMRT sequencing can easily span across highly repetitive DNA sequences (Scott and Ely, 2015; Zhu et al., 2016). Moreover, library preparation for SMRT Sequencing does not require PCR amplification, which means there would be virtually no GC bias during sequencing step (Ferrarini et al., 2013). In other words, SMRT Sequencing is able to reduce the number of gaps in the final assembled genomes while accurately determine the presence of any complex structural variations. The unbiased GC coverage and the extraordinarily long reads in SMRT Sequencing makes the platform particularly suitable for sequencing MTB genomes with unexceptionally high GC content and large numbers of repetitive regions (Ferrarini et al., 2013).

In this study, we focused on two consecutive phenotypically different MDR-TB clinical isolates collected from the same patient during the course of anti-TB regimen. Using SMRT Sequencing platform, genomic libraries prepared from the two strains were sequenced, de novo assembled and analyzed. The clonal relatedness between the two MDR-TB clinical isolates allowed an extensive removal of lineage-specific SNPs by comparative genomic analysis. Mutations responsible for the phenotypic differences could thus be revealed.

This is the first study describing the stepwise evolution of drug resistance and growth fitness in MDR-TB during in vivo challenge of anti-TB drugs. In this study, we focused on a pair of MDR-TB clinical isolates obtained from the same patient within a 4 months interval of anti-TB therapy. The post-treatment strain exhibited improved in vitro growth rate and elevated resistance toward INH, RFB, and EMB when compared to the pre-treatment strain. The clonal relationship between the two MDR-TB clinical isolates was confirmed by identical MIRU-VNTR and spoligotyping patterns, which indicated that the post-treatment strain was developed from pre-treatment strain during the course of anti-TB treatment. Since these two strains were clonally related, comparative genomic analysis could potentially filter out all the lineage-specific mutations, and genetic alterations responsible for the phenotypic differences could be revealed.

The application of SMRT Sequencing in microbial genome assembly has been evaluated extensively in recent years, and it is considered as a robust sequencing platform for MTB genomes with high GC content (Koren et al., 2013; Miyoshi-Akiyama et al., 2015; Rodríguez et al., 2015; Valafar et al., 2015). Despite its high single-pass error rate (15.4–18.7%), previous studies have proved the high consensus accuracy of SMRT Sequencing (Quality Score ≥ Q50) given a sufficient coverage depth (>20X) and the use of error-correction step prior to assembly (Koren et al., 2012). In this study, SMRT Sequencing on both MDR-TB clinical isolates achieved fully closed de novo assembled genomes represented by one single contig with uniform coverage. Since the MDR-TB strains were clonally related, it was expected that there would only be a limited number of genetic variations across the two strains. Consistent with our hypothesis, comparative genomic analysis between the two strains only identified 15 polymorphisms throughout the entire bacterial genomes. The result has further demonstrated the consensus accuracy of SMRT Sequencing platform in analyzing the genomic sequences of MTB.

To explain the antibiotic resistance patterns in pre-treatment strain and post-treatment strain, DNA sequences of resistance-related genes were examined from the assembled genomes and the identified mutations were confirmed by Sanger Sequencing. The genes were chosen based on their correlations with antibiotic resistance. Both strains showed concordant genotypic and phenotypic susceptibility profile toward RIF, PZA, STR, OFX, MOX, KAN, AMI, ETH, CYS, and PAS. Both MDR-TB strains were resistant to PZA at 100 mg/L while carrying pncA missense mutation G162D. The identical resistance profile of RIF (MIC > 16 mg/L), STR (MIC = 32 mg/L), KAN (MIC > 40 mg/L), and AMI (MIC = 4 mg/L) in both strains could be explained by the missense mutation at rpoB S531L, rpsL R43K, and rrs A1401G respectively (Tracevska et al., 2004; Chan et al., 2007; Jugheli et al., 2009; Jamieson et al., 2014). Meanwhile, the wild-type folC, gyrA, gyrB, and ethA sequences of both MDR-TB clinical isolates were concordant to the phenotypic susceptibility toward PAS, OFX, and MOX and ETH (Zheng et al., 2013; Zhang and Yew, 2015).

Interestingly, while both strains carried identical rpoB S531L missense mutation, the RFB MIC of the post-treatment strain was increased by 16-folds (MIC = 16 mg/L) without any additional mutations at the entire rpoB sequence. Promoter mutations at mabA-inhA and missense mutations at katG codon S315 were known to be the major cause for low level (MIC = 0.2–0.8 mg/L) and high level (MIC > 1 mg/L) INH resistance, respectively (Abe et al., 2008; Ando et al., 2010). Although both MDR-TB clinical isolates only carried inhA C-15T mutation, the post-treatment strain demonstrated high level of INH resistance (MIC = 4 mg/L) in the absence of katG missense mutations. Similarly, both strains carried embB M306V missense mutation in the absence of any missense ubiA mutations. The embB M306V is predominantly the cause for moderate EMB resistance (MIC = 6–14 mg/L) (Safi et al., 2008; Starks et al., 2009), and can be elevated by 2-folds in the presence of ubiA missense mutations (MIC = 16–32 mg/L) (Safi et al., 2013; He et al., 2015). While the pre-treatment strain demonstrated a susceptible to borderline resistance profile toward EMB (MIC = 8 mg/L), the post-treatment strain was resistant to high level EMB resistance with 4-fold elevation in EMB MIC (MIC = 32 mg/L). The discrepancies in the genotypic and phenotypic resistance patterns have hinted the possibility of novel antibiotic resistance development pathways. It is also notable that the post-treatment strain demonstrated an improved in vitro growth rate when compared to the pre-treatment strain, which was contrary to general belief that drug resistance development in bacteria is often associated with a fitness deficit. Multiple studies have demonstrated that bacterial fitness could be potentially restored by the presence of compensatory mutations (Comas et al., 2011; He et al., 2015). Our result has further suggested the acquisition of compensatory mutations in the post-treatment strain that restored its growth fitness during the development of additional antibiotic resistance.

Further comparative genomic analysis identified three additional polymorphisms between the two MDR-TB clinical isolates. To investigate the importance of these polymorphisms, multiple protein sequence alignment was performed using orthologs from diverse pathogenic mycobacteria strains representing different lineages of their respective species, which covered most of the pathogenic mycobacterial species commonly identified in Hong Kong. We hypothesized that residue essential for protein functions should be residing in highly conserved regions across different species. For residues that were non-essential to the protein function, polymorphisms would occur even among different lineages of the same mycobacterial species. In other words, if the polymorphisms identified in SMRT Sequencing were located within conserved regions, the mutation is very likely to have significant impact on the resultant protein function.

The first polymorphism was a single nucleotide deletion at codon 120 (nucleotide position 360) in rv0888 for pre-treatment strain, resulting in the truncation of Rv0888 from 490 residues to 136 residues. Rv0888 is an extracellular nuclease highly active in degrading various types of nucleic acids. The activity of extracellular nuclease could be essential for bacterial virulence, utilization of extracellular DNA as energy source and the degradation of neutrophil DNA extracellular traps (Dang et al., 2016). Moreover, Rv0888 could utilize sphingomyelin from host cell as an additional carbon and energy source during intracellular replication (Speer et al., 2015). Site-directed mutagenesis conducted by Dang et al. demonstrated that His353, Asp387, and Asp438 were critical to Rv0888 nuclease activity (Dang et al., 2016). Truncation of Rv0888 causes the deletion of the aforementioned functional amino acids, thus renders the protein non-functional in pre-treatment strain. The restoration of frameshift mutation back to a functional protein in post-treatment strain could potentially indicate the importance of Rv0888 during the stepwise development of MDR-TB.

The second polymorphism was lpdA V44I missense mutation in post-treatment strain. Using multiple sequence alignment result, we demonstrated that V44I was occurring at a highly conserved region in LpdA. Further structural-based analysis demonstrated that V44I was a mutation near the FAD binding site of LpdA and was able to cause protein destabilization. The lpdA gene encodes for probable lipoamide dehydrogenase (Argyrou and Blanchard, 2001; Argyrou et al., 2004). Previous reports showed that LpdA belongs to oxidoreductase family which utilizes NADP/NADPH⁺ to maintain the pool of reduced pyridine nucleotides in MTB. LpdA was also suggested to be a quinone reductase that transfers reducing equivalents from the reduced pyridine nucleotide pool to the electron transport chain and protects the bacteria from oxidative stress in in vivo conditions (Akhtar et al., 2006; Zheng et al., 2008). The effect of LpdA to drug resistance development is currently unknown. The involvement of LpdA in maintaining the redox equilibrium of MTB might suggest a potential association of LpdA with drug that requires bacterial catalase activation such as INH. Further experimental validations would be required to identify the association of lpdA V44I mutation with the elevated INH resistance in post-treatment strain.

The third polymorphism was a deleterious cobM inframe deletion at codon 69–70 in post-treatment strain. The exact functional role of CobM is still currently unknown, but it is believed that the gene encodes for Precorrin-4 C (11)-methyltransferase which catalyze the methylation of C-11 in precorrin-4 to precorrin-5 during Vitamin B12 metabolism in MTB (Young et al., 2015), meaning that the gene function could be associated with the bacterial growth rate. Interestingly, the 6-bp inframe deletion occurred at cobM codon 69–70 would removed one of the three [Ala-Asp] repeating units at codon 67–72 in post-treatment strain. The occurrence of such inframe deletion specifically in post-treatment strain strongly indicated the potential correlation between the mutation and the improvement of the survival fitness for the post-treatment strain under antibiotic selective pressure.

In this study, 11 polymorphisms identified by SMRT Sequencing were located within a 100 bp repetitive region at the center of PE-PGRS2. Due to the high GC content and the repetitive nature of the gene, we were unable to perform confirmatory Sanger sequencing to verify the mutations in PE-PGRS2. Nevertheless, the inability of confirmatory Sanger Sequencing at PE-PGRS2 in this study has further highlighted the effectiveness of SMRT Sequencing in analyzing the whole genome of MTB. Recent studies have illustrated the potential importance of PE-PGRS and PE-PPE in related to MTB pathogenicity (Kohli et al., 2012; Fishbein et al., 2015). However, the gene family of PE-PGRS were known to be highly variable with a unusually high GC content, they are prone to high frequency of sequencing errors and the effect of the mutations were often difficult to be predicted (Mcevoy et al., 2012; Golby et al., 2013; Tyler et al., 2016). When analyzing the whole genome of MTB using short read sequencing platform, gaps were often generated at regions with repetitive elements, in particularly within PE-PGRS, PE-PPE, cysA gene, and insertion elements (Gurjav et al., 2016). Unlike short read sequencing, reads generated by SMRT Sequencing were ~10–20 kbp, which can effectively covers the repetitive regions in MTB genome. The long sequencing read of SMRT Sequencing has greatly simplify the procedures for identifying mutations at PE-PGRS and PE-PPE, and thus genetic makeup of MTB could be investigated in greater details.

There are limitations in this study. Only one pair of MDR-TB clinical isolates was included for comparative analysis. The absence of additional comparison groups was due to the fact that stepwise development of MDR-TB clinical isolates with elevated drug resistance profile is extremely rare in regions adopting Direct Observed Treatment, Short Course regimen. Although stepwise elevation of drug resistance can also be achieved by inductive mutagenesis and clonal selection, these strains were derived in laboratory and might not reflect the actual situation in clinical settings (Safi et al., 2013; Smith et al., 2013). On the other hand, even when strains of the same lineage were used for comparison, it would still result in thousands of potential mutations and further extensive experimental validations would be required. The clonal relatedness of pre-treatment and post-treatment strain means that lineage-specific mutations could be extensively removed during comparative genomic analysis, and the mutations responsible for the phenotypic differences could be revealed.

In this study, growth fitness of the two MDR-TB strains were evaluated by in vitro growth curve experiment. Previous studies measured the relative fitness of TB strains using pairwise growth assays, in which the in vitro growth of bacteria were analyzed by automated liquid culture systems such as MGIT 960 System and MB/BacT ALERT 3D System (O'sullivan et al., 2010; Safi et al., 2013; Satta et al., 2016). In this study the growth rate of MDR-TB strains were continuously monitored every 48–72 h over a period of 14 days. By using the in vitro growth curve, we were able to provide a much more comprehensive evaluation on the growth fitness of MDR-TB at different time-point.

This study has demonstrated that the phenotypic resistance profiles of MDR-TB clinical isolates could be explained by the well-characterized genetic mutations at resistance-related genes. However, further experimental validations are still required to delineate the causal relationship between the newly identified mutations and MTB physiology. In addition, although we have demonstrated the different in vitro growth rate of the two MDR-TB clinical isolates, macrophage challenge experiments would be done to fully evaluate their in vivo growth fitness. The mutated genes would be transformed into both pre-treatment and post-treatment strains by electroporation. The antibiotic susceptibility profiles and the growth fitness of the transformed strains will then be examined in detail to prove the role of the mutations for the stepwise development of MDR-TB.

In conclusion, our results described the robustness of SMRT Sequencing in identifying mutations among MDR-TB clinical isolates. Comparative genomic analysis between the two clonally related MDR-TB revealed additional mutations at lpdA, cobM, and rv0888 that were developed as a result of anti-TB treatment. The underlying genetic alterations involved could explain the improved in-vitro growth fitness and the elevated drug resistance of the post-treatment strain when compared to the pre-treatment strain. Our work has suggested novel gene targets that could potentially play an important role in the physiological characteristics among MDR-TB strains.

KS-SL, GK-HS, KK-GT, SW-CT, and RR were responsible for planning and conducting the experimental works. P-LH, SS-YW provided the clinical support to this study. WZ and OC-KM provided technical support to PacBio Sequencing. W-CY is the corresponding author and supervise the project.