BACKGROUND
tuberculosis  epidemiology in new york city has undergone dramatic changes since the resurgent tb epidemic of the 1990s, when over  <dig> cases were reported each year between  <dig> and  <dig>  many in outbreak clusters among vulnerable populations  <cit> .
tb incidence is now at an all-time low   <cit> , outbreak clusters have become increasingly rare, and so-called endemic clones have become a major source of new tb infections in the us-born population  <cit> .
genotyping of mycobacterium tuberculosis  clinical isolates is a cornerstone of tb control in new york city.
however, conventional genotyping methods  typing), interrogate less than  <dig> % of the approximately 4 mb m. tuberculosis genome and thus lack the discriminatory power to detect small-scale genetic differences within closely related populations.
in these situations, genotyping will often yield little if any useful information, even in isolates with wide geographic distribution and long epidemiological histories in a given population  <cit> .
whole genome sequencing  directly overcomes these limitations and has rapidly become an important, if not central, research tool in tb epidemiology: wgs-based studies have detected previously unknown outbreak clusters among isolates with identical miru-vntr types  <cit>  and identified so-called super-spreaders responsible for multiple secondary infections in the community  <cit> .
in addition, an expanding body of work has employed wgs data to address a wide-reaching set of previously uninvestigated questions in m. tuberculosis evolution and population genomics .
next-generation wgs technologies have markedly decreased per-isolate sequencing costs, and are expected to replace many current modalities in public health microbiology  <cit> .
specific applications of interest for tb control include rapid drug resistance typing, locating cryptic outbreak clusters and transmission hotspots not identified via field epidemiology, and identification and tracking of novel m. tuberculosis strains in the community.
snp-distance based strategies have proven useful for identifying recent tb transmission  <cit>  and wgs data has allowed for unprecedented phylogenetic resolution between and within m. tuberculosis subpopulations.
population genomics studies in both m. tuberculosis and other pathogens have established important linkages between the evolutionary and epidemiological histories of endemic and/or emerging pathogen subpopulations  <cit> .
specifically, emerging m. tuberculosis subpopulations are expected to exhibit low sequence diversity, an excess number of high frequency derived alleles, and potentially harbor strain-specific patterns of positive or purifying genomic selection.
multiple m. tuberculosis strains have emerged from new york city and neighboring new jersey  over the last two decades.
for example, m. tuberculosis isolates from the s <dig> group, a low-is <dig> copy number strain first identified in new jersey, usa  <cit>  in  <dig>  circulate within the nyc-nj area, predominantly among hiv-positive and homeless populations  <cit> .
the drug-susceptible c strain was first reported in nyc, where it has caused outbreaks among at-risk populations and sporadic cases in the general population, and then spread widely across the united states  <cit> .
both c and s <dig> strains belong to m. tuberculosis lineage  <dig>  the most widely distributed and successful of the six m. tuberculosis global phylogeographic lineages  <cit>  and the most prevalent lineage in the new york city area.
this study uses wgs data from tb isolates collected in new york city and new jersey between  <dig> and  <dig>  applying both phylogenomic and population genomics methods to identify epidemiologically-relevant subpopulations within this relatively monomorphic local population.
these methods identify previously known subpopulations  retrospectively, suggest useful measures for prospective and real-time identification of newly emerging isolate groups, and yield additional information on adaptation and epidemiological success in m. tuberculosis isolates endemic to new york city.
methods
mycobacterium tuberculosis isolates
seventy one total m. tuberculosis full genome sequences were included in this study.
fourty-seven isolates from lineage  <dig> were included:  <dig> isolates from tb cases occurring in new york city and new jersey between  <dig> and  <dig>  including 9 s <dig> isolates;  <dig> additional clinical isolates from sub-saharan africa  <cit>  and north america  <cit> ; and  <dig> well-characterized laboratory strains .
two additional isolates from new york city, from lineage  <dig> and lineage  <dig>  were also sequenced for this study.
sequence data for  <dig> additional non-l <dig> isolates, plus  <dig> isolates from the m. africanum-like lineage  <dig> and the outgroup m. bovis, were obtained from publicly available sources .
sequencing, alignment, and snp calling
wgs data were obtained for  <dig> previously-unsequenced m. tuberculosis clinical isolates .
isolates were cultured on löwenstein-jensen slants and grown at 37 °c for 3– <dig> weeks.
sequencing libraries were prepared using truseq dna or nextera dna preparation kits .
raw sequencing reads were generated on the illumina hiseq  <dig> platform and aligned to the h37rv reference genome  using the burrows-wheeler aligner  <cit> .
genome assemblies for all isolates were deposited in the ncbi genbank database .
all isolates had reads covering >99% of the reference genome, and the lowest mean coverage depth for any isolate was 27×.
snps were called using a phred-scaled quality threshold of  <dig>  and annotated using snpeff v <dig>  <cit> .
we excluded from analysis all variants occuring within pe and ppe genes, a family of highly repetitive, gc-rich m.tuberculosis genes in which recombination has been observed  <cit> .table  <dig> characteristics of the isolates sequenced in this study
availability of data and materials
the dataset supporting the conclusions of this article is available in the ncbi genbank repository  and supporting sequence alignments and phylogenetic tree data are available on treebase.
phylogenetic reconstruction
phylogenetic trees were estimated using maximum likelihood methods in the posix-threads build of raxml v <dig>  <cit> .
node robustness was assessed with  <dig> bootstrap pseudoreplicates and a consensus network was calculated  <cit>  as implemented in splitstree v <dig> .
<dig>  <cit> .
a custom perl script was used to identify snps with alleles unique to a given lineage or subpopulation.
neutrality statistics and selection analysis
neutrality statistics  were calculated in dnasp v <dig> .
<dig>  <cit>  with statistical significance assessed with  <dig> – <dig>  coalescent simulations.
fay and wu’s h is particularly useful for distinguishing whether a given departure from neutrality is attributable to recent population expansion or a selective sweep  <cit> .
the gene-wise ratio of the nonsynonymous substitution rate to the synonymous substitution rate  was estimated for every gene in the m. tuberculosis genome across all phylogenetic branches using the branch-site random effects likelihood  model as implemented in hyphy v <dig>  <cit> .
this model tests for branch-specific instances of episodic diversifying selection on every internal and terminal branch on the phylogenetic tree  and, following a likelihood-ratio test and holm’s correction for multiple tests, detects branches on which a proportion of the codons have evolved under a dn/ds ratio that is significantly different from that of the rest of the codons.
the advantage of this model over other so-called branch-site models is that it does not constrain the tree on either sides of the branch being tested to be subject to diversifying selection  and purifying selection .
RESULTS
population structure and genetic diversity
maximum likelihood phylogenomic reconstruction based on  <dig>  quality-filtered snps recovered primary phylogeographic lineages 1– <dig> and identified at least six distinct subpopulations within l <dig> isolates, including s <dig>  nucleotide diversity  ranged from  <dig> e- <dig> to  <dig> e- <dig> , consistent with prior estimates of genetic diversity within coding regions of the m. tuberculosis genome  <cit> .
s <dig> strains were separated from any other l <dig> isolate by at least  <dig> snps and exhibited lower nucleotide diversity and lower mean pairwise snp distances between isolates .table  <dig> genetic diversity and neutrality test statistics by lineage
d
t
d
fl
f
fl
r
2
h
fw
hn
fw
n, number of ingroup sequences; s, number of segregating sites; π, nucleotide diversity; k, average number of nucleotide differences; d
t, tajima’s d; r
2
, ramos-onsins and rozas’ r
<dig>  d
fl and f
fl, fu and li’s d and f ; h
fw, fay and wu’s h; hn
fw, fay and wu’s normalized h. statistical significance was assess with  <dig>  coalescent simulations .
*p <  <dig> , **p <  <dig>
drug resistance-associated polymorphisms
polymorphisms at drug resistance-associated codon sites were evaluated for  <dig> known drug resistance genes .
mutations in katg, which confer resistance to isoniazid, were common among isolates from lineages 1– <dig>   <dig>  and  <dig> and l <dig> isolates from western and sub-saharan africa, and rare among l <dig> isolates from n. america and europe, occurring in only a single isolate from this group.
s <dig> isolates were found to have a strain-specific mutation in emba  previously associated with ethambutol resistance  <cit> , however the s <dig> isolates included in this study are ethambutol-sensitive.
l <dig> isolates from kwazulu-natal, south africa carried drug resistance-associated mutations in katg, rpob, pnca, and rrs, consistent with prior studies on these drug-resistant isolates  <cit> .
subgroup-defining polymorphisms
one hundred seventeen synapomorphic snps  differentiate l <dig> isolates from the non-l <dig> isolates included in this study .
seventy-five additional snps differentiate north american isolates  from non-north american isolates, and  <dig> snps differentiate s <dig> from other north american isolates.
synapomorphic snps are unequally distributed by functional category, predominantly occurring in genes associated with cell wall functions, lipid metabolism, respiration, and intermediary metabolism.
non-synonymous synapomorphic snps occur in multiple genes with known or proposed functions in virulence, growth, and/or adaptation, including known virulence factors , and also components of gene-regulatory , lipid metabolism , intermediary metabolism , and cell-wall associated pathways  with known or proposed functions in m. tuberculosis virulence .
new york city and s <dig> isolates carry a unique non-synonymous mutation in rv1290c, a conserved gene of unknown function that when disrupted causes a severe attenuation of virulence  <cit> .
additional file 3: table s <dig> lists the complete set of subgroup-defining snps.fig.
<dig> synapomorphic polymorphisms by functional category and isolate subgroup.
a virulence and adaptation, b regulatory and information pathways, c conserved proteins without known function, d cell wall and lipid metabolism, e intermediary metabolism and respiration.
l <dig> includes all  lineage  <dig> isolates included in this study, nyc-nj  includes l <dig> isolates collected in new york city or new jersey, usa, including the s <dig> outbreak cluster, and s <dig>  includes isolates belonging the new jersey outbreak cluster described in the text.
genes carrying diagnostic snps with known functions in virulence, growth, and/or adaptation are listed above each column, and of these genes, those with non-synonymous polymorphisms are highlighted in yellow.
the number of total diagnostic snps unique to s <dig>  are listed in the third column
fig.
<dig>
a consensus network of  <dig> maximum likelihood bootstrap replicates for mycobacterium tuberculosis isolates from north america, sub-saharan africa, and asia  based on  <dig>  snps.
branches are color-coded by lineage.
isolates from the s <dig> cluster, identified in new jersey in 1997– <dig>  are highlighted; b world map of isolate collection locations color-coded by lineage
neutrality test statistics and population size expansion
site frequency-based neutrality test statistics were calculated using whole-genome polymorphism data by lineage  and by subgroups within l <dig>  including the s <dig> outbreak cluster and non-s <dig> isolates from new york city and new jersey  tajima’s d  and fu & li’s d and f test statistics , were significantly negative when calculated for all l <dig> isolates as a group  and for the subgroup of non-s <dig> isolates from new york city.
negative values for d
t
, d
fl, and f
fl indicate a relative excess of low frequency alleles in a population, which can occur following recent population expansion or a selective sweep  <cit> .
fay and wu’s h, a statistic that is insensitive to population expansion but highly sensitive to selection pressure, was not significantly different from zero for all isolate subgroups, suggesting that population expansion–rather than a selective sweep–explains the relative excess of rare alleles in isolates in l1- <dig> and non-s75 l <dig> isolates.
significant values for ramos-onsins & rozas’ r
<dig> statistic, which tests for recent population size expansion based on the difference between the number of singleton mutations and the mean number of nucleotide differences between samples, were observed for all subgroups except s <dig>  all five neutrality test statistics were non-significant for the s <dig> outbreak cluster.
unlike other subgroups, the unfolded site frequency spectrum for s <dig> exhibited a lower number of low-frequency alleles  and negative values for fay and wu’s h and normalized h, consistent with a small but non-significant excess of high-frequency derived alleles in this subpopulation.fig.
<dig> unfolded site frequency spectra for isolates from the s <dig> outbreak cluster  and non-s <dig> isolates from the new york city area.
dark and light blue bars indicate the number of non-synonymous and synonymous snps at each snp allele frequency
purifying selection on genes involved in lipid metabolism and cell wall maintenance
dn/ds was significantly less than  <dig>  consistent with purifying selection, for two genes in lipid metabolism pathways and five putative transmembrane protein genes .
two lipid metabolism pathway components, phenolphthiocerol synthesis polyketide synthase a–e family  and polyketide synthase  <dig> , exhibited significantly decreased dn/ds in a specific subpopulation of new york city isolates .
evidence of episodic diversifying selection, with dn/ds significantly greater than  <dig>  was limited to three isolates in our study, the l <dig> isolate w <dig>  the l <dig> isolate t <dig>  and the m. africanum k <dig> isolate.
discussion
m. tuberculosis exhibits very low sequence diversity compared to other bacteria, minimal evidence of horizontal gene transfer , and recombination limited to known highly variable gene families  <cit> .
this lack of genetic diversity is pronounced in geographically restricted m. tuberculosis populations, such that locally endemic clone groups have posed a unique challenge to laboratory-based identification of tb outbreak clusters in new york city.
historically, these isolates have been strongly associated with homeless and at-risk populations, in which field epidemiology and contact tracing are often difficult, placing a premium on rapid and reliable laboratory identification of clustered cases.
in one case series, 52% of patients infected with c strain isolates in nyc had no phone number or address, or could otherwise not be contacted by public health investigators  <cit> .
the present study demonstrates how whole genome-based laboratory analysis can overcome these challenges, and suggests that wgs may be a particularly important tool at the local level, where genetic diversity is expected to be lower compared to more geographically diverse samples.
the results presented here provide three specific approaches for identifying outbreak clusters and emerging strain groups in local m. tuberculosis populations:  genome-based phylogenetic reconstruction;  population genetic analysis, specifically estimation of neutrality and diversity statistics within grouped samples; and  genome-wide analysis for distinguishing signatures of purifying or diversifying selection.
whole genome-based phylogenetic reconstruction yielded clearly defined population substructure among locally-endemic isolates in the new york city area, and identified s <dig> isolates as distinct clade in the phylogeny.
s <dig> isolates also exhibited poorly differentiated terminal branching patterns, and relatively lower bootstrap support at individual nodes, which may reflect the limits of phylogenetic resolution inherent to available genome sequencing data.
although this approach allows for robust retrospective identification of outbreak clusters and emerging strain groups, it is perhaps less well suited for rapid identification of clustered transmission among new tb cases, in which low levels of genetic differentiation may preclude high-confidence phylogenomic resolution between isolates.
however, as wgs-based technologies replace conventional genotyping methods, phylogenetic reconstruction will likely become an important tool in public health microbiology, providing a “phylogenetic reference” for tb isolates sequenced within a given geographic area or tb control program  <cit> .
in addition, clustering of incident isolates in a specific phylogenetic branch could suggest ongoing transmission within a specific at-risk population.
snp distance-based inference of recent transmission, in which the pairwise snp distance is used to infer whether two isolates were transmitted directly between their respective hosts, is likely to become an important epidemiological tool in tb control  <cit> .
although the distribution of pairwise snp differences is expected to vary between low- and high-transmission areas   <cit> , emerging m. tuberculosis subpopulations are still expected to exhibit relatively few snp differences between isolates.
identification of subpopulation-defining synapomorphic polymorphisms can support this approach by identifying unique snps shared between isolates in emerging subpopulations.
the two additional methods used in this study  are likely to have more value in retrospective analyses, where they can yield useful information about the epidemiological and evolutionary history of circulating m. tuberculosis subpopulations.
subgroup-defining polymorphisms can provide useful genetic markers for m. tuberculosis strain identification, similar to other minimal snp sets used in m. tuberculosis phylogenetics  <cit> .
s <dig> isolates in this study could be distinguished from other north american isolates using only  <dig> snps, and determination of similar subgroup-defining snp sets could provide a straightforward tool for rapidly determining if a given tb isolate belongs to an existing outbreak cluster or endemic strain group.
more broadly, subgroup-defining polymorphisms also provide interesting, if limited, insight into the evolutionary history of lineage 4 m. tuberculosis isolates in north america and the specific l <dig> populations endemic to new york city.
isolates in these populations have only a minimal number of drug resistance-associated mutations, and instead have acquired multiple non-synonymous snps in virulence- and growth-associated pathways.
mutations in pks <dig> and yrbe2b are of particular interest, first because of their well-characterized roles in m. tuberculosis virulence and second because they may both influence tnf-mediated host immune responses  <cit> .
s <dig> isolates strains are known to induce higher levels of tnf-α in vitro  <cit> , which may help explain why s <dig> isolates have spread preferentially in immunocomproised patients.
although the functional consequences of these mutations are still unknown, these findings suggests that overall pathoadaptive fitness, rather than the acquisition of drug resistance mutations, may have played an important role in the evolutionary history of l <dig> m. tuberculosis populations in new york city.
selection analysis identified two loci in m. tuberculosis lipid metabolism pathways, ppsa and pks <dig>  with significantly decreased dn/ds ratios consistent with evolution under strong purifying selection.
observing signatures of purifying selection localized to a single subpopulation , may suggest adaptation to a particular subpopulation or transmission niche, and thus provide useful information about risk factors for acquisition of infection with an emerging strain group.
ppsa and other pps family genes are involved in the synthesis of phthiocerol and phenolphthiocerol, two components of cell wall lipids unique to pathogenic mycobacteria that likely participate in host-pathogen interactions  <cit>  and virulence  <cit> .
interestingly, farhat et al.
identified ppsa and pks <dig> among  <dig> genes that exhibit signatures of convergence and possible positive selection in multidrug-resistant m. tuberculosis isolates  <cit> .
although these loci may exhibit signatures of positive selection in drug-resistant populations, it is not unexpected that ppsa and pks <dig> would exhibit signatures of purifying selection in populations without a similar history of drug selection pressure.
consistent with this hypothesis we observed relatively fewer drug resistance-associated mutations in the same subpopulations where ppsa and pks <dig> exhibit signatures of purifying selection.
furthermore we found that dn/ds ratios at known drug-resistance loci were not significantly greater than one in our sample, consistent with prior studies in drug-susceptible m. tuberculosis isolates  <cit> .
the dn/ds ratio has limited power to detect positive selection in recently diverged intraspecific sequences and may underestimate the magnitude of negative selection in genes under strong purifying selection  <cit> .
however, because the dn/ds ratio is expected to underestimate the magnitude of the selection coefficient in this context, our analysis is likely conservative, and the true magnitude of negative selection on ppsa and pks <dig> may be larger than we have reported.
lastly, estimation of multiple neutrality statistics yielded evidence for past population expansion across multiple subpopulations, consistent with prior studies on demographic expansion in m. tuberculosis populations  <cit> , with the notable exception of s <dig>  this finding, in conjunction with the negative but nonsignificant h values estimated for s <dig> isolates , is consistent with the epidemiological history of this recently diverged group of closely related isolates.
however, it is important to acknowledge that factors such as sample size and time since demographic expansion can influence the power of statistics that draw from the site frequency spectrum to detect past population growth.
specifically, site frequency spectrum-based statistics may fail to detect population expansion if the elapsed time since an expansion is either too small or too large, or with small sample sizes  <cit> , and thus may be less useful for identification of emerging strain groups, as illustrated here.
importantly, the retrospective sample used in this study includes less than  <dig> % of all m. tuberculosis infections occurring in new york city between  <dig> and  <dig>  <cit> .
nevertheless, this study demonstrates that even a small sample of isolates can yield meaningful information about the epidemiological and evolutionary history of endemic m. tuberculosis isolate groups in low-transmission settings.
CONCLUSIONS
wgs-based technologies are likely to replace many conventional genotyping methods currently used in public health microbiology and tb epidemiology.
how to maximize the public health value of this paradigm shift, and the large quantities of genomic data it will soon make available, is still an open question.
whole genome-based drug resistance profiling, snp distance-based methods to identify ongoing transmission, and phylogenetic reconstruction will likely yield the most direct, practical benefits, and the wgs data collected during these activities will provide an important resource for ongoing research in tb epidemiology and pathogen evolution.
additional files
additional file 1: table s <dig>  complete list of non-synonymous and synonymous synapomorphic polymorphisms for all l <dig> isolates, l <dig> isolates from new york city and new jersey, and isolates from the s <dig> outbreak cluster.
additional file 2:
figure s <dig>  whole-genome maximum likelihood phylogenetic reconstruction of mycobacterium tuberculosis isolates from north america, sub-saharan africa, and asia .
values at the nodes indicate branch support based on  <dig> bootstrap replicates.
letter labels denote branches with genes under purifying selection  or lineage-defining polymorphisms.
green boxes in the adjoining matrix indicate snps at drug resistance-associated codon sites in known drug resistance gene.
inh: isonaizid; rif: rifampin; pza: pyrazinamide; eth: ethionamide; emb: ethambutol; ami: amikacin; flq: fluoroquinolones; pas: para-aminosalicylic acid; sm: streptomycin.
additional file 3:
table s <dig>  statistically significant dn/ds values, estimated by gene and across phylogenetic branches.
p
holm, holm-corrected p-value.
location of nodes a–e are indicated in fig.
<dig>
abbreviations
mirumycobacterial interspersed repetitive units
snpsingle nucleotide polymorphism
tbtuberculosis
wgswhole genome sequencing

 tuberculosis epidemiology in new york city has undergone dramatic changes since the resurgent tb epidemic of the 1990s, when over cases were reported each year between and many in outbreak clusters among vulnerable populations. tb incidence is now at an all-time low, outbreak clusters have become increasingly rare, and so-called endemic clones have become a major source of new tb infections in the us-born population. genotyping of mycobacterium tuberculosis clinical isolates is a cornerstone of tb control in new york city. however, conventional genotyping methods typing), interrogate less than % of the approximately 4 mb m. tuberculosis genome and thus lack the discriminatory power to detect small-scale genetic differences within closely related populations. in these situations, genotyping will often yield little if any useful information, even in isolates with wide geographic distribution and long epidemiological histories in a given population. whole genome sequencing directly overcomes these limitations and has rapidly become an important, if not central, research tool in tb epidemiology: wgs-based studies have detected previously unknown outbreak clusters among isolates with identical miru-vntr types and identified so-called super-spreaders responsible for multiple secondary infections in the community. in addition, an expanding body of work has employed wgs data to address a wide-reaching set of previously uninvestigated questions in m. tuberculosis evolution and population genomics. next-generation wgs technologies have markedly decreased per-isolate sequencing costs, and are expected to replace many current modalities in public health microbiology. specific applications of interest for tb control include rapid drug resistance typing, locating cryptic outbreak clusters and transmission hotspots not identified via field epidemiology, and identification and tracking of novel m. tuberculosis strains in the community. snp-distance based strategies have proven useful for identifying recent tb transmission and wgs data has allowed for unprecedented phylogenetic resolution between and within m. tuberculosis subpopulations. population genomics studies in both m. tuberculosis and other pathogens have established important linkages between the evolutionary and epidemiological histories of endemic and/or emerging pathogen subpopulations. specifically, emerging m. tuberculosis subpopulations are expected to exhibit low sequence diversity, an excess number of high frequency derived alleles, and potentially harbor strain-specific patterns of positive or purifying genomic selection. multiple m. tuberculosis strains have emerged from new york city and neighboring new jersey over the last two decades. for example, m. tuberculosis isolates from the s group, a low-is copy number strain first identified in new jersey, usa in circulate within the nyc-nj area, predominantly among hiv-positive and homeless populations. the drug-susceptible c strain was first reported in nyc, where it has caused outbreaks among at-risk populations and sporadic cases in the general population, and then spread widely across the united states. both c and s strains belong to m. tuberculosis lineage the most widely distributed and successful of the six m. tuberculosis global phylogeographic lineages and the most prevalent lineage in the new york city area. this study uses wgs data from tb isolates collected in new york city and new jersey between and applying both phylogenomic and population genomics methods to identify epidemiologically-relevant subpopulations within this relatively monomorphic local population. these methods identify previously known subpopulations retrospectively, suggest useful measures for prospective and real-time identification of newly emerging isolate groups, and yield additional information on adaptation and epidemiological success in m. tuberculosis isolates endemic to new york city. methods mycobacterium tuberculosis isolates seventy one total m. tuberculosis full genome sequences were included in this study. fourty-seven isolates from lineage were included: isolates from tb cases occurring in new york city and new jersey between and including 9 s isolates; additional clinical isolates from sub-saharan africa and north america; and well-characterized laboratory strains. two additional isolates from new york city, from lineage and lineage were also sequenced for this study. sequence data for additional non-l isolates, plus isolates from the m. africanum-like lineage and the outgroup m. bovis, were obtained from publicly available sources. sequencing, alignment, and snp calling wgs data were obtained for previously-unsequenced m. tuberculosis clinical isolates. isolates were cultured on löwenstein-jensen slants and grown at 37 °c for 3– weeks. sequencing libraries were prepared using truseq dna or nextera dna preparation kits. raw sequencing reads were generated on the illumina hiseq platform and aligned to the h37rv reference genome using the burrows-wheeler aligner. genome assemblies for all isolates were deposited in the ncbi genbank database. all isolates had reads covering >99% of the reference genome, and the lowest mean coverage depth for any isolate was 27×. snps were called using a phred-scaled quality threshold of and annotated using snpeff v. we excluded from analysis all variants occuring within pe and ppe genes, a family of highly repetitive, gc-rich m.tuberculosis genes in which recombination has been observed.table characteristics of the isolates sequenced in this study availability of data and materials the dataset supporting the conclusions of this article is available in the ncbi genbank repository and supporting sequence alignments and phylogenetic tree data are available on treebase. phylogenetic reconstruction phylogenetic trees were estimated using maximum likelihood methods in the posix-threads build of raxml v. node robustness was assessed with bootstrap pseudoreplicates and a consensus network was calculated as implemented in splitstree v.. a custom perl script was used to identify snps with alleles unique to a given lineage or subpopulation. neutrality statistics and selection analysis neutrality statistics were calculated in dnasp v. with statistical significance assessed with – coalescent simulations. fay and wu’s h is particularly useful for distinguishing whether a given departure from neutrality is attributable to recent population expansion or a selective sweep. the gene-wise ratio of the nonsynonymous substitution rate to the synonymous substitution rate was estimated for every gene in the m. tuberculosis genome across all phylogenetic branches using the branch-site random effects likelihood model as implemented in hyphy v. this model tests for branch-specific instances of episodic diversifying selection on every internal and terminal branch on the phylogenetic tree and, following a likelihood-ratio test and holm’s correction for multiple tests, detects branches on which a proportion of the codons have evolved under a dn/ds ratio that is significantly different from that of the rest of the codons. the advantage of this model over other so-called branch-site models is that it does not constrain the tree on either sides of the branch being tested to be subject to diversifying selection and purifying selection.  population structure and genetic diversity maximum likelihood phylogenomic reconstruction based on quality-filtered snps recovered primary phylogeographic lineages 1– and identified at least six distinct subpopulations within l isolates, including s nucleotide diversity ranged from e- to e-, consistent with prior estimates of genetic diversity within coding regions of the m. tuberculosis genome. s strains were separated from any other l isolate by at least snps and exhibited lower nucleotide diversity and lower mean pairwise snp distances between isolates.table genetic diversity and neutrality test statistics by lineage d t d fl f fl r 2 h fw hn fw n, number of ingroup sequences; s, number of segregating sites; π, nucleotide diversity; k, average number of nucleotide differences; d t, tajima’s d; r 2, ramos-onsins and rozas’ r d fl and f fl, fu and li’s d and f; h fw, fay and wu’s h; hn fw, fay and wu’s normalized h. statistical significance was assess with coalescent simulations. *p < , **p <  drug resistance-associated polymorphisms polymorphisms at drug resistance-associated codon sites were evaluated for known drug resistance genes. mutations in katg, which confer resistance to isoniazid, were common among isolates from lineages 1– and and l isolates from western and sub-saharan africa, and rare among l isolates from n. america and europe, occurring in only a single isolate from this group. s isolates were found to have a strain-specific mutation in emba previously associated with ethambutol resistance, however the s isolates included in this study are ethambutol-sensitive. l isolates from kwazulu-natal, south africa carried drug resistance-associated mutations in katg, rpob, pnca, and rrs, consistent with prior studies on these drug-resistant isolates. subgroup-defining polymorphisms one hundred seventeen synapomorphic snps differentiate l isolates from the non-l isolates included in this study. seventy-five additional snps differentiate north american isolates from non-north american isolates, and snps differentiate s from other north american isolates. synapomorphic snps are unequally distributed by functional category, predominantly occurring in genes associated with cell wall functions, lipid metabolism, respiration, and intermediary metabolism. non-synonymous synapomorphic snps occur in multiple genes with known or proposed functions in virulence, growth, and/or adaptation, including known virulence factors, and also components of gene-regulatory, lipid metabolism, intermediary metabolism, and cell-wall associated pathways with known or proposed functions in m. tuberculosis virulence. new york city and s isolates carry a unique non-synonymous mutation in rv1290c, a conserved gene of unknown function that when disrupted causes a severe attenuation of virulence. additional file 3: table s lists the complete set of subgroup-defining snps.fig. synapomorphic polymorphisms by functional category and isolate subgroup. a virulence and adaptation, b regulatory and information pathways, c conserved proteins without known function, d cell wall and lipid metabolism, e intermediary metabolism and respiration. l includes all lineage isolates included in this study, nyc-nj includes l isolates collected in new york city or new jersey, usa, including the s outbreak cluster, and s includes isolates belonging the new jersey outbreak cluster described in the text. genes carrying diagnostic snps with known functions in virulence, growth, and/or adaptation are listed above each column, and of these genes, those with non-synonymous polymorphisms are highlighted in yellow. the number of total diagnostic snps unique to s are listed in the third column fig. a consensus network of maximum likelihood bootstrap replicates for mycobacterium tuberculosis isolates from north america, sub-saharan africa, and asia based on snps. branches are color-coded by lineage. isolates from the s cluster, identified in new jersey in 1997– are highlighted; b world map of isolate collection locations color-coded by lineage neutrality test statistics and population size expansion site frequency-based neutrality test statistics were calculated using whole-genome polymorphism data by lineage and by subgroups within l including the s outbreak cluster and non-s isolates from new york city and new jersey tajima’s d and fu & li’s d and f test statistics, were significantly negative when calculated for all l isolates as a group and for the subgroup of non-s isolates from new york city. negative values for d t, d fl, and f fl indicate a relative excess of low frequency alleles in a population, which can occur following recent population expansion or a selective sweep. fay and wu’s h, a statistic that is insensitive to population expansion but highly sensitive to selection pressure, was not significantly different from zero for all isolate subgroups, suggesting that population expansion–rather than a selective sweep–explains the relative excess of rare alleles in isolates in l1- and non-s75 l isolates. significant values for ramos-onsins & rozas’ r statistic, which tests for recent population size expansion based on the difference between the number of singleton mutations and the mean number of nucleotide differences between samples, were observed for all subgroups except s all five neutrality test statistics were non-significant for the s outbreak cluster. unlike other subgroups, the unfolded site frequency spectrum for s exhibited a lower number of low-frequency alleles and negative values for fay and wu’s h and normalized h, consistent with a small but non-significant excess of high-frequency derived alleles in this subpopulation.fig. unfolded site frequency spectra for isolates from the s outbreak cluster and non-s isolates from the new york city area. dark and light blue bars indicate the number of non-synonymous and synonymous snps at each snp allele frequency purifying selection on genes involved in lipid metabolism and cell wall maintenance dn/ds was significantly less than consistent with purifying selection, for two genes in lipid metabolism pathways and five putative transmembrane protein genes. two lipid metabolism pathway components, phenolphthiocerol synthesis polyketide synthase a–e family and polyketide synthase, exhibited significantly decreased dn/ds in a specific subpopulation of new york city isolates. evidence of episodic diversifying selection, with dn/ds significantly greater than was limited to three isolates in our study, the l isolate w the l isolate t and the m. africanum k isolate. discussion m. tuberculosis exhibits very low sequence diversity compared to other bacteria, minimal evidence of horizontal gene transfer, and recombination limited to known highly variable gene families. this lack of genetic diversity is pronounced in geographically restricted m. tuberculosis populations, such that locally endemic clone groups have posed a unique challenge to laboratory-based identification of tb outbreak clusters in new york city. historically, these isolates have been strongly associated with homeless and at-risk populations, in which field epidemiology and contact tracing are often difficult, placing a premium on rapid and reliable laboratory identification of clustered cases. in one case series, 52% of patients infected with c strain isolates in nyc had no phone number or address, or could otherwise not be contacted by public health investigators. the present study demonstrates how whole genome-based laboratory analysis can overcome these challenges, and suggests that wgs may be a particularly important tool at the local level, where genetic diversity is expected to be lower compared to more geographically diverse samples. the results presented here provide three specific approaches for identifying outbreak clusters and emerging strain groups in local m. tuberculosis populations: genome-based phylogenetic reconstruction; population genetic analysis, specifically estimation of neutrality and diversity statistics within grouped samples; and genome-wide analysis for distinguishing signatures of purifying or diversifying selection. whole genome-based phylogenetic reconstruction yielded clearly defined population substructure among locally-endemic isolates in the new york city area, and identified s isolates as distinct clade in the phylogeny. s isolates also exhibited poorly differentiated terminal branching patterns, and relatively lower bootstrap support at individual nodes, which may reflect the limits of phylogenetic resolution inherent to available genome sequencing data. although this approach allows for robust retrospective identification of outbreak clusters and emerging strain groups, it is perhaps less well suited for rapid identification of clustered transmission among new tb cases, in which low levels of genetic differentiation may preclude high-confidence phylogenomic resolution between isolates. however, as wgs-based technologies replace conventional genotyping methods, phylogenetic reconstruction will likely become an important tool in public health microbiology, providing a “phylogenetic reference” for tb isolates sequenced within a given geographic area or tb control program. in addition, clustering of incident isolates in a specific phylogenetic branch could suggest ongoing transmission within a specific at-risk population. snp distance-based inference of recent transmission, in which the pairwise snp distance is used to infer whether two isolates were transmitted directly between their respective hosts, is likely to become an important epidemiological tool in tb control. although the distribution of pairwise snp differences is expected to vary between low- and high-transmission areas, emerging m. tuberculosis subpopulations are still expected to exhibit relatively few snp differences between isolates. identification of subpopulation-defining synapomorphic polymorphisms can support this approach by identifying unique snps shared between isolates in emerging subpopulations. the two additional methods used in this study are likely to have more value in retrospective analyses, where they can yield useful information about the epidemiological and evolutionary history of circulating m. tuberculosis subpopulations. subgroup-defining polymorphisms can provide useful genetic markers for m. tuberculosis strain identification, similar to other minimal snp sets used in m. tuberculosis phylogenetics. s isolates in this study could be distinguished from other north american isolates using only snps, and determination of similar subgroup-defining snp sets could provide a straightforward tool for rapidly determining if a given tb isolate belongs to an existing outbreak cluster or endemic strain group. more broadly, subgroup-defining polymorphisms also provide interesting, if limited, insight into the evolutionary history of lineage 4 m. tuberculosis isolates in north america and the specific l populations endemic to new york city. isolates in these populations have only a minimal number of drug resistance-associated mutations, and instead have acquired multiple non-synonymous snps in virulence- and growth-associated pathways. mutations in pks and yrbe2b are of particular interest, first because of their well-characterized roles in m. tuberculosis virulence and second because they may both influence tnf-mediated host immune responses. s isolates strains are known to induce higher levels of tnf-α in vitro, which may help explain why s isolates have spread preferentially in immunocomproised patients. although the functional consequences of these mutations are still unknown, these findings suggests that overall pathoadaptive fitness, rather than the acquisition of drug resistance mutations, may have played an important role in the evolutionary history of l m. tuberculosis populations in new york city. selection analysis identified two loci in m. tuberculosis lipid metabolism pathways, ppsa and pks with significantly decreased dn/ds ratios consistent with evolution under strong purifying selection. observing signatures of purifying selection localized to a single subpopulation, may suggest adaptation to a particular subpopulation or transmission niche, and thus provide useful information about risk factors for acquisition of infection with an emerging strain group. ppsa and other pps family genes are involved in the synthesis of phthiocerol and phenolphthiocerol, two components of cell wall lipids unique to pathogenic mycobacteria that likely participate in host-pathogen interactions and virulence. interestingly, farhat et al. identified ppsa and pks among genes that exhibit signatures of convergence and possible positive selection in multidrug-resistant m. tuberculosis isolates. although these loci may exhibit signatures of positive selection in drug-resistant populations, it is not unexpected that ppsa and pks would exhibit signatures of purifying selection in populations without a similar history of drug selection pressure. consistent with this hypothesis we observed relatively fewer drug resistance-associated mutations in the same subpopulations where ppsa and pks exhibit signatures of purifying selection. furthermore we found that dn/ds ratios at known drug-resistance loci were not significantly greater than one in our sample, consistent with prior studies in drug-susceptible m. tuberculosis isolates. the dn/ds ratio has limited power to detect positive selection in recently diverged intraspecific sequences and may underestimate the magnitude of negative selection in genes under strong purifying selection. however, because the dn/ds ratio is expected to underestimate the magnitude of the selection coefficient in this context, our analysis is likely conservative, and the true magnitude of negative selection on ppsa and pks may be larger than we have reported. lastly, estimation of multiple neutrality statistics yielded evidence for past population expansion across multiple subpopulations, consistent with prior studies on demographic expansion in m. tuberculosis populations, with the notable exception of s this finding, in conjunction with the negative but nonsignificant h values estimated for s isolates, is consistent with the epidemiological history of this recently diverged group of closely related isolates. however, it is important to acknowledge that factors such as sample size and time since demographic expansion can influence the power of statistics that draw from the site frequency spectrum to detect past population growth. specifically, site frequency spectrum-based statistics may fail to detect population expansion if the elapsed time since an expansion is either too small or too large, or with small sample sizes, and thus may be less useful for identification of emerging strain groups, as illustrated here. importantly, the retrospective sample used in this study includes less than % of all m. tuberculosis infections occurring in new york city between and. nevertheless, this study demonstrates that even a small sample of isolates can yield meaningful information about the epidemiological and evolutionary history of endemic m. tuberculosis isolate groups in low-transmission settings.  wgs-based technologies are likely to replace many conventional genotyping methods currently used in public health microbiology and tb epidemiology. how to maximize the public health value of this paradigm shift, and the large quantities of genomic data it will soon make available, is still an open question. whole genome-based drug resistance profiling, snp distance-based methods to identify ongoing transmission, and phylogenetic reconstruction will likely yield the most direct, practical benefits, and the wgs data collected during these activities will provide an important resource for ongoing research in tb epidemiology and pathogen evolution. additional files additional file 1: table s complete list of non-synonymous and synonymous synapomorphic polymorphisms for all l isolates, l isolates from new york city and new jersey, and isolates from the s outbreak cluster. additional file 2: figure s whole-genome maximum likelihood phylogenetic reconstruction of mycobacterium tuberculosis isolates from north america, sub-saharan africa, and asia. values at the nodes indicate branch support based on bootstrap replicates. letter labels denote branches with genes under purifying selection or lineage-defining polymorphisms. green boxes in the adjoining matrix indicate snps at drug resistance-associated codon sites in known drug resistance gene. inh: isonaizid; rif: rifampin; pza: pyrazinamide; eth: ethionamide; emb: ethambutol; ami: amikacin; flq: fluoroquinolones; pas: para-aminosalicylic acid; sm: streptomycin. additional file 3: table s statistically significant dn/ds values, estimated by gene and across phylogenetic branches. p holm, holm-corrected p-value. location of nodes a–e are indicated in fig.  abbreviations mirumycobacterial interspersed repetitive units snpsingle nucleotide polymorphism tbtuberculosis wgswhole genome sequencing 
in addition, we observed that isolates from new york city endemic clone groups have acquired multiple non-synonymous snps in virulence- and growth-associated pathways, and relatively few mutations in drug resistance-associated genes, suggesting that overall pathoadaptive fitness, rather than the acquisition of drug resistance mutations, has played a central role in the evolutionary history and epidemiology of m. tuberculosis subpopulations in new york city. m. tuberculosis subpopulations in new york city and new jersey can be distinguished via phylogenomic reconstruction, evidence of demographic expansion and subpopulation-specific signatures of selection, and by determination of subgroup-defining nucleotide substitutions. wgs-based methods may be particularly useful in areas with less diverse mycobacterium tuberculosis populations, such as new york city, where conventional genotyping is often uninformative and field epidemiology often difficult. neutrality statistics based on the site frequency spectrum were less useful for identifying m. tuberculosis subpopulations, likely due to the low levels of informative genetic variation in recently diverged isolate groups. our results demonstrate that some but not all wgs-based methods are useful for detection of emerging m. tuberculosis clone groups, and support the use of phylogenomic reconstruction in routine tuberculosis laboratory surveillance, particularly in areas with relatively less diverse m. tuberculosis populations. whole genome sequencing has rapidly become an important research tool in tuberculosis epidemiology and is likely to replace many existing methods in public health microbiology in the near future. our study also supports the use of wider-reaching phylogenomic and population genomic methods in tuberculosis public health practice, which can support tuberculosis control activities by identifying genetic polymorphisms contributing to epidemiological success in local m. tuberculosis populations and possibly explain why certain isolate groups are apparently more successful in specific host populations. mycobacterium tuberculosiswhole genome sequencingphylogenomicssurveillanceissue-copyright-statement© the author 2016 this study applies four candidate strategies for wgs-based identification of emerging m. tuberculosis subpopulations, employing both phylogenomic and population genetics methods.
 whole genome sequencing has rapidly become an important research tool in tuberculosis epidemiology and is likely to replace many existing methods in public health microbiology in the near future. wgs-based methods may be particularly useful in areas with less diverse mycobacterium tuberculosis populations, such as new york city, where conventional genotyping is often uninformative and field epidemiology often difficult. this study applies four candidate strategies for wgs-based identification of emerging m. tuberculosis subpopulations, employing both phylogenomic and population genetics methods.  m. tuberculosis subpopulations in new york city and new jersey can be distinguished via phylogenomic reconstruction, evidence of demographic expansion and subpopulation-specific signatures of selection, and by determination of subgroup-defining nucleotide substitutions. these methods identified known historical outbreak clusters and previously unidentified subpopulations within relatively monomorphic m. tuberculosis endemic clone groups. neutrality statistics based on the site frequency spectrum were less useful for identifying m. tuberculosis subpopulations, likely due to the low levels of informative genetic variation in recently diverged isolate groups. in addition, we observed that isolates from new york city endemic clone groups have acquired multiple non-synonymous snps in virulence- and growth-associated pathways, and relatively few mutations in drug resistance-associated genes, suggesting that overall pathoadaptive fitness, rather than the acquisition of drug resistance mutations, has played a central role in the evolutionary history and epidemiology of m. tuberculosis subpopulations in new york city.  our results demonstrate that some but not all wgs-based methods are useful for detection of emerging m. tuberculosis clone groups, and support the use of phylogenomic reconstruction in routine tuberculosis laboratory surveillance, particularly in areas with relatively less diverse m. tuberculosis populations. our study also supports the use of wider-reaching phylogenomic and population genomic methods in tuberculosis public health practice, which can support tuberculosis control activities by identifying genetic polymorphisms contributing to epidemiological success in local m. tuberculosis populations and possibly explain why certain isolate groups are apparently more successful in specific host populations. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords mycobacterium tuberculosiswhole genome sequencingphylogenomicssurveillanceissue-copyright-statement© the author 2016 