Tuberculosis (TB) is a global public health emergency. Increasingly drug resistant strains of Mycobacterium tuberculosis (M. tb) continue to emerge and spread, highlighting adaptability of this pathogen. Most studies of M. tb evolution have relied on ‘between-host’ samples, in which each person with TB is represented by a single M. tb isolate. However, individuals with TB commonly harbor populations of M. tb numbering in the billions. Here, we use analyses of M. tb genomic data from within and between hosts to gain insight into influences shaping genetic diversity of this pathogen. We find that the amount of M. tb genetic diversity harbored by individuals with TB can vary dramatically, likely as a function of disease severity. Surprisingly, we did not find an appreciable impact of TB treatment on M. tb diversity. In examining genomic data from M. tb samples within and between hosts with TB, we find that genes involved in the regulation, synthesis, and transportation of immunomodulatory cell envelope lipids appear repeatedly in the extremes of various statistical measures of diversity. Many of these genes have been identified as possible targets of selection in other studies employing different methods and data sets. Taken together, these observations suggest that M. tb cell envelope lipids are targets of selection within hosts. Many of these lipids are specific to pathogenic mycobacteria and, in some cases, human-pathogenic mycobacteria. We speculate that rapid adaptation of cell envelope lipids is facilitated by functional redundancy, flexibility in their metabolism, and their roles mediating interactions with the host.
Mycobacterium tuberculosis (M. tb) causes over nine million new cases of tuberculosis (TB) per year and is estimated to infect one-third of the world’s population [1]. The emergence of increasingly drug resistant strains of M. tb [2] demonstrates the bacterium’s ability to adapt to antibiotic pressures, despite limited genetic diversity [3]. Prior research has identified the influence of bottlenecks, population sub-division, and purifying selection on genetic diversity of M. tb circulating among human hosts [4–7]. In these studies, each TB patient was represented by a single M. tb strain isolated in pure culture. However, individuals with TB harbor a large population of M. tb cells for a period of months to years, which raises the possibility of significant diversification of bacterial populations over the course of individual infections. There are few studies of within-host evolution of M. tb. One example is a study of the transposable element IS6110 marker that found multiple lines of evidence suggestive of positive selection on M. tb populations within hosts [8]. Advances in sequencing technologies have since enabled detailed, genome-wide studies of the evolution of intra-host populations of both pathogenic and commensal microbes [9–22]. Whole-genome sequencing (WGS) studies of natural populations of M. tb have focused primarily on the emergence of drug resistance [23–26]. Here, we use analyses of genetic data from within and between patients with TB to characterize M. tb variation across evolutionary scales. We find that overall diversity of M. tb populations within hosts can vary dramatically and identify candidate genetic loci for M. tb adaptation during infection.
We quantified genetic diversity of five within-host populations of M. tb (Table 1) using data from three published studies. The data we selected from the original studies were directly comparable with respect to sequencing platform and other factors we observed to affect diversity (see Methods and S1 Table). The data set includes samples susceptible to all first line drugs (A0 and E0) as well as INH-monoresistant, multidrug-resistant (MDR) and extensively drug-resistant (XDR) samples (Table 1). HIV status was not reported in any of the primary studies. M. tb lineage 2 (East Asian, Patients A, B, C, and D) and lineage 4 (Euro-American, Patient E) are included in this sample; patient origin was reported as China (Patients A, B, and C), Abkhazia (Patient D) and “Eastern Europe” (Patient E). We used PoPoolation software [27] to estimate two standard measures of genetic diversity, nucleotide diversity (π) and Watterson’s theta (ƟW) (see Box 1 for descriptions), from the pooled within-host M. tb genomic data. It is challenging to distinguish rare genetic variants from sequencing errors in pooled genomic data: PoPoolation implements methods that account for the effects of sequencing errors on low frequency variant calls [27,28]. In order to eliminate effects of coverage differences on diversity, we sub-sampled the genomic data to a uniform 50X coverage (see Methods). We compared the pooled M. tb within-host data to a globally diverse sample of 201 ‘between-host’ genomes from the seven major lineages, where each individual with TB is represented by a single isolate of M. tb (S2 Table). These data were originally reported in Comas et al [29]. Drug resistance phenotypes were not reported in the original study. We examined the between-host alignment at 1004 loci with a reported association with drug resistance [30] and found that 37 of these loci had segregating polymorphisms. We conclude that the between-host dataset is likely to include both susceptible and drug resistant isolates. Estimates of M. tb diversity within hosts were very sensitive to coverage and platform; measures of diversity from some samples were also sensitive to base quality score and minor allele count threshold (S1 File and S3 Table). Patients D and E were similar in that both developed progressive drug resistance during treatment, and therapy was tailored to results of extended drug susceptibility testing (DST). Patient E’s treatment outcome was eventually successful; treatment outcome of Patient D was not reported. For a given set of parameters, estimates of diversity from Patients D (M. tb lineage 2) and E (M. tb lineage 4) were very similar to each other, and stable across serial samples (Fig 1 and S3 Table). Estimates from Patient E were slightly more sensitive to the minor allele count threshold, suggesting this patient harbored more rare variants than Patient D. Patient A (M. tb lineage 2) was described as non-adherent during standard therapy of drug-susceptible disease; the second and third samples from this patient exhibited INH-monoresistance. Patient A was treated with two months of four-drug treatment, followed by 27 months of INH and RIF (i. e. only one drug to which the bacteria were susceptible). Measures of diversity from the M. tb population of Patient A were similar overall to those from Patients D and E, but showed an increasing trend over the sampling period (Fig 1). Increasing diversity could indicate loss of control of the infection. Patient A defaulted treatment, and no further samples were reported after the sample collected at 24 months into treatment. Patients B and C had MDR TB. Extended DST results were not reported for these patients and both died. M. tb diversity of pre-terminal samples from Patients B and C was extremely high relative to the other intra-host samples (Fig 1). This difference in diversity could be due to a technical issue, such as relatively high sequencing error rates in these data. In this case, we would expect that application of more stringent quality filters would decrease observed diversity. Application of a higher base quality threshold did not result in loss of relative diversity of these samples. Differences between samples from Patients B, C, and the others were, however, less marked when the minor allele count threshold was increased (S3 Table). This suggests that high diversity of B36 and C11 was driven by rare variants. Excess rare variants are a hallmark of an expanding population; diversity of pre-terminal samples from Patients B and C could reflect expansion of the M. tb population in the final phases of an uncontrolled infection. Another possibility is that terminal progression of their TB infections involved extensive breakdown of lung tissue allowing sampling of previously inaccessible M. tb sub-populations. TB patients included in our within-host data set were culture positive for prolonged periods: this is atypical for settings with well-functioning TB control programs and these patients are unlikely to be broadly representative. Further studies of within-host diversity of TB patients from a range of settings and with a variety of clinical presentations are needed to fully characterize M. tb adaptation within hosts. Nonetheless, these data delineate some interesting patterns of within-host diversity. It is striking that M. tb within-host diversity is so similar across patients from three independent studies (Patients A, D, and E). In vitro estimates of M. tb’s mutation rate vary substantially according to lineage [31], yet we did not observe an obvious effect of lineage on diversity in comparisons of M. tb populations from lineage 2 (Patients A and D) and lineage 4 (Patient E). This could reflect differences between in vitro and in vivo mutation rates, or perhaps other parameters are more important in shaping overall patterns of diversity within hosts. Diversity of M. tb lineage 2 was lower than lineage 4 at the between-host scale (Fig 1), whereas lineage 2 isolates have been observed to evolve more quickly in the lab. Our data set includes four within-host samples that were collected prior to the initiation of treatment (A0, C1, D0, and E0). We did not observe a substantial decrease in M. tb diversity following initiation of treatment in any of these patients (A19, C11, D9, and E8). It could be that TB treatment does not affect the average number of pairwise differences and segregating sites of resident M. tb populations. Alternatively, measures of M. tb diversity may change in response to treatment in TB patients whose sputum cultures convert more quickly than the patients in our study. Our study included two M. tb samples with susceptibility to all first line agents (A0 and E0). Diversity of these samples was not markedly different from samples with a drug resistant phenotype. In summary, our results suggest that M. tb lineage, initiation of TB treatment, and drug resistance do not have strong impacts on diversity measures for within-host M. tb populations. Disease severity, on the other hand, appears to have marked effects on M. tb diversity. Several studies have identified distinct patterns of variation at M. tb genetic loci associated with different functions [4,32]. Such patterns of variation may reflect distinct regimes of natural selection, heterogeneity of mutation rates, or other influences. We sought to identify bacterial genetic loci with extreme patterns of variation at within- and between-host scales; these are candidate loci for M. tb adaptation during transmission and infection. For each gene with data that passed our quality thresholds (see Methods), we quantified π and ƟW using methods that account for sequencing error [27]. We used the same approach and estimated two statistics designed to capture deviations from neutral patterns of variation: Tajima’s D (TD, [33]) and the ratio of non-synonymous changes per non-synonymous site (πN) to synonymous changes per synonymous site (πS) (see Box 1). We disregarded genes lacking either synonymous or non-synonymous variation in comparisons of πN/πS. Since neutral forces such as population growth can affect patterns of variation, we compared relative values of statistics in each within-host sample, where all genes are likely to have the same demographic history. Natural selection can lead to population differentiation when the relative fitness of genotypes varies among environments; empirical outlier analyses for loci with extreme measures of population differentiation are commonly used to identify candidates of positive selection [34–36]. Treating each serial sputum sample from an individual patient as a distinct M. tb population, we calculated pairwise FST values for polymorphic sites covered by at least 10 reads with a minimum minor allele count of 6 (pooled across all samples from a patient) with PoPoolation2 [37]. In order to reduce biases introduced by variable coverage, we conditioned our analysis on the ability to detect a significant change in allele frequency between samplings using a two-sided Fisher’s exact test as previously proposed [38] (see Methods). We have shown that numerous genes affecting PKS associated lipids (mycolic acids, PDIM, PAT/DAT, PGL, POL, and MPM) exhibit extreme patterns of variation across a range of statistical measures and independent data sets. We hypothesize that these patterns of variation are due to adaptation within M. tb’s natural human host. Several characteristics of PKS-associated lipids make them likely targets of selection during infection. PKS associated lipids are important mediators of host-pathogen interactions. For example, mycolic acids are known to modify host cell phenotype [64] and to affect virulence in animal models of TB [65]. PDIM has major impacts on M. tb’s virulence in animal models of TB [66,67] via numerous mechanisms (reviewed in [68]) including protection from innate immune responses [69]. Like PDIM, the PKS synthesized lipids PAT/DAT, SL, MPM, POL and PGL also appear to mediate interactions between M. tb and the immune system [66,70–86]. The capacity for adaptation is predicated on phenotypic variation, which has been demonstrated for PKS associated lipids. Production of PAT/DAT has been shown to vary among clinical isolates; some isolates do not make them at all [87]. Production of PGL also varies among clinical isolates, and variant forms of PDIM have been identified in clinical isolates [54]. Based on their roles in pathogenesis and their variability, it has been proposed previously that modifications of these lipids could arise during natural infection in order to optimize host cell manipulation [54]. The observation of unusual patterns of variation at genes involved in the synthesis, transportation, and regulation of these lipids provides evidence in support of this hypothesis. Several features are conducive to efficient selection on genes controlling the synthesis, regulation, and transportation of PKS associated lipids. There are connections between biosynthetic pathways controlling production of distinct PKS synthesized lipids (e. g. PDIM, SL, and mycolic acids) such that metabolites can be shuttled between them [88,89]. This flexibility allows M. tb to respond rapidly to environmental fluctuations. It may also allow more efficient selection for adaptive mutations, since single mutations can potentially affect multiple lipid products. In addition, since intermediate metabolites could potentially be shuttled down multiple pathways, they need not accumulate with potentially toxic effects if one pathway is affected by a harmful mutation. PDIM, PAT/DAT, and SL are not essential for growth of M. tb in artificial media, and PDIM is in fact frequently lost during passage of M. tb in the laboratory [87,90]. This suggests that functions of these lipids are specific to natural, within-host environments. The expression of M. tb immunomodulatory lipids is responsive to physiological conditions (hypoxia, starvation) encountered during natural infection [91,92]. In tissue culture and animal models of TB, the shift from axenic growth to infection is accompanied by alterations in expression of genes involved in lipid metabolism with consequent changes to the bacterial cell envelope [88,93–95]. Further supporting the idea that these lipids are important for adaptation to the pathogenic niche, MPM, PGL, and PDIM are only found among pathogenic mycobacteria, and PAT/DAT and SL are specific to the Mycobacterium tuberculosis complex (MTBC) (Fig 6) [54,62,79]. It was shown recently that even among members of the MTBC, production of PAT/DAT and SL is specific to human-pathogenic mycobacteria [62]. Recent work has shown that PAT/DAT, PDIM, and SL all impair phagosomal acidification and thereby improve M. tb survival within macrophages; there appears to be significant flexibility in how these functions are performed, and distinct lipid moieties may compensate for each other [87]. By providing cover during exploration of the fitness landscape, both functional redundancy and metabolic flexibility may increase the potential for rapid adaptation. For example, they would enable progress toward fitness peaks that are otherwise unreachable as a result of a preceding ‘valley’ created by mutations that are deleterious except in certain combinations. This may be particularly important for a clonally reproducing organism such as M. tb. We have shown that genetic diversity of M. tb populations can vary dramatically between individuals with TB. Variation correlated with clinical severity, rather than treatment status, drug resistance phenotype or M. tb lineage. Further studies are needed to explore the replicability of patterns observed here across different types of patients with TB and treatment settings. We have also observed distinct patterns of variation among genes associated with M. tb virulence lipids: these patterns carry through functional links among genes, and across different statistical measures and datasets in our own and previously published studies. These lipids are good candidates for adaptation during infection based on the phenotypic variation observed among clinical isolates, their roles mediating host interactions, their specific association with the pathogenic niche, and their flexible functional architecture. It is important to note uncertainties in the data and its interpretation. We cannot, for example, assert that all genetic variants described here were generated within hosts. In all the studies whose data we analyzed, clinical specimens were manipulated in vitro prior to M. tb genome sequencing. Some mutations may have occurred during culture. Genetic variation in M. tb is skewed to rare allele frequencies, at both the within- and between-host scale. While we have endeavored to ameliorate the impacts of spurious variants introduced by sequencing errors and errors introduced during processing of sequencing read data, rare variants are difficult to disentangle from sequencing errors, particularly in pooled sequencing data. As is standard in studies of M. tb genomics, our analyses excluded genomic segments (totaling ~11%, see Methods) that are poorly resolved with short read sequencing technologies. Additional measures to reduce the effects of sequencing errors included trimming low quality bases from sequencing reads, implementing mapping and base quality thresholds, and excluding indels as well as SNPs demonstrating strand-bias and/or tail-distance-bias. We used PoPoolation, which explicitly accounts for the effects of sequencing error, to estimate population genetic parameters in our study. Validation studies of PoPoolation have shown that implementation of PoPoolation’s recommended minor allele count and quality thresholds decreases the error rate from 1% to ~ 0. 01% post processing [27]. These studies also demonstrated that PoPoolation parameter estimates from pooled sequencing data are highly correlated with estimates from Sanger sequencing data. However, there are differences between PoPoolation validation datasets and those analyzed in our study: PoPoolation was validated for sequencing reads generated on the Illumina GAIIx from a Drosophila melanogaster population (and data simulated to mimic these reads), whereas we have analyzed sequencing reads generated on the Illumina HiSeq from M. tb populations. GAIIx and HiSeq differ in their biases [96], and GC content differs significantly between M. tb and D. melanogaster. We note that despite these potential limitations of PoPoolation in correcting errors in our pooled data, we observed similar patterns of genetic variation in the between-host dataset, which is based on single-colony isolates and a different variant-discovery pipeline. Interpretation of genetic data across evolutionary scales is also not straightforward. The consistent identification of specific functional categories in extremes of various statistics suggests that they are due to natural selection, but the type of selection and where it operates is not clear. The patterns we observed may be due to selection imposed by drug therapy, the immune response and/or transmission, with further modification as a result of complex neutral influences on variation. In addition, our hypotheses need to be investigated at the level of phenotypes. Based on the analyses presented here and their contextualization with published data, we propose that further characterization of individual-level, phenotypic variation in M. tb PKS associated lipids is likely to be fruitful. In addition, several of the loci described here are worthy of investigation for a potential role mediating or compensating for drug resistance.
For the within-host M. tb data set, we trimmed low-quality bases from the FASTQ data using a threshold quality score of 20, and reads of length less than 35bp were discarded using Trim Galore! (http: //www. bioinformatics. babraham. ac. uk/projects/trim_galore/) —a wrapper tool around Cutadapt [99] and FastQC (http: //www. bioinformatics. babraham. ac. uk/projects/fastqc/). We mapped reads to H37Rv (NC_000962. 3) [100] using the default parameters of the BWA MEM algorithm with the–M flag enabled for downstream compatibility [101], and we removed duplicates using Picard Tools (http: //picard. sourceforge. net). Local realignment was performed with the Genome Analysis ToolKit (GATK) [102] in a sample aware manner, and aligned reads with mapping quality > 20 were converted to mpileup format using Samtools software [103]. The resulting reference-guided assembly of each sample spanned over 98. 9% of the H37Rv genome, with a mean depth of coverage per site ranging from 56X to 192X (S1 Table). Loci at which indels were present in a sample were removed along with 5bp of flanking sequence using the PoPoolation package [27]. For the between-host M. tb data set, we trimmed low-quality bases from FASTQ data using a threshold quality of 15, and reads resulting in less than 20bp length were discarded using Trim Galore! (http: //www. bioinformatics. babraham. ac. uk/projects/trim_galore/) —a wrapper tool around Cutadapt [99] and FastQC (http: //www. bioinformatics. babraham. ac. uk/projects/fastqc/). Reads were mapped to H37Rv (NC_000962. 2) [100] using the suggested algorithm of BWA for the given sequencing strategy (e. g. paired-end/single-end, read length) [101,104], and duplicates were removed using Picard Tools (http: //picard. sourceforge. net). We used GATK to perform local realignment and variant calls using a minimum phred-scaled confidence threshold of 20 [102]. Variants were filtered with the following expression “QD < 2. 0 || FS > 60. 0 || MQ < 40. 0 || MQRankSum < -12. 5 || ReadPosRankSum < -8. 0” as described in the GATK best practices. Genome alignments were generated with scripts that can be found at https: //github. com/tracysmith/RGAPepPipe. Transposable elements, phage elements, and repetitive families of genes (PE, PPE, and PE-PGRS gene families) that are poorly resolved with short read sequencing were removed from the mpileup files and alignment prior to subsequent analyses plus 5bp up- and down-stream of the genes. We additionally removed regions that were found to have poor mapping quality using the CallableLoci tool of the GATK: for each within-host sample, any region reported as poorly mapped using the following flags was removed from all datasets (including between-host) plus 5bp up- and down-stream: -frlmq 0. 04 –mmq 20 -mlmq 19. As an additional measure, we also removed sites with the lowest 2% average mapping quality in the global alignment from all datasets prior to subsequent analysis. Regions removed from all datasets based on these criteria can be found in S8 Table. Polymorphisms in the within-host dataset demonstrating strand-bias or tail-distance-bias identified by a previously described method [22] (https: //github. com/tamilieberman/IntrasamplePolymorphismCaller) were removed along with 5bp up- and down-stream from all within-host samples (S8 Table). We used the PoPoolation package [27] to estimate nucleotide diversity (π), Watterson’s theta (ƟW), and Tajima’s D (TD) in sliding windows across the M. tb genome. Sensitivity analyses of these statistics to input parameters are described in the S1 File. Parameter estimates were strongly influenced by sequencing coverage (S1 File). In order to alleviate biases in parameter estimates caused by variable coverage among samples, we randomly sub-sampled read data without replacement to a uniform coverage of 50X; this process was performed 9 times for each sample to reduce potential biases introduced by sampling of rare alleles. We used default equations in the PoPoolation package to estimate π, ƟW, and TD [27,28]. In order to reduce the effects of sequencing errors on parameter estimates, PoPoolation implements modified calculations of the classical π, ƟW, and TD that are only evaluated on SNPs above a designated minor allele count. PoPoolation estimators account for the truncated allele frequency spectrum and re-sequencing of the same chromosomes that occurs during pool-seq. PoPoolation’s recommended minor allele count and quality thresholds are based on simulation studies indicating that their implementation decreases the error rate of 1% in the raw data to ~0. 01% in the processed data. PoPoolation parameter estimates from pooled data have also been found highly concordant with estimates from Sanger sequencing data. Unless otherwise noted, all calculations performed with the PoPoolation package were implemented using the recommended minimum minor allele count of 2 for 50X coverage, and a pool-size of 10,000 (see S1 File for justification). Calculations performed under differing parameters can be found in the Supplemental Data. Sliding-window analyses were performed using a window-size of 100K and a step-size of 10K; data presented are the mean of each window across 9 replicate sub-sampled mpileups (Fig 1). Genome-wide averages are expressed as the mean across all windows sufficiently covered (≥ 60% coverage of 50X) (S3 Table). To calculate gene-by-gene estimates of π, ϴW, and TD from the 50X sub-sampled mpileups, we again used the PoPoolation package [27]. Gene coordinates were obtained from the standardized Gene Transfer Format (. gtf) for the H37Rv annotation on the TB database (tbdb. org) (S2 File) [105]. Excluding genes with inadequate coverage (<50% of the gene), we calculated the mean value of each statistic across the 9 replicate sub-samplings for each gene of each sample, and compared it to all other genes within the sample. Using the same. gtf file and the PoPoolation package, we also calculated the average number of nonsynonymous differences per nonsynonymous site (πN) and the average number of synonymous differences per synonymous site (πS) for each gene in the 50X sub-sampled mpileups. Recognizing that chance samplings of very rare mutations in one replicate sub-sampled mpileup would lead to skewed distributions, we took the median values of πN and πS across 9 replicate sub-samplings for each gene of each sample and calculated πN/πS. Excluding genes with inadequate coverage (<50% of the gene) and genes with πN and/or πS equal to zero, πN/πS values were compared relative to all genes passing these criteria within the sample. To identify candidate SNPs under selection, we treated each temporal sample of a patient as a population, and estimated pairwise FST for each variable site in the genome with PoPoolation2 [37]. Unlike the sliding-window and gene-wise estimates, all sequencing data for which the minimum base quality was ≥ 20 was considered (i. e. no sub-sampling). We excluded all sites with a coverage ˂ 10, and only considered those sites with a minimum allele count ≥ 6 (pooled across all samples for a patient). FST estimates were subjected to an empirical outlier analysis. To reduce biases resulting from variable coverage, we conditioned our analysis on the ability to detect a significant change in allele frequency between samplings using a two-sided Fisher’s exact test as previously proposed [38]. To account for the large number of tests performed we used a false discovery rate (FDR) of 5% and calculated adjusted p-values (q-values) (R stats Package [106]). SNPs with extreme FST values (> 0. 10) and a q-value < 0. 01 were deemed outliers and are listed in S5 Table. S9 Table shows SNPs passing all filtering criteria that were used in the above described population genetic estimates for the within-host populations. We used PoPoolation [27] with the “disable-corrections” flag enabled (calculations are performed using the classical equations) to generate sliding-window estimates of π and ƟW from whole genome alignments of all 201 globally extant M. tb strains, M. tb isolates from lineage 2 (East Asian lineage, n = 37), and M. tb isolates from lineage 4 (Euro-American lineage, n = 53) (S2 Table) [29]. We required a minimum of 75% of the strains have non-missing data for a site to be included in the analysis (“minimum coverage”), set the “minimum minor allele count” to one, and the “pool-size” to the number of strains being analyzed. We calculated gene-by-gene estimates of π, ϴ, TD, πN, and πS for the global dataset in the same way as was done for within-host samples, save for a few exceptions: rather than use the default PoPoolation estimators that apply corrections for sequencing errors [27,28], we employed the “disable-corrections” flag (calculations are performed using the classical equations), set the “minimum minor allele count” to one, and the “pool-size” was set to 201 to reflect the number of strains in the dataset; the “minimum coverage” was set to 151 to exclude any genes that were not represented by at least 75% of strains in the dataset. Finally, no averaging was performed. For πN/πS calculations, we again excluded genes with πN and/or πS equal to zero, and compared values relative to all genes passing these criteria within the sample. For each within-host sample, between 58. 7–90. 1% of annotated genes in the H37Rv genome had sufficient coverage to calculate TD; for the between-host sample, TD was calculated for 90. 2% of annotated genes in the H37Rv genome, as this was the fraction covered by at least 75% of the strains. Genes with TD values in the top and bottom 5% of the distribution for a given sample were deemed candidate genes of selection. The significance of enrichment for functional categories in candidate genes of selection was assessed with a two-sided Fisher’s exact test. To account for multiple hypothesis testing, we used a false discovery rate of 5% and calculated q-values (R stats Package, [106]). We used the following annotation categories to classify M. tb genes: computationally predicted Clusters of Orthologous Groups (COG) (n = 21 categories) [40]; essential and nonessential genes for growth in vitro as determined by transposon site hybridization (TraSH) mutagenesis [107]; genes essential for growth in a murine model of TB [108];" targets of independent mutation" associated with drug resistance [44]; and the M. tb-specific, manually curated functional annotation lists from TubercuList (n = 7) [39]. COG annotations for the H37Rv genome were obtained from the TB database (tbdb. org) [105], as were the TraSH “in vitro essential”, “in vivo essential”, and “nonessential” gene annotations. TubercuList functional annotations were obtained from tuberculist. epfl. ch [39] and reflect the most up-to-date annotations when the database was accessed (12/01/2013). We did not include the TubercuList categories “PE/PPE”, “stable RNA”, and “unknown” in our analyses. S3 File contains the genes included in each category. Functional enrichment of genes in the top 5% of π and ϴW for each sample were performed in the same manner. Variable proportions of annotated genes in the H37Rv genome had sufficient coverage and contained both nonsynonymous and synonymous variation in a given within-host sample. To look for commonalities among patients, we examined all genes found to have a positive πN/πS value in at least one sample from multiple patients; genes found to have a πN/πS > 1 in three patient M. tb populations are shown in Table 2. Upon noticing a preponderance of genes in the superpathway of mycolate biosynthesis on the cellular overview tool of the TB database (tbdb. org) [105] we formally tested for an enrichment of genes in this category using a two-sided Fisher’s exact test. We used the Python package pysamstats (https: //github. com/alimanfoo/pysamstats) to calculate the root-mean-square (RMS) value of base qualities for variant alleles and the RMS mapping quality for reads aligned at such polymorphic sites in the reference-guided-assemblies of within-host samples. Variants occurring in genes categorized as CHP, LIP, or COG: Q were subject to a one-sided Student’s T-test to determine whether the mean RMS base qualities of the category was significantly lower than that occurring in any gene. We downloaded finished genomes of mycobacterial species from NCBI (S7 Table). Reference guided assemblies for species where only sequencing reads were available were performed as described above for the between-host dataset. We used Prokka v 1. 7 [109] for genome annotation. Protein sequences output by Prokka were clustered into orthologous groups using OrthoMCL [110]. The core proteins (those found only once in every genome) were aligned using MAFFT [111], trimmed with trimAl [112], and concatenated. Scripts used for core genome analysis can be found at https: //github. com/tatumdmortimer/core-genome-alignment. We used RAxML v 8 [113] for maximum likelihood phylogenetic analysis of the core alignment. We used Dendroscope [114] for tree viewing and editing.
Tuberculosis (TB) is a grave threat to global public health and is the second leading cause of death due to infectious disease. The causative agent, Mycobacterium tuberculosis (M. tb), has emerged in increasingly drug resistant forms that hamper our efforts to control TB. We need a better understanding of M. tb adaptation to guide development of more effective TB treatment and control strategies. The goal of this study was to gain insight into M. tb evolution within individual patients with TB. We found that TB patients harbor a diverse population of M. tb. We further found evidence to suggest that the bacterial population evolves measurably in response to selection pressures imposed by the environment within hosts. Changes were particularly notable in M. tb genes involved in the regulation, synthesis, and transportation of lipids and glycolipids of the bacterial cell envelope. These findings have important implications for drug and vaccine development, and provide insight into TB host pathogen interactions.
Diversity of Mycobacterium tuberculosis across Evolutionary Scales