Signatures of Selection at Drug Resistance Loci in Mycobacterium tuberculosis

Mycobacterium tuberculosis, the causative agent of tuberculosis (TB), is a significant burden on global health. Antibiotic treatment imposes strong selective pressure on M. tuberculosis populations. Identifying the mutations that cause drug resistance in M. tuberculosis is important for guiding TB treatment and halting the spread of drug resistance. Whole-genome sequencing (WGS) of M. tuberculosis isolates can be used to identify novel mutations mediating drug resistance and to predict resistance patterns faster than traditional methods of drug susceptibility testing. We have used WGS from natural populations of drug-resistant M. tuberculosis to characterize effects of selection for advantageous mutations on patterns of diversity at genes involved in drug resistance. The methods developed here can be used to identify novel advantageous mutations, including new resistance loci, in M. tuberculosis and other clonal pathogens.

M ycobacterium tuberculosis, the causative agent of tuberculosis (TB), is estimated to have caused 1.4 million deaths in 2015, making it the leading cause of death due to an infectious disease. The proportion of TB due to multidrug-resistant (MDR) M. tuberculosis resistant to the first-line antituberculosis drugs isoniazid (INH) and rifampin (RIF) is increasing (1), which poses a major threat to global public health. Unlike most bacteria, M. tuberculosis evolves clonally, so resistance cannot be transferred among strains or acquired from other species of bacteria: drug resistance in M. tuberculosis results from de novo mutation within patients and transmission of drug-resistant clones (2)(3)(4). The relative contributions of de novo emergence and transmitted drug resistance vary across sampling locations (5)(6)(7)(8)(9). Another potential variable in the emergence of drug-resistant TB is M. tuberculosis's lineage structure: seven distinct lineages have been identified among globally extant populations of M. tuberculosis. Of these lineages, lineage 2 (L2) has been associated with relatively high rates of drug resistance, and it has been postulated that the acquisition of resistance is a result of higher rates of mutation in this lineage (10). Studies of M. tuberculosis evolution within hosts with TB have shown that the emergence of drug resistance is associated with clonal replacements that lead to reductions in genetic diversity of the bacterial population (11,12).
Many of the methods developed to identify advantageous mutations, such as those conferring antibiotic resistance, depend on recombination to differentiate target loci from neutral variants (13). However, in clonal organisms like M. tuberculosis, neutral and deleterious mutations that are linked to advantageous variants will evolve in tandem with them. This linkage among sites can also cause competition between genetic backgrounds with beneficial mutations, decreasing the rate of fixation for beneficial alleles, while deleterious alleles are purged less efficiently (14)(15)(16).
While the majority of the M. tuberculosis genome is subject to purifying selection (i.e., selection against deleterious mutations) (3), antibiotic pressure exerts strong selection for advantageous variants that confer resistance. M. tuberculosis drug resistance has been the focus of extensive investigation, and a variety of resistance mutations have been characterized for commonly used antituberculosis drugs (17). Drug resistance mutations can be associated with fitness costs (18)(19)(20), and compensatory mutations that ameliorate these fitness costs have been identified in the context of rifampin resistance (21,22). Resistance mutations found to have lower fitness costs in vitro-as measured by competition assays-are found at higher frequencies among M. tuberculosis clinical isolates and appear to be transmitted at higher rates relative to mutations with high in vitro fitness costs (18,23). Candidate loci involved in resistance and compensation for its fitness effects have been identified previously by screening for homoplastic variants (i.e., mutations that emerge more than once in the phylogeny) that are significantly associated with drug-resistant phenotypes (24) and genes with higher dN/dS (ratio of nonsynonymous mutations to synonymous mutations) in resistant compared to sensitive isolates (25). Application of these methods to wholegenome sequence data from M. tuberculosis clinical isolates has recovered known drug resistance loci, as well as loci associated with cell surface lipids and biosynthesis, DNA replication, and metabolism.
The goal of the present study was to use patterns of genetic diversity at known drug resistance loci to identify the population genomic signatures of positive selection in natural populations of M. tuberculosis. Using whole-genome sequence data from two populations for which phenotypic resistance data were available, we have identified several distinct signatures associated with these loci under selection. On the basis of these results, we propose methods of identifying loci under positive selection, including novel drug resistance loci, in clonal bacteria such as M. tuberculosis.

RESULTS
We inferred the phylogeny of 1,161 M. tuberculosis isolates from Russia and South Africa (see Materials and Methods and see Table S1 in the supplemental material) using the approximate maximum likelihood method implemented in FastTree 2 (Fig. 1). The majority of the isolates belong to lineage 2 (L2) (n ϭ 667) and L4 (n ϭ 481). M. tuberculosis nucleotide diversity was similar to previous estimates from a globally distributed sample (26). We identified lineage-specific patterns in overall diversity, with L4 having higher diversity than L2 (average number of pairwise differences of L2 [ L2 ], 3.6 ϫ 10 Ϫ5 , L4 , 1.5 ϫ 10 Ϫ4 ). Previously published analyses of whole-genome sequence data from L2 indicate that the majority of L2 isolates worldwide belong to a sublineage that has undergone relatively recent expansion (27,28). In this sample from Russia and South Africa, the majority of L2 isolates belong to this sublineage, while the L4 isolates are associated with deeper branching sublineages. This likely contributes to the observed differences in diversity.
The overall diversity of L2 was lower than that of L4 in our sample ( Fig. 2, P Ͻ 2.2 ϫ 10 Ϫ16 ). A total of 762 of the isolates in our sample (66%) are resistant to one or more antituberculosis drugs (Table 1). Drug-resistant TB can be acquired as a result of de novo mutations within a patient or by infection with a resistant strain. When resistance is primarily mediated by de novo mutations, diversity should be similar in susceptible and resistant populations, as resistance will arise on multiple genetic backgrounds. In contrast, if resistance develops primarily via transmission of resistant strains, the resistant subpopulation should be less diverse than the susceptible subpopulation. We compared the nucleotide diversity of susceptible and resistant subpopulations and found genome-wide estimates of nucleotide diversity to be higher in isolates susceptible to a range of drugs for which phenotyping data were available (P ϭ 0.029 by paired t test). In comparisons of gene-wise diversity in susceptible and resistant populations, we found that resistant isolates had a greater number of genes with no diversity, but levels of diversity within genes in which it was measurable were similar between resistant and susceptible populations (Fig. 3).
Of the 3,162 genes included in our analyses, 109 (3%) were invariant across all isolates in our sample. This is likely due to strong purifying selection on these genes. An additional 136 genes harbored variation in the drug-susceptible populations but were invariant across all of the drug-resistant populations. We did not observe the converse, i.e., genes that were invariant in susceptible isolates specifically, which supports the conclusion from genome-wide diversity estimates that resistant isolates represent a subset of the diversity found in susceptible populations and suggests that there may be purifying selection that is specific to the setting of drug resistance. In order to evaluate  Among genes in which nucleotide diversity is measurable, it is similar between resistant and susceptible isolates even when drug resistance-associated genes and targets of independent mutation identified by Farhat et al. (24)  whether the observed pattern was likely to arise by chance, we performed weighted random sampling of genes. The weighting was based on diversity in susceptible populations, assuming that genes with low diversity in susceptible populations are more likely to be invariant in resistant populations. After randomly sampling genes in each drug-resistant population 1,000 times, we found that no samples contained shared genes among all resistant populations (first-and second-line drugs). This suggests that specific genes tend to lose diversity in the setting of drug resistance, which could result from purifying selection specific to this setting. A potentially important caveat is that in our data set, drug-resistant populations are not independent, and the same isolates are often resistant to multiple drugs. Since resistance to second-line drugs frequently arises on genetic backgrounds already resistant to one or more first-line drugs, we repeated the sampling with only first-line drugs and found that the maximum number of sampled genes shared across all populations was 2 (median, 0). Overall, these results suggest that certain genes are more likely to lose diversity as drug resistance evolves, but we cannot completely rule out the possibility that the pattern arose as a result of overlapping membership in resistant populations. We compared the diversity of drug resistance-associated genes ( Table 2) with the rest of the genome using two measures of diversity: average number of pairwise differences () and number of segregating sites (). We found the resistance genes gid, rpsL, and pncA to be in the top fifth percentile of gene-wise and/or values. rrs and ethA are in the top fifth percentile of , but not . Surprisingly, despite katG being a target of multiple drug resistance mutations (Table 2), we did not identify extreme levels of diversity in katG (80th and 82nd percentiles of and , respectively).
We also examined gene-wise diversity values within each lineage to look for lineage-specific high-diversity genes. In both L2 and L4, gid, rpsL, pncA, ethA, and thyA were in the top fifth percentile of diversity ( and/or ). In L2, rpoB, embB, Rv1772, and folC were additionally in the top fifth percentile of gene-wise and/or values. In L4, Rv0340 was in the top fifth percentile of gene-wise and/or . While rpoB and embB were not in the top 5th percentile of gene-wise in L4, they still had high diversity (91st and 82nd percentiles, respectively). The lineage-specific differences in diversity of Rv1772, folC, and Rv0340 suggest that there are interactions between these loci and loci that differentiate L2 and L4.
We used gene-wise estimates of Tajima's D values to investigate gene-specific skews in the site frequency spectrum that could result from selection, where negative values indicate an excess of rare variants and positive values indicate an excess of intermediate frequency variants. We previously identified a relationship between gene length and gene-wise estimates of Tajima's D values for M. tuberculosis (26), and this finding was corroborated here (R 2 ϭ 0.3 after log 2 transformation). In order to identify genes with extreme values of Tajima's D values-out of proportion with their length-we performed linear regression on log 2 -transformed gene lengths and Tajima's D values and identified genes with the largest residuals (Fig. 4). pncA, ethA, and embC all had Tajima's D values lower than expected based on their length (fifth percentile of residual values). This indicates that these genes contain an excess of rare variants compared to other genes in the genome. Excess rare variants can result from a population expansion, a selective sweep, or purifying selection.
We calculated the ratio of and of resistance-associated genes in isolates susceptible and resistant to first-line drugs and identified genes with markedly different diversities in resistant and susceptible subpopulations (Fig. 5A). Among resistance genes in the top fifth percentile of gene-wise and overall, diversity of pncA and ethA is relatively high among resistant isolates, whereas diversity of gid is similar in resistant and susceptible populations. We also examined differences in this ratio between isolates in L2 and L4 (Fig. 5B). Rv1772 and embR were more diverse in resistant isolates in L2, and kasA and tlyA were more diverse in resistant isolates in L4.
We used fixation index (F ST ) outlier analysis to identify single nucleotide polymorphisms (SNPs) and indels that exhibited extreme differences in frequency between susceptible and resistant populations. Our a priori expectation was that variants me-diating resistance would be at markedly higher frequency in the drug-resistant subpopulation and that drug targets would be enriched among genes harboring variants with high F ST . After removing SNPs in regions corresponding to indels and variants at sites missing data for greater than 5% of isolates, the highest F ST SNPs in comparisons of resistant and susceptible subpopulations to first-line drugs are in katG (2155168, F ST ϭ 0.89, INH), rpoB (761155, F ST ϭ 0.72, RIF), and rpsL (781687, F ST ϭ 0.37, streptomycin [STR]). These SNPs were also F ST outliers in the lineage-specific analyses. We used a randomization procedure to assess the significance of observed F ST values and found the maximum F ST values after randomly assigning resistant and susceptible designations to be 0.023 for INH, 0.019 for RIF, and 0.018 for STR. In addition to SNPs within known drug resistance-associated genes, we identified F ST outliers in genes that may be novel targets for drug resistance (Table 3).
Homoplastic SNPs-i.e., SNPs that evolve more than once in a phylogeny-are candidate loci under positive selection and have previously been used to identify resistance-associated mutations in M. tuberculosis (24). Of the 235 genes with homoplastic SNPs that we identified in our sample, 13 are known to be associated with drug resistance (Fig. 6), and resistance genes were significantly enriched among genes with homoplastic SNPs (P ϭ 1.2 ϫ 10 Ϫ4 by Fisher's exact test). pncA had the largest number of homoplastic SNPs of any gene in the genome (27 distinct SNPs that appear Ͼ1 in the phylogeny). The SNPs identified in F ST analysis were also identified as homoplastic (Fig. 6). Our results suggest that complementary approaches based on homoplasy and F ST outlier analysis can be used to identify SNPs associated with a trait of interest (in this case drug resistance). In addition to genic SNPs, we observed homoplastic SNPs that are also F ST outliers in intergenic regions upstream of drug resistance-associated genes (Table 3). These SNPs are candidate resistance and compensatory mutations with a regulatory mechanism of action. In our analyses of indels, we controlled for the possibility that indels affecting the same gene may not be called in exactly the same position by considering indels within the same gene as identical. We identified four drug resistance-associated genes with homoplastic indels: gid, ethA, rpoB, and pncA. F ST values for the deletion in gid were in the top fifth percentile for populations resistant to capreomycin (CAP), ethambutol (EMB), ethionamide (ETO), kanamycin (KAN), ofloxacin (OFL), and pyrazinamide (PZA), but interestingly, the deletion was not associated with STR resistance (F ST ϭ 0.04). Unlike homoplastic SNPs, homoplastic indels were not significantly enriched for drug resistance-associated loci (P ϭ 1).
We recovered 20 out of 40 known drug targets by identifying genes with extreme values of diversity, homoplastic SNPs, or SNPs that are F ST outliers in comparisons of resistant and susceptible subpopulations. All genes with both extremely high diversity (top fifth percentile) and homoplastic mutations were associated with drug resistance (i.e., gid, ethA, pncA, and rpsL). We identified 67 genes with high diversity and Tajima's D values more negative than expected based on gene length; only two of these were associated with drug resistance (i.e., ethA and pncA). Twenty out of 51 homoplastic SNPs that are also F ST outliers fall within or upstream of known drug resistance-associated genes. The remaining SNPs may be false-positive results or novel drug resistance mutations.

DISCUSSION
Highly virulent bacterial pathogens such as M. tuberculosis, Yersinia pestis (29), Francisella tularensis (30), and Mycobacterium ulcerans (31) appear to evolve clonally, i.e., with little to no evidence of lateral gene transfer. It is important to identify advantageous mutations in these and other organisms, as they are likely to be associated with phenotypes such as drug resistance, heightened transmissibility, or host adaptation. However, few methods are available for identifying loci under positive selection in the setting of clonal evolution. We adopted an empirical approach to this problem and used natural population data to characterize patterns of diversity at loci known to be under positive selection in M. tuberculosis.
In this analysis of clinical isolates from settings where drug resistance is endemic, we found genome-wide diversity to be higher in susceptible M. tuberculosis subpopulations than in the subpopulations resistant to first-and second-line drugs (with the exception of populations resistant to prothionamide [PRO] and moxifloxacin [MOX]). The observation of higher diversity in drug-susceptible populations is consistent with in South Africa concluded that XDR cases result primarily from transmission of resistance, rather than de novo evolution of resistance mutations during infection (9). The primary studies for the sequence data analyzed here also identified clusters of drugresistant isolates (5,6), suggesting that resistant isolates were being transmitted. Our results, along with these previously published observations, suggest that the fitness of drug-resistant isolates can be high enough to allow them to circulate in regions where drug-resistant M. tuberculosis is endemic. As discussed below, the fitness effects of M. tuberculosis drug resistance mutations appear to vary substantially; the finding of transmitted resistance in this and other studies suggests that the fitness of isolates harboring low-cost mutations is comparable to that of susceptible M. tuberculosis. The populations in our study have a high burden of drug-resistant TB, and the role of transmitted drug resistance may differ in other settings. An alternative-but not mutually exclusive-explanation for the observation of higher diversity in susceptible populations is that drug-resistant M. tuberculosis is under distinct evolutionary constraints that reduce average genome-wide levels of diversity. In support of this hypothesis, we identified a specific subset of genes that were invariant across drug-resistant populations. Interestingly, while average diversity was lower for resistant subpopulations, the gene-wise diversity distributions had heavier tails, indicating there were more genes with extreme levels of diversity.
We found the genetic architecture of resistance to vary among targets, and resistance-associated genes tended to fall within categories that we term "sloppy," "tight," and "hybrid" targets of selection (the latter has a combination of tight and sloppy features and applies to rpsL, embB, and rpoB). "Sloppy" resistance genes are characterized by high levels of diversity. Genes associated with PZA, EMB, ETO, and STR resistance (i.e., pncA, gid, rpsL, rrs, and ethA) have high levels of diversity; some also had an excess of rare variants (pncA, ethA, and embC). The finding that these genes accumulate multiple, individually rare mutations implies that there is a large target for resistance and/or compensatory mutations within the gene: that is, resistance can result from multiple different variants acting individually or in concert.
In addition to its numerous rare mutations, pncA also contains the highest number of homoplastic SNPs (27 SNPs emerged more than once in the phylogeny) of any gene in the data set. Among the 62 nonsynonymous pncA mutations in our data set, 55 have been previously reported in association with drug resistance (TB Drug Resistance Mutation Database [32]). The newly described SNPs may mediate drug resistance or compensation for the fitness effects of other variants. Relaxed purifying selection is likely to play a role in concert with selection for diverse advantageous resistance mutations in the accumulation of diversity in pncA and other sloppy targets. The fact that numerous mutations are segregating in a natural population suggests that alterations to these genes are generally associated with negligible fitness costs. An M. tuberculosis strain harboring a deletion in pncA conferring resistance to PZA was estimated to be endemic in Quebec, Canada, by 1800, long before the use of PZA for the treatment of TB (33)(34)(35). This supports the idea that purifying selection on pncA is relatively weak, which would contribute to its exceedingly high diversity and broaden the adaptive paths to resistance.
In contrast to pncA, gid, which is associated with low-level STR resistance (36), does not appear to have the signatures of a "sloppy" target for resistance despite its high diversity. We identified just three homoplastic SNPs within gid, and previous studies have found that STR-resistant isolates do not carry the same gid mutations (37). This could indicate that a multitude of mutations within gid confer resistance, but levels of diversity in the gene were similar in resistant and susceptible isolates. Previous studies of sequence polymorphism in gid have identified high diversity in this gene in both resistant and susceptible isolates (37)(38)(39): gid appears to be subject to relaxed purifying selection in the presence and absence of antibiotic pressure. Since gid mutations confer low-level resistance, it is also possible that misclassification of resistance phenotypes contributed to the lack of differentiation we and others have observed between putatively STR-resistant and -susceptible subpopulations. In addition, mutations in rpsL, which cause high-level resistance, could mask the contribution of gid to STR resistance.
We found some drug targets to be highly diverse in resistant subpopulations of either L2 or L4 (but not both), suggesting that resistance mutations in these genes interact with the genetic background; the fitness effects of mutations in these genes could, for example, vary on different genetic backgrounds. Lineage-specific F ST outliers are another category of candidate locus with lineage-dependent roles in drug resistance (Table 3). Epistatic interactions between drug resistance mutations and M. tuberculosis lineage have been reported previously: for example, specific mutations in the inhA promoter have been associated with the L1 and M. africanum genetic backgrounds (40,41).
In contrast to "sloppy" targets, we discovered individual homoplastic SNPs associated with drug-resistant subpopulations (i.e., with high F ST ) representing "tight" targets of selection in genes conferring resistance to INH, RIF, and STR. Numerous resistance mutations have been reported for katG, rpoB, rpsL, embB, and gyrA, but we find drug-resistant subpopulations to be defined by a specific subset of mutations in these genes. This suggests that certain mutations are strongly favored relative to others conferring resistance to the same drugs when M. tuberculosis is in its natural environment. Antibiotic resistance can impose fitness costs on M. tuberculosis during in vitro growth, with the range of fitness costs differing among mutations, even within the same gene (18). Mutations can also have different fitness effects depending on the genetic background, but the most fit mutants were the same across M. tuberculosis lineages in a study of RIF resistance (18).
In our analyses, we found the dominant INH resistance mutation in katG to affect the serine at position 315. This change reduces affinity to INH but preserves catalase activity (42) and is associated with lower fitness costs than other katG mutations, both in vitro and in a mouse model (43,44). This mutation was recently shown to precede mutations conferring resistance to other drugs during accumulation of resistance in the evolution of multidrug-resistant M. tuberculosis (45). The dominant mutations we identified in rpoB (codon 450) and rpsL (codon 43) have also been found to have lower fitness costs in vitro compared to other mutations conferring resistance to RIF and STR in these genes (18,44,46). These results suggest that many of the findings regarding the relative fitness costs of M. tuberculosis resistance mutations in vitro and in animal models are relevant to the pathogen's natural environment.
The fitness effects of mutations in gyrA (codon 94) and embB (codon 306) have not been measured; on the basis of our homoplasy and F ST results, we propose that they have lower fitness costs than other mutations in these genes and that they represent "tight" targets of selection. Mutations at gyrA codon 94 were previously found to be the most prevalent in a survey of gyrA and gyrB mutations in fluoroquinolone-resistant M. tuberculosis clinical isolates (47). Interestingly, the mutation in embB codon 306 has been previously associated with acquisition of multiple resistances (48), and we find that this position is an F ST outlier for all first-line drugs in L4. This mutation is not an F ST outlier in L2 (i.e., top fifth percentile), with the percentiles for F ST values ranging from 0.07 to 0.68 for first-line drugs in this lineage. These observations suggest that the genetic background affects interactions among resistance mutations and that embB codon 306 is important for acquisition of multidrug resistance in L4, but not in L2.
We searched for indels with the signature of a "tight" target, i.e., homoplastic mutations segregating at markedly different frequencies in drug-susceptible and -resistant subpopulations. Unlike the pattern observed with SNPs, genes associated with drug resistance were not significantly enriched among the genes harboring homoplastic indels. We identified one homoplastic indel that was also an F ST outlier-a deletion in gid that causes a frameshift. Patterns of variation in gid are complex and suggest a role for relaxation of purifying selection (i.e., in the accumulation of excess SNPs in both resistant and susceptible isolates) and perhaps a tight target associated with multiresistance (i.e., this homoplastic/F ST outlier deletion that was associated with resistance to CAP, EMB, ETO, KAN, OFL, and PZA).
Our finding that, save for the frameshift mutation in gid, indels in resistance genes do not have the signature of "tight" targets suggests that they are generally associated with higher fitness costs than SNPs. Fifteen drug targets have been found in transposon mutagenesis experiments to be essential for M. tuberculosis growth in vitro, including rpoB and rpsL; deletions in these genes are likely to interrupt important functions (49). Deletions in nonessential genes could also have fitness costs. Deletions in katG, which is nonessential, can result in INH resistance, but they are not observed as frequently in clinical isolates as the KatG S315 SNP is, particularly among transmitted INH-resistant strains (23).
There are several limitations to our study. Resistance to multiple drugs was common in our sample, and in some cases, it was difficult to identify patterns of diversity and population differentiation that were specific to individual drugs. Our results are also limited by the accuracy with which drug resistance phenotypes were determined and a lack of phenotypic data for some drugs (particularly second-line drugs). Our sample was heavily skewed to lineages 2 and 4, and the results are not necessarily applicable to other M. tuberculosis lineages. Finally, the data analyzed here were generated with short sequencing read technologies, and we were thus limited to characterizing diversity in regions of the M. tuberculosis genome that can be resolved with these methods: regions that were masked from analysis (e.g., due to sequence repeats) may include unknown resistance targets. We also used an L4 genome (H37Rv) as a reference, and gene content specific to L2 may not have been identified.
We were not able to recover all drug resistance-associated genes using the analyses performed here. This is likely a result of limited phenotypic data for some drugs and their associated targets (e.g., thyA and folC, which are associated with aminosalicylic acid [PAS] resistance). Our list of drug targets was dominated by genes associated with INH resistance, and signatures in genes that harbor rare resistance-associated alleles may be subtle compared to the KatG S315 mutation found at high frequency in drug-resistant populations.
We identified 31 SNPs that are not on the list of known drug resistance genes, which both emerged more than once in the phylogeny (homoplasies) and were segregating at markedly different frequencies in resistant and susceptible subpopulations (F ST outliers). These SNPs may be novel resistance determinants; notably, all nonsynonymous SNPs within this group are in genes linked with drug resistance in other studies (i.e., they are in genes encoding efflux pumps, genes differentially regulated in resistant isolates or in response to the presence of drug, potential drug targets, or genes in the same pathways as drug targets or resistance determinants) (50)(51)(52)(53)(54). In addition to a direct, previously unrecognized role in resistance, these SNPs could compensate for fitness costs of drug resistance. For example, we identified a homoplastic F ST outlier in rpoC, and mutations in rpoC have been shown to compensate for RIF resistance in experimental evolution studies (22). Intriguingly, we found lipid metabolism genes to be enriched among genes harboring homoplastic SNPs (P ϭ 0.013). We have previously shown that these genes have extreme values of diversity in a global sample of M. tuberculosis isolates and within individual hosts (26), suggesting that lipid metabolism genes may also be under positive selection in M. tuberculosis populations. The results presented here could be extended by phenotypic characterization of lipid profiles and identification of homoplastic variants that are at markedly different frequencies in isolates with distinct lipid profiles.
Here we have used drug resistance loci in M. tuberculosis to identify the signatures of positive selection in a clonal bacterium. We found these loci to be associated with distinct patterns of diversity that likely reflect differing genetic architectures underlying the traits under selection. The evolutionary path to resistance is broad for some drugs with "sloppy targets," whereas for drugs with "tight targets," the means of acquiring resistance appear more limited. This is likely due to fitness effects of resistance mutations in M. tuberculosis's natural environment, as numerous resistance mutations have been identified in tight target genes. We also found evidence suggesting that there are important interactions among loci during the evolution of resistance. Our results suggest that purifying selection on a subset of genes intensifies in the setting of resistance, which could reflect epistatic interactions and/or a response to the metabolic milieu imposed by antimycobacterial agents. The results presented here can be used to create more realistic models of resistance evolution in M. tuberculosis and to develop novel strategies of preventing or mitigating the acquisition of resistance. For example, the narrow path to resistance for drugs with tight targets reveals potentially exploitable vulnerabilities, as does the finding of interdependencies among specific loci and the genetic background in the evolution of resistance and multiresistance. As new TB drugs become available for clinical use, the approach outlined here can be extended to investigate their architectures of resistance.
Efforts are under way to sequence and perform drug susceptibility testing on thousands of M. tuberculosis isolates with the goal of creating an exhaustive catalogue of drug resistance mutations and eventually using whole-genome sequencing (WGS) to diagnose drug resistance in clinical settings (CRyPTIC project, http://modmedmicro .nsms.ox.ac.uk/cryptic/, last accessed 24 May 2017). We found that loci under positive selection can be identified using relatively simple methods: "tight" targets are highly differentiated in their allele frequencies across phenotypic groups (i.e., F ST outliers) and appear as homoplasies on the phylogeny; "sloppy" targets are characterized by high diversity and/or low Tajima's D values, as well as homoplasies. Extrapolating from patterns observed among known resistance variants, we have discovered new candidate regulatory and genic resistance variants. The methods used in this study are widely available and should scale to analysis of the large collections of genomic and phenotypic data that are currently being generated. This approach can be extended to identify novel resistance loci in bacteria for which drug susceptibility phenotypes are defined, as well as other positively selected loci in clonal bacterial populations.

MATERIALS AND METHODS
Reference guided assembly. We downloaded sequencing read data from two large surveys of drug-resistant M. tuberculosis in Russia (5) and South Africa (6). We used FastQC (55) and TrimGalore! (56) for quality assessment and adaptor trimming of the reads. Trimmed reads were mapped to M. tuberculosis H37Rv (accession no. NC_000962.3) using BWA-MEM v 0.7.12 (57). We used Samtools v 1.2 (58) and Picard Tools (https://broadinstitute.github.io/picard/) for sorting, format conversion, and addition of read group information. Variants were identified using Pilon v 1.16 (59). A detailed description of the reference guided assembly pipeline is available at https://github.com/pepperell-lab/RGAPepPipe. We removed isolates with mean coverage less than 20ϫ, isolates with percentage of the genome covered at 10ϫ less than 90%, isolates where a majority of reads did not map to H37Rv, and isolates where greater than 10% of sites were unknown after mapping. The final data set contains 1,161 M. tuberculosis isolates (see Table S1 in the supplemental material). The alignment was masked to remove repetitive regions, including proline-glutamate (PE)/proline-proline-glutamate (PPE) genes.
Phylogenetic analysis. We estimated the approximately maximum likelihood phylogeny using the masked alignment from reference guided assembly with FastTree 2.1.9 (60). We compiled FastTree using the double precision option to accurately estimate branch lengths of closely related isolates. We used FigTree (http://tree.bio.ed.ac.uk/software/figtree/) for tree visualization. SNP annotation. A VCF (variant call format) file of single nucleotide variants was created from the masked alignment using SNP-sites v 2.3.2 (61). Single nucleotide polymorphisms (SNPs) were annotated Selection at Drug Resistance Loci in M. tuberculosis January/February 2018 Volume 3 Issue 1 e00108-17 msystems.asm.org 13 using SnpEff v 4.1j (62) to identify synonymous, nonsynonymous, and intergenic SNPs based on the annotation of M. tuberculosis H37Rv. Indel identification. Insertions and deletions were identified during variant calling with Pilon. We used Emu (63) to normalize indels across multiple isolates. We used a presence/absence matrix for the normalized indels for further analyses of indel diversity.
Population genetics statistics. Whole-genome and gene-wise diversity ( and ) and neutrality (Tajima's D) statistics were calculated using Egglib v 2.1.10 (64) for whole-genome alignments and gene-wise alignments. Isolates were further divided by lineage and drug resistance phenotype. Sites with missing data due to indels or low-quality base calls in more than 5% of isolates in the alignment were not included in calculation of statistics. Tajima's D values showed a correlation with gene length in our sample. To find genes with extreme values of Tajima's D, we performed linear regression in R (65) on log-transformed Tajima's D values and gene length and identified genes with large residual values. To identify alleles with marked differences in frequency in resistant and susceptible isolates, Weir and Cockerham's F ST (66) was calculated using populations of resistant and susceptible isolates for each drug using vcflib v1.0.0-rc0-262-g50a3 (https://github.com/vcflib/vcflib). For nonbiallelic SNPs, we calculated F ST for the two most common variants.
Homoplasy. We used TreeTime (67) to perform ancestral reconstruction and place SNPs and indels in the phylogeny. We identified homoplastic SNPs and indels as those arising multiple times in the phylogeny.
Data availability. Unless otherwise noted, all data and scripts associated with this study are available at https://github.com/pepperell-lab/mtbDrugResistance.