Whole-Genome Sequencing Exhibits Better Diagnostic Performance than Variable-Number Tandem Repeats for Identifying Mixed Infections of Mycobacterium tuberculosis

ABSTRACT Mixed infections of Mycobacterium tuberculosis, defined as the coexistence of multiple genetically distinct strains within a single host, have been associated with unfavorable treatment outcomes. Different methods have been used to detect mixed infections, but their performances have not been carefully evaluated. To compare the sensitivity of whole-genome sequencing (WGS) and variable-number tandem repeats (VNTR) typing to detect mixed infections, we prepared 10 artificial samples composed of DNA mixtures from two strains in different proportions and retrospectively collected 1,084 clinical isolates. The limit of detection (LOD) for the presence of a minor strain was 5% for both WGS and VNTR typing. The overall clinical detection rate of mixed infections was 3.7% (40/1,084) for the two methods combined, WGS identified 37/1,084 (3.4%), and VNTR typing identified 14/1,084 (1.3%), including 11 also identified by WGS. Multivariate analysis demonstrated that retreatment patients had a 2.7 times (95% confidence interval [CI], 1.2 to 6.0) higher risk of mixed infections than new cases. Collectively, WGS is a more reliable tool to identify mixed infections than VNTR typing, and mixed infections are more common in retreated patients. IMPORTANCE Mixed infections of M. tuberculosis have the potential to render treatment regimens ineffective and affect the transmission dynamics of the disease. VNTR typing, currently the most widely used method for the detection of mixed infections, detects mixed infections only by interrogating a small fraction of the M. tuberculosis genome, which necessarily limits sensitivity. With the introduction of WGS, it became possible to study the entire genome, but no quantitative comparison has yet been undertaken. Our systematic comparison of the ability of WGS and VNTR typing to detect mixed infections, using both artificial samples and clinical isolates, revealed the superior performance of WGS at a high sequencing depth (~100×) and found that mixed infections are more common in patients being retreated for tuberculosis (TB) in the populations studied. This provides valuable information for the application of WGS in the detection of mixed infections and the implications of mixed infections for tuberculosis control.

can complicate drug-susceptibility testing (DST) because the resistance present in the minority strain may not be detected (2,3). As a result, a drug regimen based on DST results that only detect the majority strain could be inappropriate and lead to treatment failure. Cohen et al. (4) used mathematical models to show that mixed infections could also compromise the expected benefits of new antituberculosis vaccines if they do not adequately cover all circulating strain types. Mixed infections can also cause M. tuberculosis transmission links to be missed if an undetected strain in a patient with a mixed infection is transmitted to another patient (5,6).
Current approaches to identify mixed infections employ molecular genotyping methods and whole-genome sequencing (WGS). The classical molecular genotyping methods include spoligotyping (7), IS6110 restriction fragment length polymorphism (RFLP) typing (8), and variable-number tandem repeat (VNTR) typing (9). Among these, VNTR typing is currently the most widely used method because of its efficiency and simplicity (1,10). However, VNTR typing examines heterogeneity within only a limited set of loci, which limits its sensitivity. WGS offers superior resolution compared to molecular typing methods (11) and is increasingly used to detect mixed M. tuberculosis infections (12,13). In this report, we compared the accuracy of WGS and VNTR for detecting mixed infections in artificial and clinical samples and identified risk factors for mixed infections. This work should provide theoretical background for further understanding mixed infections and developing effective control strategies.

RESULTS
The limit of detection (LOD) of WGS and VNTR typing for mixed infections. We first investigated the sensitivity of WGS and VNTR typing for detecting mixed infections in artificial samples containing DNA mixtures in various proportions from two M. tuberculosis strains belonging to the same (group A) or different (group B) phylogenetic lineages (Fig. 1A). Sequencing libraries were constructed using the minimum genomic DNA (50 ng, ;1.1 Â 10 7 genomes) recommended by the manufacturer. WGS analysis revealed that the average sequencing depth of these artificial samples was 134Â (105 to 171Â), with an average genomic coverage of 99.5%. The LOD for mixed infections was the presence of each strain in a proportion of at least 5% (;5.5 Â 10 5 genomes). This 5% LOD was found in mixtures from both groups A and B, and the ratio of the two different strains in the WGS data was consistent with the actual mixing ratio ( Fig. 2A). VNTR typing, using the same quantity of genomic DNA as used for WGS, similarly had a LOD for mixed infection of 5% in the artificial samples from both groups A and B ( Competition in PCR affected the LOD of VNTR typing for mixed infections. Previous studies have shown that the minimum amount of template DNA for PCR amplification of the VNTR loci is about 440 copies of the genome (14), but we found that the minor strain was undetectable when it was present in the artificial samples at ratios less than 5%, which is calculated to represent about 5.5 Â 10 5 genome copies. To investigate the minimum amount of template DNA required for VNTR to identify mixed infections, we performed 10-fold dilutions of the artificial samples. We found that when the minor strain was present in a proportion of at least 5%, it was still detected when the 50 ng total DNA in artificial mixtures A and B was diluted 10 23 and 10 22 , respectively, although the number of heterozygous loci decreased as the quantity of template DNA was reduced (Fig. 2B, Fig. S2).
PCR amplification of different DNA fragments flanked by the same primer sites can result in competition between the targets (15). We therefore speculated that competition in PCR might lead to a decreased LOD for detecting mixed infections using VNTR. To test this possibility, we determined the minimum amount of DNA required for VNTR amplification when only the minor strain was present in group A and B samples. PCR amplification of the minor strains alone was achieved with a minimum of 0.005 ng of DNA (1.1 Â 10 3 genomes), which corresponds to 0.01% (0.005 ng/50 ng) of the total DNA present in the mixed artificial samples (Fig. 2C).
Identification of mixed infections in clinical M. tuberculosis isolates. To compare the ability of WGS and VNTR typing to detect mixed infections in clinical isolates, we analyzed the isolates grown from pretreatment sputum samples of 1,084 tuberculosis (TB) patients. Phylogenetic analysis based on their WGS data (Fig. 3A Fig. 3B). The minor strain was present in a proportion of 25% or less in 94.6% (35/37) of the mixed infections (Fig. 3C), and no heteroresistance was detected in any of these 37 isolates. VNTR typing detected 14 isolates (1.3%, 14/1,084) with mixed infections: 1 isolate had 4 heterogeneous VNTR loci, 1 had 3 heterogeneous VNTR loci, and the remaining 12 isolates had 2 heterogeneous loci. Combining WGS and VNTR results, mixed infections were detected in a total of 40 isolates (3.7%, 40/1,084). WGS exhibited a higher detection rate than VNTR typing (3.4% versus 1.3%, P = 0.001), but only 78.6% (11/14) of mixed infections identified by VNTR typing were also detected by WGS (Fig. 3D).
Analysis of risk factors for mixed infections. To identify the risk factors associated with mixed infections, we compared the demographic, clinical, and bacteriological characteristics of cases with and without mixed infections ( Table 1). The proportion of retreatment cases was much higher among the 40 patients with mixed infections than in the patients infected with a single strain (20.0% versus 8.2%, P = 0.046). Beijing strains were slightly, but not significantly, more frequent in patients with mixed infections than in other patients (90.0% versus 79.0%, P = 0.081). In the multivariate logistic analysis, adjusted by Beijing genotype, retreatment patients had 2.7 times (95% confidence interval [CI], 1.2 to 6.0) the risk of mixed infections as new cases.

DISCUSSION
In this study, we compared the capacity of whole-genome sequencing and VNTR typing to detect mixed infections. Based on artificial samples, we found that both WGS and VNTR typing could detect mixed infections when the genomic DNA from the minority strain was present in a proportion of at least 5%. Combining the results from both methods analyzing 1,084 clinical isolates, the total proportion of mixed infections was 3.7% (40/1,084), with WGS detecting more mixed infections than VNTR typing (37/ 1084, 3.4% versus 14/1,084, 1.3%). The concordance between the two methods was 97.3% (1,055/1,084). Multivariable analysis showed that retreated patients were at significantly higher risk of having a mixed infection than new cases.
The main factors affecting the sensitivity of VNTR and WGS for detecting mixed infections were the genomic detection range and the completeness of the reference evolutionary tree. WGS can provide comprehensive genetic information and therefore has a higher resolution than VNTR, which is based on limited genomic regions (1, 10), and this higher resolution implies higher sensitivity to detect genomic differences (16). Many studies of molecular epidemiology have shown that the clustering rate determined by WGS is lower than that determined by VNTR because WGS is better able to discriminate between similar but distinct strains (17,18). However, because the identification of mixed infections by WGS depends on identifying two distinct mutational phylogenetic paths that are present in the database (12), if the divergent paths of the two strains are not included in the database, the mixed infection will not be detected (Fig.  S1). The three cases of mixed infections identified by VNTR but missed by WGS illustrate that the phylogenetic-based WGS method used in this study is not perfect. The sensitivity of this method for detecting mixed infectious could perhaps be improved if the reference evolutionary tree incorporated more local strains.
Recently, many other WGS-based methods have been developed for the detection of mixed infections. MixInfect (19) is a method for mixed M. tuberculosis infection estimation that uses a Bayesian model-based clustering technique. While this method can estimate mixture proportions, it does not provide any functionality for resolving the constituent strains. QuantTB, described by Anyansi et al. (13), quantifies the multiplicity and abundance of mixed M. tuberculosis infections based on a reference database. Similar to our phylogenetic method, its performance is highly dependent on the database's representation of the common strains in the relevant local context. Gabbassov et al. (20), developed a tool called SplitStrains that is grounded in a rigorous statistical framework. This tool can estimate the M. tuberculosis mixture proportions and partially reconstruct their sequences using an expectation-maximization (EM) algorithm. However, the EM algorithm has a very slow convergence, and it makes the convergence to the local optima only. In addition, several metagenomic tools exist to classify mixed populations of strains within a single species (21), but these methods were not designed to be able to discriminate between strains of highly clonal species such as M. tuberculosis.
Although we determined that the LOD for mixed infections using both WGS and VNTR  (22). The average sequencing depth in this study was about 100Â, and single nucleotide polymorphisms (SNPs) with frequencies below 5% were eliminated to exclude possible sequencing errors (23). In contrast, Gan et al. (12) used deep WGS with a depth of ;1,000Â and reported that the LOD for mixed infections was 0.64%. Although this depth (1,000Â) is currently not routine when performing WGS of many M. tuberculosis isolates, it could become feasible with improved sequencing technology, thereby reducing the LOD for detecting mixed infections. For PCR-based VNTR typing, despite the high sensitivity of PCR amplification for a single strain, the PCR target competition in mixed infections appears to be an obstacle to lowering the LOD. Previous studies have reported different prevalences of mixed infections, but they employed different methods and study designs. Our finding that WGS is more sensitive for detecting mixed infections suggests that reports based only on VNTR typing may have underestimated the prevalence of mixed infections (8). In contrast, Fang et al. (24), using VNTR typing, reported mixed infections in 5.6% of clinical isolates in Shanghai, China, higher than the 3.7% we found using both WGS and VNTR. The difference may be because their definition of a mixed infection required only one heterogenous locus, while our definition required two (25). Peng et al. (26) analyzed two sputum samples per patient and detected mixed infections in 11.2% of TB cases in Heilongjiang, China, but the proportion of mixed infections decreased to 5.6% when the analysis was based on only one sample. This is consistent with recent findings that sequential sputum samples can have different representations of in-host bacterial diversity (27).
WGS, although not perfect, was superior to VNTR for identifying mixed infections. Another advantage of WGS is that the data can be easily stored, subjected to additional analyses in the future, and compared with results from different laboratories (11). This will provide a foundation for further research not only on the detection of mixed infections but also on their significance and why they occur. For example, most previous studies have focused on mixed infections caused by strains with different resistance profiles, but more research is needed on the clinical implications of mixed infections with strains having the same resistance profile (9,28). For example, do the two strains have distinct transmission potentials that could be explained by identifying determinants in the bacterial genome (29)? Are there host or environmental factors that predispose patients to mixed infections (30,31)? These questions may be amenable to analysis with WGS. Furthermore, we detected no mixed infections involving different lineages, which may indicate that different lineage types could vary in their susceptibility to mixed infections. The simultaneous interaction between multiple strains and the host immune system in a mixed infection scenario remains to be explored.
This study had several limitations. The inconsistencies in the performance of WGS and VNTR for detecting mixed infections in laboratory versus clinical samples may be due to insufficient strain genotype combinations in the artificial samples; e.g., there were no combinations of strains that cannot be classified by VNTR typing. Mixed infections detected by WGS should be further confirmed using additional methods, such as single-colony sequencing.
In conclusion, our systematic comparison of the ability of WGS and VNTR to detect mixed infections, using both artificial samples and clinical isolates, revealed the superior performance of WGS at a high sequencing depth (;100Â) and found that mixed infections are more common in patients being retreated for TB in the populations studied.

MATERIALS AND METHODS
Preparation of artificial samples and selection of clinical M. tuberculosis isolates. Genomic DNA was extracted in the biosafety level 2 plus (BSL-21) laboratory from M. tuberculosis strains using the cetyltrimethyl-ammonium bromide-lysozyme (CTAB) method (32). Each artificial sample was composed of mixtures of DNA from two strains in different proportions (Table S1). Two groups of DNA mixtures were prepared (Fig. 1A): group A contained DNAs from two strains belonging to different lineages (L2 and L4) that were separated by a genetic distance of 1,159 single nucleotide polymorphisms (SNPs); group B contained DNAs from two strains belonging to the same lineage (L2) that had a genetic distance of 392 SNPs. DNAs from the two strains were diluted separately in double-distilled water (ddH 2 O) to a final concentration of 50 ng/mL and mixed in each of the following proportions: 25%:75%, 5%:95%, 0.5%:99.5%, 0.1%:99.9%, and 0.01%:99.99%. We also prepared 10 additional samples consisting of only the minor strain in group A (M1) and B (M2) samples, with the same concentrations of the strain DNA as in the artificial mixtures.
The Songjiang District is located in the southwest of Shanghai, China, and has an incidence of pulmonary tuberculosis (TB) of approximately 33.6/100,000. All suspected TB cases in Songjiang are referred to the local TB-designated hospitals, and all culture-positive M. tuberculosis isolates are collected and stored in the Shanghai Municipal Centre for Disease Control and Prevention (SCDC). The SCDC provided our study with all 1,084 stored isolates cultured from the sputum samples of TB patients diagnosed from 2011 to 2015 in the Songjiang District, before initiating treatment.
Analysis of WGS data. Paired-end sequencing libraries were constructed with 300-bp inserts, and 150-bp paired-end reads were generated on an Illumina HiSeq 2500 platform. A previously described pipeline was used for calling single nucleotide polymorphisms (SNPs) (27). Briefly, the Sickle tool was used to trim WGS data to keep reads with a Phred base quality above 20 and lengths greater than 30 nucleotides. Reads were mapped to the M. tuberculosis H37Rv reference strain (GenBank AL123456) with Bowtie 2 software (version 2.2.9), and then SAMtools (version 1.3.1) was used for SNP-calling with a mapping quality greater than 30. Varscan (version 2.3.9) was used to identified fixed (frequency, $ 95%) and unfixed (,95%) SNPs with at least 5 supporting reads and the strand bias filter option on. SNPs in repetitive regions of the genome (e.g., the proline-glutamic acid (PE)/proline-proline-glutamic acid (PPE) polymorphic guanine-cytosine-rich sequence (PE/PPE-PGRS) family genes, phage sequences, insertion, or mobile genetic elements) were excluded. The drug resistance profiles and lineages of the M. tuberculosis strains were predicted from WGS data using SAM-TB (33). Multidrug-resistant tuberculosis (MDR-TB) was defined as isolates with mutations conferring resistance to at least isoniazid and rifampicin. Samples with unfixed drug-resistant mutations (DRMs) in at least one locus were considered heteroresistant. Phylogenetic trees were constructed using the maximum-likelihood method (RAxML-NG) (34) and visualized on the Interactive Tree of Life platform (https://itol.embl.de/). Files containing sequencing reads were deposited in the National Institutes of Health Sequence Read Archive under BioProject PRJNA869190.
Detection of mixed infections. WGS data were used to identify mixed infections using the phylogenetic method of Gan et al. (12), which is based on detecting two evolutionary paths. The evolutionary path of an M. tuberculosis strain can be determined by mapping its SNPs onto a reference phylogeny computed from published M. tuberculosis genome sequences (Fig. 1B). The proportion of mixed infections was obtained by calculating the arithmetic mean of the mutation frequencies based on the unique mutations matched in each pathway (Fig. S1).
VNTR typing was performed using a 913 locus combination that has a high resolution for Beijing lineage strains (25). After extraction of bacterial genomic DNA by boiling, nine conventional loci (QUB11b, QUB18, Mtub21, Miru26, QUB26, Mtub04, Miru31, Miru40, and VNTR 2372) and three hypervariable loci (VNTR 3820, 3232, and 4120) were separately amplified by PCR using the primers listed in Table S2. The PCRs for the nine conventional loci were performed in a volume of 10 mL containing 1Â Taq PCR master mix (CoWin Biotech Co. Ltd., Beijing, China), 0.4 mM each primer, and 1 mL DNA template. The reactions for the three hypervariable loci were performed in a volume of 20 mL containing 1Â GC buffer I (TaKaRa Biotech Co. Ltd., Dalian, China), 200 mM each deoxynucleoside triphosphate (dNTP), 0.5 U of Taq (TaKaRa Biotech Co. Ltd.), 0.4 mM each primer, and 1 mL DNA template. The thermocycling conditions were as follows: 95°C for 5 min, followed by 30 cycles at 94°C for 30 s, 58°C (64°C for loci MIRU 40 and VNTR 4120) for 30 s, and 72°C for 30 s (1.5 min for loci VNTR 3232, VNTR 3820, and VNTR 4120), with a final extension at 72°C for 7 min. Based on the PCR results, a locus was considered heterozygous if the PCR product showed at least 2 bands on agarose gel electrophoresis (Fig. 1C). If a sample had two or more heterozygous loci, it was considered a mixed infection (9).
Statistical analysis. The basic demographic and clinical data for the enrolled TB patients were extracted from the routine TB surveillance system. The chi-square test was used to assess possible associations of mixed infections. Univariable and multivariate logistic regressions were used to estimate the odds ratios (ORs) and 95% confidence intervals (CIs). All statistical analyses were performed and visualized in R (version 3.6.0), with P values of ,0.05 considered statistically significant.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 1 MB.