The Next Frontier in Tuberculosis Investigation: Automated Whole Genome Sequencing for Mycobacterium tuberculosis Analysis

1. Introduction

Mycobacterium tuberculosis (MTB) is a preventable and usually curable disease, and yet it causes tens of millions of new cases of tuberculosis (TB) resulting in more than one (1) million deaths every year worldwide [1]. Multidrug-resistant (MDR) tuberculosis is a major health problem and seriously threatens worldwide TB control and prevention initiatives. Since 2020, there have been an estimated 465,000 Rifampicin-resistant (RR) tuberculosis cases, with 78% of these cases classified as multidrug-resistant (MDR). This indicates that the MTB strains causing these cases are resistant to both Rifampicin and Isoniazid. [2]. There are two main challenges with MDR TB: underdiagnosis of clinical cases and poor treatment outcomes, both of which lead to more serious clinical complications. Collectively, the failure to identify and successfully treat these cases leads to further exacerbation and transmission of MDR TB strains.

The WHO has emphasized that accurate diagnosis and characterization of TB is the key strategy for successful treatment and eradication of complicated TB cases [3]. Rapid genotypic-based diagnostic tools, such as GeneXpert, have been widely adopted as they are faster and cheaper than traditional culture-based diagnostic susceptibility testing. However, outbreaks caused by MDR strains with novel antimicrobial resistance (AMR) profiles may not be picked up by these assays, which focus on a limited number of targets. This underscores the importance of employing assays with a genome-wide approach that cover a broader range of resistance determinants [4]. Although molecular targeted assays have shown great value in TB AMR characterization [2], these methods are still dependent on known genotypic and mutation profiles [1]. Failure to detect novel AMR genotypes and mutations within the MTB genome can lead to inaccurate prescription of patient treatment regimes, posing a risk of treatment failure. [3]. Whole genome sequencing (WGS) is an agnostic approach to screen the entire pathogen genome for the detection of specific single nucleotide polymorphisms (SNPs) to discriminate known MTB complex (MTBC) (sub)lineages. This approach not only enables presence/absence detection of AMR genes but also the agnostic screening and deep characterization of AMR loci. WGS has significant potential to replace traditional culture-based phenotypic diagnostic tests, which can typically take weeks to months [5]. By reducing turnaround time to less than two (2) weeks, it will improve antibiotic stewardship for first-line drugs in TB treatment.

Widespread adoption and implementation of this methodology remains challenging due to a lack of necessary infrastructure, the high cost of entry and the requirement for technical expertise to process the samples and analyze the generated sequencing data [4]. Integrating next generation (NGS) technologies with automation will alleviate these challenges and has been proven successful for application in the food safety and infectious diseases space [6,7,8,9]. Automation imparts multiple benefits to existing manual workflows. It can improve results consistency and reproducibility by reducing the number of human touch points, eliminating human/random errors in the process. Researcher efficiency can also be improved via the removal of the need for manual and routine tasks. Automation enables the performance of workflows round the clock without human supervision and intervention, further reducing turnaround time from sample input to results output. And finally, it aids in improving workplace safety by eliminating the exposure of laboratorians to potentially harmful chemicals [10]. In this study, we evaluated the accuracy and specificity of a fully automated bacteria WGS assay for the characterization of MDR TB isolates from reference collections as well as clinical isolates.

2. Results

2.1. WGS Characterization of MTB and NTM Isolates

A total of 86 samples (44 unique isolates) were sequenced in this study using the fully automated bacteria WGS assay. MTBs and NTMs isolates were assembled, with their taxon correctly identified down to the species level and their sub-lineages characterized (Table 1). AMR genotypic profiles and associated drug resistance predictions were determined and summarized in Table 2.

2.2. Robustness and Reproducibility of the Fully Automated Bacteria WGS Assay

Three (3) strains of clinical MTB isolates provided by NSPHL, together with four (4) strains each of MTB and NTM isolates procured from reference repositories (ATCC and BEI Resources), were used to assess the reproducibility and robustness of the fully automated bacteria WGS assay. Replicates were sequenced on three (3) different systems by two (2) different operators and results are summarized in Table 3. Lyve-SET analysis of NSPHL Strains 8, 12 and 17 reported no SNP distance between replicates (see Table S1).

2.3. Study Site-Specific Clinical Isolate Relatedness Analysis

Whole genome phylogenetic analysis of clinical isolates was performed to predict MTB clusters from NSPHL and SFPHL (Figure 2a and Figure 2b, respectively). Phylogenetic groupings for the genomes within each study site were consistent with pairwise SNP distances as computed by Lyve-SET for the same genomes (Tables S1 and S2). Isolates within the same cluster shared the same node and showed lower SNP distances, while isolates from different clusters had higher SNP distances. For example, NSPHL Strains 10 and 17 that were within a single cluster had a SNP distance of 229, while NSPHL Strains 10 and 13 that were between clusters had a SNP distance of 1,979 (see Table S1). SFPHL clinical isolates 36359772, 36360344 and Z008267 that were within the same cluster had zero (0) SNP distance, while isolates Z008272 and Z008273 that were from different clusters had a SNP distance of 2,013 (see Table S2).

2.4. Comparison of Genotypic and Phenotypic AMR Results

Genotypic AMR profiles from the WGS data of the clinical isolates were compared with phenotypic data from AST provided by the study sites. The predicted AMR drug classes, based on the AMR gene profiles, perfectly matched the TB-specific AST phenotypic results from NSPHL clinical isolates (Table 4). The vast majority of clinical MTB isolates from SFPHL were sensitive to Streptomycin, Isoniazid, Rifampicin, Ethambutol, and Pyrazinamide, aligning perfectly with the genotypic AMR profiles. Three (3) isolates exhibited non-concordance between their genotypic AMR profiles and AST results. Isoniazid resistance was predicted for 36360377 while Pyrazinamide resistance was predicted for 36360355 and Z008272. AST results indicated that all three (3) isolates were sensitive to all TB first-line drugs (data not shown).

3. Discussion

Although bacteriological confirmation of TB is still the gold standard and necessary for the testing of resistance to anti-TB drugs, genetic sequencing at reference-level laboratories is slowly gaining popularity [11]. WGS is, without a doubt, a more powerful tool and can play a major role in the diagnosis of MTB drug resistance when compared with PCR-based methods. It has the advantage of enabling the identification of “off targets” or new candidate resistance mutations when facing AST discrepancies or resistance to new and/or repurposed drugs, for which resistance catalogs are still being developed [12]. In recent years, application of WGS for TB diagnosis has rapidly progressed from academic research only perspective to routine patient care in the clinics, population-level surveillance, and formulation of public health intervention strategies [13,14,15,16]. These strategies in turn help with more effective antibiotic stewardship, as well as the implementation of pathogen-based precision medicine treatments for TB [17].

Integrating WGS with automation will help reduce the complexity and hands-on time required for the experimental workflow, further reducing turnaround time for diagnosis to increase actionability. This would open up the potential of providing a turnkey solution for widespread adoption, as well as democratizing WGS technology for TB diagnosis to laboratories, hospitals, and clinics with limited personnel and technical expertise. Here, we showed that a fully automated bacteria WGS assay, using extracted DNA from isolates as starting material, is capable of generating highly robust and reproducible results with a turnaround time of approximately 24.5 hours as compared to the 6 - 10 days needed for manual workflow. The fully automated assay was able to achieve the short turnaround time by enabling around-the-clock operation without the need of human supervision and intervention, which is not the case for manual workflows restricted by stipulated laboratorian operating hours. Sample registrations, and the loading of samples and consumables took approximately 30 minutes in total. Data from sequencing the entire MTB genome provided high resolution information, conferring unique insights into AMR genotypic profiles of individual isolates, as well as how they were related to one another.

Greater depth of information helped tease out the intricate relationships of isolates classified within the same sub-lineage. For example, in the case of SFPHL clinical isolates, although isolates 36360342, 36360361, 36360364, Z008274, Z008277 and Z008278 are all from sub-lineage 1.1.1.1, isolates 36360342, 36360361 and 36360364 sub-clustered together and had a SNP distance of zero (0) (see Figure 2b and Table S2). Therefore, we hypothesized that these three (3) isolates might be the same strain or from the same origin. Isolates 36360347, 36360388 and Z008270 had a SNP distance of five (5) when compared to one another while isolate Z008279 had a SNP distance of between 173 and 194 when compared against these three (3) isolates (see Figure 2b and Table S2). We therefore hypothesized that isolates 36360347, 36360388 and Z008270 might be derived from the same outbreak while isolate Z008279 is from a separate event, even though they share the same AMR profile and 4.6.2.2 sub-lineage classification.

Correlations between genotypic and phenotypic profiles of AMR have been widely debated, with varying concordances reported dependent on drug classes and organisms tested [18,19,20,21]. Similarly high, but imperfect concordance has been reported for MTB [12,22]. Findings in this study are consistent with those in the literature. Clinical isolates from NSPHL revealed 100% concordance between WGS genotypic AMR profiles and AST phenotypes, further supporting the use case of applying WGS for comprehensive TB diagnosis under clinical settings. Out of the 24 clinical isolates from SFPHL, 21 showed agreement between genotypic and phenotypic AMR profiles, while three (3) isolates had some discrepancies. These three (3) isolates had AMR genotypes predicting resistance against Isoniazid or Pyrazinamide, despite AST data indicating susceptibility to both drugs in all three cases. The authors emphasize that achieving consistent data without human error is a significant achievement, far outweighing the few discrepancies. This outcome strongly endorses the use of a fully automated WGS assay in TB diagnosis. While the results from the genotypic AMR profile may impact treatment options, in this case, the “false positive” results would encourage clinicians to proceed with caution in treatment selection. This cautious approach further promotes antibiotic stewardship efforts in choosing and tailoring appropriate first-line drug treatments for specific TB cases.

4. Conclusion

In conclusion, the present study highlighted the suitability and robustness of a fully automated bacteria WGS assay for seamless, accurate detection and in-depth characterization of MTB clinical isolates and NTM reference isolates with high consistency. The assay demonstrated 100% concordance in species and sub-lineage classification, and 92% (35 out of 38 isolates) concordance in genotypic vs phenotypic AMR characterization. The WGS data generated will further propel research in the investigation of different MTB virulence factor mechanisms, such as secretion factors, cell surface components, enzymes involved in cellular metabolism, and transcriptional regulators [23]. This fully automated WGS solution, with a short turnaround time, can further aid in the surveillance and resolution of TB outbreaks in a timely manner. A significantly shorter (less than 2-day) turnaround time will increase actionability and patient safety by reducing and preventing MDR TB transmission and related fatalities, as well as reducing healthcare costs [24,25]. A fully automated bacteria WGS assay will be a valuable technology and tool to enhance the efficiency of diagnostic laboratories while supporting WHO’s vision in the elimination of this deadly disease. This study has implications for both government and nongovernmental organization resource allocation decisions and policy formulation, thereby establishing a favorable proposition in the adoption of a fully automated turnkey WGS solution for the decentralization of precision medicine and TB diagnostics.

5. Materials and Methods

5.1. Sample Source

Bacterial isolates or extracted DNA from reference MTB and non-tuberculosis Mycobacterium (NTM) were procured from ATCC or BEI Resources, NIAID, NIH. DNA extracted from de-identified clinical isolates were obtained from two (2) study sites: the Nevada State Public Health Laboratory (NSPHL) and the San Francisco Public Health Laboratory (SFPHL). Details of isolates and strains used in this study are summarized in Table 5.

5.2. Isolate Cultures and DNA Extraction

Non-tuberculosis Mycobacterium (NTM) and an attenuated MTB strain (NR-122) extracted at Clear Labs, Inc. were first grown on Löwenstein–Jensen (LJ) Agar slants (Hardy Diagnostics, Santa Maria, CA) incubated at 35°C ± 2°C for 7-10 days. Bacterial isolates were resuspended in Clear Labs’ Resuspension Buffer v2.0 (Clear Labs, Inc., San Carlos, CA) and heat inactivated at 95°C for 30 minutes. MTB cell suspensions were lysed and DNA extracted using a proprietary Clear Labs inhouse-developed extraction protocol. Extracted DNA was stored at -80°C until use.

MTB clinical samples provided by NSPHL were grown on LJ agar slants, resuspended in HPLC water and inactivated by boiling for 30 minutes. DNA was then extracted from heat inactivated cultures using the Maxwell RSC Cultured Cell Kit (Promega, Madison, WI) or DNeasy Blood & Tissue QIAcube Kit (QIAGEN, Venlo, Netherlands). Extracted DNA was stored at -80°C until use.

MTB clinical samples provided by SFPHL were grown on LJ agar slants, resuspended in HPLC water and inactivated by boiling for one (1) hour. DNA extractions were then performed using the MagNA Pure 24 Total NA Isolation Kit (Roche Life Sciences, Indianapolis, IN) or the EZ1&2 Virus Mini Kit v2.0 (Qiagen, Germantown, MD) as detailed in Table 1. Extracted DNA was stored at -80°C until use.

5.3. DNA Quantification

Extracted DNA was quantified on a Qubit™ 4 Fluorometer using the Qubit 1x dsDNA High Sensitivity Kit according to the manufacturer’s instructions (Thermo Fisher Scientific, Waltham, MA).

5.4. Workflow of a Fully Automated Platform for WGS of Extracted DNA from MTB and NTM

Extracted DNA from all isolates were first diluted to a concentration of 0.5 ng/µL in 30 µl volume (total DNA input yield of 15 ng). Samples were prepared, sequenced and analyzed using the “Low DNA Input” protocol of the Clear Dx™ Microbial Surveillance WGS v2.0 application on the fully automated and integrated Clear Dx™ platform, which comprised liquid handling robotics, thermal cycler and sequencers (Figure 1). Briefly, samples, reagents, and consumables were first loaded onto the instrument according to the manufacturer’s instructions (Clear Labs, San Carlos, CA). Up to 12 samples per run were then registered using the ClearView App with “Mycobacterium (4.4 MB)” selected as “Organism” and “30 – 40x” selected as “Coverage”. Upon pushing the “Start Run” button, the fully automated system performs library preparation using the 2 x 150 bp chemistry, followed by the loading of sequencing cartridges onto the Illumina iSeq100 sequencers (Illumina, Inc., San Diego, CA) as a part of the end-to-end workflow. After completion of sequencing, raw sequencing data was automatically uploaded onto the cloud for bioinformatics analysis. Final output generated by the workflow comprised compressed FASTQ files, sample identity (genus and/or species level), and quality metrics for each isolate. The fully automate end-to-end workflow took approximately 24.5 hours in total.

5.5. Data Analysis

The TheiaProk bioinformatics pipeline was used for bioinformatics analysis [26]. Briefly, paired-end Illumina reads in FASTQ format generated by the fully automated bacteria WGS assay were first trimmed and quality filtered using Trimmomatic and fastq_scan pipelines. Cleaned reads were then de novo assembled using the Shovill algorithm, assembly quality control (QC) checked by QUAST and BUSCO pipelines, and taxonomic assignment performed using GAMBIT. Further downstream analysis includes genome annotation with Prokka, AMR characterization with AMRFinderPlus, and TB-specific analysis using TBProfiler. MashTree was used for site-specific phylogenetic analysis of isolates while Lyve-SET [27] was used for SNP analysis with M. tuberculosis H37Rv (NC_000962.3) set as the reference genome. Antimicrobial susceptibility test (AST) results of MTB isolates were provided by NSPHL and SFPHL.