High Functional Diversity in Mycobacterium tuberculosis Driven by Genetic Drift and Human ...

14 Pages·2008·2.55 MB·English

by  System Administrator

PLoS
BIOLOGY
High Functional Diversity in Mycobacterium
tuberculosis Driven by Genetic Drift and
Human Demography
Ruth Hershberg1[, Mikhail Lipatov1[, Peter M. Small2,3, Hadar Sheffer2, Stefan Niemann4, Susanne Homolka4,
Jared C. Roach5, Kristin Kremer6, Dmitri A. Petrov1, Marcus W. Feldman1, Sebastien Gagneux2,7*
1DepartmentofBiology,StanfordUniversity,Stanford,California,UnitedStatesofAmerica,2InstituteforSystemsBiology,Seattle,Washington,UnitedStatesofAmerica,3
BillandMelindaGatesFoundation,Seattle,Washington,UnitedStatesofAmerica,4ForschungszentrumBorstel,NationalReferenceCenterforMycobacteria,Borstel,
Germany,5SeattleChildren’sHospitalResearchInstitute,Seattle,Washington,UnitedStatesofAmerica,6MycobacteriaReferenceUnit(CIb-LIS),NationalInstituteofPublic
HealthandtheEnvironment,Bilthoven,TheNetherlands,7MRCNationalInstituteforMedicalResearch,London,UnitedKingdom
Mycobacteriumtuberculosisinfectsonethirdofthehumanworldpopulationandkillssomeoneevery15seconds.For
more than a century, scientists and clinicians have been distinguishing between the human- and animal-adapted
membersoftheM.tuberculosiscomplex(MTBC).However,allhuman-adaptedstrainsofMTBChavetraditionallybeen
consideredtobeessentiallyidentical.Wesurveyedsequencediversitywithinaglobalcollectionofstrainsbelongingto
MTBCusingsevenmegabasepairsofDNAsequencedata.WeshowthatthemembersofMTBCaffectinghumansare
more genetically diverse than generally assumed, and that this diversity can be linked to human demographic and
migratoryevents.Wefurtherdemonstratethattheseorganismsareunderextremelyreducedpurifyingselectionand
that,asaresultofincreasedgeneticdrift,muchofthisgeneticdiversityislikelytohavefunctionalconsequences.Our
findings suggest that the current increases in human population, urbanization, and global travel, combined with the
populationgeneticcharacteristicsofM.tuberculosisdescribedhere,couldcontributetotheemergenceandspreadof
drug-resistant tuberculosis.
Citation:HershbergR,LipatovM,SmallPM,ShefferH,NiemannS,etal.(2008)HighfunctionaldiversityinMycobacteriumtuberculosisdrivenbygeneticdriftandhuman
demography.PLoSBiol6(12):e311.doi:10.1371/journal.pbio.0060311
tuberculosiscomplex (MTBC)havetraditionally beenassumed
Introduction
to be essentially identical. This notion was primarily driven
Mycobacterium tuberculosis is a gram-positive bacterium and bytheresultsofearlystudiesthatrevealedverylowlevelsof
the causative agent of human tuberculosis. The worldwide DNA sequence variation in human MTBC [14,15]. More
emergence of multidrug-resistant strains of M. tuberculosis is recent surveys of global strain collections show that in fact
threatening to make tuberculosis incurable [1]. Although human MTBC consists of separate strain lineages associated
renewedeffortsarebeingdirected towardsthedevelopment with different regions of the world [16–20]. However, all of
of new tools to better control tuberculosis [2], much about thesestudieshaveimportantlimitationssuchthattheactual
the evolution of this obligate human pathogen remains phylogeneticdistancesand relativegeneticdiversitieswithin
unknown [3]. and between mycobacterial lineages have not been deter-
In 1898, Harvard pathologist Theobald Smith demonstra- mined[21,22].Specifically,thestudybyBrudeyetal.[17]used
ted that tubercle bacilli isolated from humans differed the standard molecular epidemiological method known as
significantlyfrombacilliisolatedfromcattleintheircapacity spoligotypingtodeterminetheglobalpopulationstructureof
tocausediseaseindifferentanimalspecies[4].Eventually,the M. tuberculosis. However, because this technique indexes
two bacilli were granted separate species status, with M. genetic diversity based on the presence or absence of a
tuberculosis designating the typical human pathogen, and repetitive sequence at a single locus (the ‘‘direct repeat
Mycobacterium bovis referring to the bovine form [5]. Because region’’ of M. tuberculosis), which is prone to convergent
M.bovishasthecapacitytocausediseaseinavarietyofanimal
species,includinghumans,itwasoriginallythoughttoexhibit
AcademicEditor:MartinJ.Blaser,NewYorkUniversitySchoolofMedicine,United
a much broader host range than M. tuberculosis. However,
StatesofAmerica
recent comparative genomic analyses have revealed such a
Received July 10, 2008; Accepted October 31, 2008; Published December 16,
high degree of genetic diversity in M. bovis that modern
2008
population geneticists now consider the species to be
Copyright:(cid:2)2008Hershbergetal.Thisisanopen-accessarticledistributedunder
comprised of several ecotypes, each of which is adapted to
thetermsoftheCreativeCommonsAttributionLicense,whichpermitsunrestricted
particularanimalhostspecies[6–10].Someoftheseecotypes use,distribution,andreproductioninanymedium,providedtheoriginalauthor
andsourcearecredited.
have been given distinct species designations. For example,
Mycobacteriummicrotiisapathogenofvoles[11],Mycobacterium Abbreviations:MTBC,Mycobacteriumtuberculosiscomplex;SNP,singlenucleotide
polymorphism
pinnipedii a pathogen of seals and sea lions [12], and
Mycobacteriumcaprae apathogen ofgoats[13]. *Towhomcorrespondenceshouldbeaddressed.E-mail:[email protected]
By contrast, the human-adapted members of the M. [Theseauthorscontributedequallytothiswork.
PLoSBiology|www.plosbiology.org 0001 December2008|Volume6|Issue12|e311
FunctionalDiversityinM.tuberculosis
Author Summary froma representative collection of108 MTBC strains (Table
S1). This collection included 99 human-adapted strains that
Tuberculosis remains a worldwide public health emergency. The were selected to represent the broadest geographic and
emergenceofdrug-resistantformsoftuberculosisinmanypartsof genetic diversity from a global collection of 875 strains
the world is threatening to make this important human disease characterized previously by the analysis of deletions across
incurable.Eventhoughmanyresourcesarebeinginvestedintothe thegenome[19].Anadditionalsevenstrainswereselectedto
developmentofnewtuberculosiscontroltools,westilldonotknow
represent four animal-adapted ecotypes, including M. bovis,
the extent of genetic diversity in tuberculosis bacteria, nor do we
M. microti, M. pinnipedii, and M. caprae. We also included the
understand the evolutionary forces that shape this diversity. To
vaccine strain M. bovis BCG Pasteur and one strain of
address these questions, we studied a large collection of human
tuberculosis strains using DNA sequencing. We found that strains Mycobacterium canettii as our predicted outgroup. M. canettii is
originating in different parts of the world are more genetically formallyconsideredpartofMTBC.However,inthisstudywe
diverse than previously recognized. Our results also suggest that use ‘‘MTBC’’ to refer to all other members of MTBC,
muchofthisdiversityhasfunctionalconsequencesandcouldaffect excluding M. canettii. M. canettii strains have been shown to
the efficacy of new tuberculosis diagnostics, drugs, and vaccines. be more distantly related to the remaining MTBC than any
Furthermore, we found that the global diversity in tuberculosis two other MTBC strains are to each other [28]. For each of
strainscanbelinkedtotheancienthumanmigrationsoutofAfrica,
the108strainsincludedinthestudy,wedeterminedtheDNA
aswellastomorerecentmovementsthatfollowedtheincreasesof
sequenceof89genes,whichtogethercorrespondedto65,829
humanpopulationsinEurope,India,andChinaduringthepastfew
base pairs per strain, or 1.5% of the ;4.4 Mbp genome of
hundred years. Taken together, our findings suggest that the
MTBC (Tables S2 and S3, Figure S1). These 89 genes
evolutionarycharacteristicsoftuberculosisbacteriacouldsynergize
with the effects of increasing globalization and human travel to comprised housekeeping genes and antigens analyzed in the
enhancetheglobalspreadofdrug-resistanttuberculosis. early sequencing studies mentioned above [14,15], as well as
genesofspecialinterest,includingputativenewdrugtargets
[29–31], genes believed to be involved in latency and
evolution, this technique is of limited use for phylogenetic reactivation [32,33], DNA repair genes [34], genes encoding
andpopulationgeneticanalyses[16,18,20,21,23].Thestudyby anovelbacterialsecretionsystem[35,36],andnovelantigens
Baker et al. [16] used a multilocus sequencing approach to [37,38]. Although the selection of these 89 genes was not
study M. tuberculosis diversity, but because just seven genes random,thegenesaredistributedfairlyuniformlyaroundthe
were analyzed, only a small number of phylogenetically MTBC chromosome, thus covering all the main parts of the
informative single nucleotide polymorphisms (SNPs) were MTBCgenome (FigureS1).
identified.InthestudiesbyGutackeretal.[20]andFilliolet We first used the concatenated DNA sequences of the 89
al. [18], the authors used a very similar approach: they genesforeachofthe108strainsandconductedamaximum
compared the full genome sequences of MTBC strains parsimony analysis. Our analysis produced a single phyloge-
available at the time and identified a series of synonymous netic tree with a homoplasy index of 0.0043 (Figure 1). This
SNPs,whichtheyusedtogenotypelargecollectionsofstrains. phylogenetic tree is completely congruent with the one we
However, such approaches are known to lead to so-called constructedusingtheneighbor-joiningmethod(FigureS2)as
phylogenetic discovery bias and distorted phylogenetic wellaswiththedeletion-basedanalysiswereportedpreviously
inference [22,24,25]. In our previous study [19], we used (Figure1,TableS1)[19].Thenegligibledegreeofhomoplasy
genomicdeletions(largesequencepolymorphisms)toanalyze observed in our sequence-based tree, and the fact that our
aglobalcollectionofstrains.Eventhoughwewereabletouse sequence-based and deletion-based trees are congruent
these deletions to classify strains unambiguously, genetic further supports the highly clonal population structure of
distances based on genomic deletions are difficult to MTBC [39,40]. The primary branches of our sequence-based
interpret [3,21]. Finally, because of the inherent limitations tree are also consistent with earlier studies that classified M.
of the molecular markers used in all of the studies reviewed tuberculosis into ‘‘ancient’’ and ‘‘modern’’ forms based on the
above, the evolutionary processes that shape strain diversity presenceorabsenceofagenomicdeletionknownasTbD1[8].
in MTBC have not been adequately investigated; generally, However, because our new sequence data allow us to better
actualDNAsequencedataarepreferredforphylogeneticand interpretgeneticdistancesbetweenlineages,wefindthatthe
populationgenetic analyses [26,27]. difference between ancient and modern MTBC is more
Here we report our in-depth analyses of a large set of pronounced than one would assume on the basis of the
coding sequence data from a global collection of MTBC presenceofasinglegenomicdeletion(Figure1).
strains. These analyses reveal that the human-adapted
Human-Adapted MTBC Is More Genetically Diverse than
members of MTBC are more genetically diverse than
Assumed
generally recognized. We also demonstrate that genetic drift
is likely to be an important evolutionary force generating We further compared the genetic distances between
different strain lineages using our new sequence data.
diversity in MTBC, and that this diversity can be linked to
Overall, our analysis reveals greater genetic diversity of the
changes in human demography and to both ancient and
human-adapted organisms than previously appreciated.
recent humanmigrations.
Specifically, our sequence-based phylogeny shows that all
the animal-adapted members of MTBC form an ingroup
Results and Discussion
relative to the rest of the phylogeny. This suggests that even
The Global Phylogeny of M. tuberculosis though these animal strains belong to four distinct ecotypes
We investigated the genetic diversity within MTBC using adaptedtodistinctanimalhostspecies,theyrepresentonlya
seven megabases of DNA sequence data that we generated proportionofthegeneticdiversityfoundinallofthehuman-
PLoSBiology|www.plosbiology.org 0002 December2008|Volume6|Issue12|e311
FunctionalDiversityinM.tuberculosis
Figure1.MaximumParsimonyPhylogenyofM.tuberculosisComplexUsing89ConcatenatedGeneSequencesin108Strains
Thebranchesarecoloredaccordingtothemainlineagesdefinedpreviouslybasedonourgenomicdeletionanalysis,exceptfortheanimal-adapted
strains,whichareindicatedinorange.Themaincladesarelabeledaccordingtotheirdominanceinparticulargeographicareas.Thebranchleadingto
PLoSBiology|www.plosbiology.org 0003 December2008|Volume6|Issue12|e311
FunctionalDiversityinM.tuberculosis
M.canettiihasbeentruncatedinthefigurebecauseofthelargenumbersofchangesthatseparatethishypothesizedoutgroupfromtherestofthe
phylogeny(TableS4).Ancientandmodernstrainlineagesareindicated.ThegreenandbrownlineagescorrespondtostrainstraditionallyknownasM.
africanum[21].
doi:10.1371/journal.pbio.0060311.g001
adapted MTBC (Figure 1). In fact, the average pairwise This is substantially higher than that in most other bacteria
geneticdistanceamonganimal-adaptedstrainsisequaltothe [41] and is significantly higher than that observed for the
averagedistanceamonghuman-adaptedstrains(23.61versus SNPsspecific to M. canettii(G-test, p, 0.0001).
23.65 differences, respectively; Figure 2). Additionally, the The high dN/dS in MTBC could be a consequence of
distribution of all possible pairwise genetic distances within diversifyingselectiondrivenbyhostimmunity.Weaddressed
the human-adapted group was not significantly different thispossibilityby comparingtheratioofnonsynonymous to
from that within the animal-adapted group (Mann–Whitney synonymous SNPs in surface-exposed or excreted, virulence,
ranksumtest,p¼0.74). Takentogether, thisevidence shows and housekeeping genes (Table S2). Contrary to what we
that the genetic diversity among the human-adapted mem- would expect if immune selection was responsible for the
bers of MTBC is just as pronounced as that among the highdN/dSinMTBC,wefoundthatallthreegeneclasseshad
animal-adapted strains tested. Although the human MTBC comparable ratios of nonsynonymous to synonymous SNPs
strains used in this study are not a true population survey (Table S6), similar to the situation reported for Salmonella
because they were selected to maximize either diversity or typhi [42]. It is therefore unlikely that diversifying selection
geographicaldistribution,theseobservationssuggestthatthe explains the high ratio of dN/dS unless housekeeping
genetic diversity among the human-adapted members of proteins,virulenceproteins,andsurface-exposedorsecreted
MTBC is just as pronounced as that among the animal- proteins are all equally and strongly exposed to immune
adapted strains. This is especially surprising given the fact surveillance and areallinvolved inthe creation ofantigenic
that the latter belong to four distinct ecotypes adapted to variation.Thislatterpossibilityseemsunlikelybecause,even
distinctanimal hostspecies. though M. tuberculosis is an intracellular pathogen, during
various stages of the natural history of human tuberculosis
Purifying Selection Is Severely Reduced in MTBC
(for example transmission and infection processes, and
Next,weusedoursequencedatasettoestimatetheroleof cavitary and disseminated disease) bacilli are in fact located
purifying selection in the evolution and genetic diversity of extracellularly. Furthermore, there is increasing evidence
MTBC.Onecommonlyusedmethodofexaminingthedegree that M. tuberculosis interacts with both the humoral and the
ofpurifyingselection acting onsequences is tocalculate the cellular immune system as evidenced by antibody-based
ratio of the rates of nonsynonymous and synonymous serological recognition of M. tuberculosis antigens in patient
changes (dN/dS). In the absence of selection this ratio is sera[38].Thissuggeststhatsurface-exposedorexcretedgene
expected to near unity. Purifying selection is expected to products in MTBC should, in principle, be under stronger
reduce this ratio, while positive selection is expected to immunesurveillance than structuralgenes.
increaseit. Inotherorganisms,highdN/dSratiosareoftenconsidered
Wediscoveredatotalof488SNPsinour108strains.TheM. to indicate a reduction in purifying selection [26,27].
canettiistraindifferedfromeachoftheotherMTBCstrainsat However, such high dN/dS values may also stem from the
129–145sites(0.2%oftheexaminedsites,TableS4),whilethe close relatedness of the MTBC strains. Rocha et al. pointed
maximum number of SNPs between any two other MTBC outthatdN/dSisoftenhigherincasesinwhichtheorganisms
strainswas46(0.07%oftheexaminedsites).Thissupportsthe beingcomparedareverycloselyrelated[43].Forsuchclosely
useofM.canettiiasacloselyrelatedoutgroup.Thecomparison related organisms the dN/dS ratio may be elevated because
oftheproteinsintheM.canettiistrainwiththemajority-rule slightly deleterious nonsynonymous mutations that are
consensus of the proteins in the remaining MTBC strains destined to be lost during long periods of time may not yet
revealed 125 nucleotide differences, of which 43 were have been removed by purifying selection [43]. If dN/dS
nonsynonymous. This corresponds to the ratio of the rates within MTBC is inflated because of nonsynonymous muta-
ofnonsynonymousandsynonymouschanges(dN/dS)of0.18. tions being slightly deleterious, we expect the frequency
This ratio is similar to that previously observed between M. distribution for nonsynonymous SNPs to be skewed towards
tuberculosisandtherelativelymoredistantlyrelatedMycobacte- low frequencies when compared with synonymous SNPs.
riumavium(0.17)[41].Itisalsosimilartoratiosweobtainedby Using our large dataset of sequences it was possible to test
pairwise comparisons of the two fully sequenced M. avium thisdirectly.Wefoundthatalthoughthefrequencydistribu-
strains: M. avium 104 and M. avium paratuberculosis (dN/dS¼ tionofnonsynonymousSNPswashighlyskewedtowardslow
0.17),andforthetwofullysequencedstrains:Mycobacteriumsp. frequencies, this was equally true for the synonymous SNPs
JLSandMycobacteriumsp.MCS(dN/dS¼0.15).Thisappearsto (Figure3),withnodifferenceintheproportionofsingletons
representthegeneraldN/dSwithintheActinobacteria,aswe inthetwotypesofSNPs(G-test,p¼0.66).Similarly,whenwe
conclude from an examination of pairwise genome-wide examined the frequency of 68 genomic deletions in 100
comparisons of different fully sequenced Actinobacteria clinical MTBC strains reported earlier [44], we found that
genomes(unpublisheddata).Itthusappearsthatthestrength these deletions were also not more likely to occur as
ofselectionactingonM.canettiimaybesimilartothatacting singletonswhencomparedwiththesynonymousSNPsfound
onotherMycobacteriaandActinobacteria. here (G-test, p ¼ 0.32). These results indicate that neither
Outofthe370SNPsthatsegregatedamongtheremaining nonsynonymous SNPs nor deletions show any evidence of
107 MTBC strains (i.e., excluding M. canettii; Table S5), 231 being slightly deleterious within MTBC (i.e., they behave as
(62%)werenonsynonymousand139(38%)synonymous.The selectively neutral). Second, if slightly deleterious mutations
average pairwise dN/dS ratio for the MTBC strains was 0.57. contributetotheincreaseddN/dSinMTBC,weexpectdN/dS
PLoSBiology|www.plosbiology.org 0004 December2008|Volume6|Issue12|e311
FunctionalDiversityinM.tuberculosis
Figure2.FrequencyofPairwiseGeneticDistancesAmong108StrainsoftheMTBC
Pairwisecomparisonsofgeneticdistancesbetweenhuman-adaptedstrainsareindicatedinyellow,thosebetweenanimal-adaptedstrainsinblue,and
thosebetweenhuman-andanimal-adaptedstrainsinred.Forbetterillustration,thefrequenciesofgeneticdistancesintheanimal-to-animalstrain
comparisonsandtheanimal-to-humanstraincomparisonsweremultipliedby20and3,respectively.However,forthestatisticalanalysis(seemaintext),
theactualfrequencieswereused.
doi:10.1371/journal.pbio.0060311.g002
to be lower among common SNPs, as slightly deleterious and36%werevariable(Table1).Weexpectthatmutationsat
polymorphisms are less likely to reach high frequencies. the conserved positions should have stronger functional
However, we found that the proportion of nonsynonymous effects and be more deleterious on average than mutations
SNPsamongthe75SNPsthatwerepresentinmorethanfive at the variable positions. Consistent with this view, non-
strains was no different from the proportion in the entire synonymouschangesspecifictoM.canettiipredominantlyfall
dataset(63%comparedwith62%,G-test,p¼0.97).Tofurther intovariablepositions(72%,G-test,p,0.0001,Table1).We
test our hypothesis, we compared the ratio of nonsynon- expect that weakened purifying selection in MTBC should
ymous to synonymous SNPs in different parts of our allowmoreaminoacidchangesinMTBCtobeobservedatthe
phylogenetic tree (Figure 1). We found no significant differ- conserved (i.e., constrained) positions. As expected, in con-
enceinthisratiobetweeninternalandexternalphylogenetic trasttothechangesinM.canettii,themajority(58%)ofamino
branchesorbetweentheancientandmodernMTBClineages acid mutations in MTBC fall into the conserved positions
(Table S6). These results together indicate that, contrary to (Table1).Infact,theproportionsofaminoacidmutationsin
the suggestion of Rocha et al. [43], the high dN/dS values MTBC that fall into variable and conserved positions is not
within MTBC are not an artifact of the close relatedness of significantly different from that expected if purifying selec-
the MTBC strains, but are in fact likely to be the result of a tion in MTBC was no longer making a distinction among
reductionin selectiveconstraint. mutations in these two classes of sites (G-test, p ¼ 0.1). A
To further probe the observed reduction in selective differentwaytolookatthisistocalculateameasuresimilarto
constraint,weclassifieddifferentpositionswithinthestudied dN/dS that compares the rates of variable and conserved
proteinsaccordingtotheirpatternsofconservationinallten amino acid changes. P can be defined as the number of
var
fully sequenced mycobacterial species that are distantly changesthatfallatvariablesitesdividedbythetotalnumber
related to the MTBC/M. canettii strains. We first searched for of variable sites. Similarly P is defined as the number of
cons
orthologsofthe89sequencedgenesinthesetenmycobacte- changes at conserved sites divided by the total number of
rialspecies.For62ofthegeneswecouldfindorthologsinat conservedsites.P /P equals0.22fortheM.canettiichanges
cons var
leastfiveofthetenspecies(TableS2).Wealignedtheprotein and 0.77 for the MTBC changes. This analysis further
sequences of these 62 genes in the mycobacterial species in underscoresthelowlevelofpurifyingselectioninMTBC.
which they were present (excluding MTBC/M. canettii ortho-
logs)andexaminedthelevelofconservationateachalignment TheObservedReductioninSelectionIsNotSpecifictothe
position.Theaminoacidpositionswerethendividedintotwo 89-Gene Dataset Used
groups: (i) ‘‘conserved’’ positions—positions that either have Toevaluatewhetherthe89genesweselectedwereindeed
identicalaminoacidsinalltheexamineddistantmycobacte- representative of the whole M. tuberculosis genome, we
rial species or vary only among biochemically similar amino repeated some of our analyses using completed MTBC
acids, and (ii) ‘‘variable’’ positions—the least constrained genomes. Currently, full genomic sequences are available
positions that show some radical amino acid changes (see for six MTBC strains: M. tuberculosis H37Rv, M. tuberculosis
Methods).64%ofpositionsintheseproteinswereconserved H37Ra, M. tuberculosis F11, M. tuberculosis CDC1551, M. bovis,
PLoSBiology|www.plosbiology.org 0005 December2008|Volume6|Issue12|e311
FunctionalDiversityinM.tuberculosis
Figure3.FrequencyDistributionofNonsynonymousMutations,SynonymousMutations,andGenomicDeletionsinMTBCStrains
ForeachsitethatunderwentamutationwithinMTBCweobservedtwostates:thederived(ormutated)state,andtheancestralstate.Forsynonymous
andnonsynonymousmutationswededucedtheancestralstatebasedonthesequenceoftheoutgroupstrainM.canettii.Forthegenomicdeletions
(originallydescribedin[44]),theancestralstateisthenon-deletedstateandthedeletedstateisderived.Thefrequenciesdepictedherearethoseofthe
mutatedstates.Allthreetypesofmutationsshowsimilarfrequencydistributions,whichsuggeststhatmostnonsynonymousmutationsanddeletions
arenotslightlydeleteriousbutratherselectivelyneutralwithinMTBC.Blue,nonsynonymousmutations;red,synonymousmutations;yellow,genomic
deletions.
doi:10.1371/journal.pbio.0060311.g003
and M. bovis BCG Pasteur 1173P2. We conducted genome- of recombination or horizontal gene transfer in MTBC, for
wide pairwise comparisons of the DNA sequences of all whichsomeevidencehasbeenpresented[47,48].Thesegenes
orthologousproteinpairsbetweentheH37Rvstrainandthe do not seem to be representative of the genome as a whole,
remainingfivestrains.Thisallowedustoobtainthenumber given the small number of mutations in all other genes and
of synonymous and nonsynonymous differences for all thelarge number ofmutations within thesegenes; including
orthologs across the genome. As strains belonging to MTBC them may skew the analysis. We therefore removed these
show very low levels of diversity, it is impossible to use genes prior to summing the mutations. For all five compar-
pairwise comparisons to calculate dN and dS for specific isonstheratioofnonsynonymoustosynonymousdifferences
genes.However,itispossibletosumthemutationsacrossthe wasverysimilartotheratiofoundwithinour89genedataset
genome and examine whether the ratios of nonsynonymous (Table 2, G-test, P . 0.92 for all five comparisons). It thus
andsynonymousmutationsthatwefoundforthe89genesin appearsthatthereducedselectionweobservedisnotlimited
our datasetare representativeofthe entire genome. Forthe tothe genes withinour dataset.
vast majority of orthologs pairs we found zero to two This genome-wide analysis also provides a further indica-
differences between the H37Rv strain and any of the other tionthattheobservedhighratioofdN/dSisnottheresultof
five strains (.97% for all comparisons). However, in each immunesurveillance.Similarlytothemajorityofthegenesin
comparisontherewasasmallnumberofgenes(lessthanten) thegenome,andinsharpcontrasttopks12andsomegenesof
thatwereclearlyoutliers.Thesegenesshowedamuchhigher thePEandPPEfamiliesthatareknowntobeunderimmune
number of differences than the genome in general. Many of surveillance,the89geneswesequencedshowonlyzerototwo
these genes, such as pks12 and members of the PE and PPE differences between the H37Rv strain and any of the other
gene families, were previously implicated as involved in five fully genome sequenced strains. This, in addition to our
antigenicvariation[45,46]andmaybesubjecttodiversifying finding that the high dN/dS ratio is not specific to the 89
selection. Alternatively, they could represent rare examples analyzedgenes,butratherischaracteristicofthegenomeasa
Table1. Classification ofMutationsAccording to Degreeof ConservationinDistantly RelatedMycobacterial Species
Position Number(%)ofAll Number(%)of Number(%)ofDifferences
Classification PositionsExamined MutationsinMTBC betweenMTBCandM.canettii
Conserved 10,385(64%) 97(58%) 9(28%)
Variable 5,859(36%) 71(42%) 23(72%)
doi:10.1371/journal.pbio.0060311.t001
PLoSBiology|www.plosbiology.org 0006 December2008|Volume6|Issue12|e311
FunctionalDiversityinM.tuberculosis
Much of MTBC Diversity Has Functional Consequences
Table2.SynonymousandNonsynonymousMTBCDifferencesin There are two possibilities regarding the nature of the
the89GeneDatasetandinGenome-WidePairwiseComparisons observed reduced constraint. One possibility is that the
reduction in purifying selection is pathway specific and that
Comparison Number(%)of Number(%)of itaffectsonlyspecificmutationsthatnolongerhavefunctional
Synonymous Nonsynonymous consequencesasaresultofchangesintheecologicalrequire-
Differences Differences
mentsoftheMTBCstrains.Ifthisweretrue,onewouldexpect
to see a relaxation of purifying selection in only a subset of
89genedataset 139(37.6%) 231(62.4%)
genes, and not in housekeeping genes that are generally
H37Rv/F11 223(37.2%) 376(62.8%)
H37Rv/CDC1551 240(37.1%) 407(62.9%) universally required. However, the reduction of purifying
H37Rv/H37Ra 29(38.2%) 47(61.8%) selection within the MTBC appears to be genome-wide and
H37Rv/M.bovis 607(37.3%) 1,020(62.7%) occurs across different classes of genes, including house-
H37Rv/M.bovisBCG 571(37.2%) 962(62.8%)
keeping genes. This suggests the alternative possibility, in
which a genome-wide reduction in the efficacy of purifying
doi:10.1371/journal.pbio.0060311.t002 selection has occurred. The weakened selective constraint in
MTBC appears to allow amino acid changes that would be
whole, indicates that immune surveillance is not likely to
removed by purifying selection in M. canettii and other
explain the observed high dN/dS values, unless it affects the
mycobacteriatopersistwithinMTBC.Suchchangesarelikely
entiretuberculosisgenome.Furthermore,thefactthatsome
to have deleterious functional consequences in M. canettii, as
genesthatareknowntobeunderimmunesurveillanceshow
otherwisetheywouldnothavebeenremovedbyselection.Most
much higher levels of diversity than most other genes in the
of these mutations are still likely to have similar functional
genomefurtherindicatesthatmostgenesarenotaffectedto
effects in the other MTBC strains. However, because of the
suchahigh extent byimmune surveillance.
reducedefficiencyofselectiontheyaremorelikelytopersist.
WefurtherusedthesixfullysequencedMTBCgenomesto
It is possible to roughly estimate the proportion of the
calculate P /P at a genome-wide level. We created
cons var differences within MTBC that are likely to carry such
multiplesequencealignmentsofalloftheannotatedproteins
functionalconsequences.Onthebasisofour89genedataset,
oftheH37Rvstrainthathaveclearorthologsintheotherfive
we found that in M. canettii, 23 mutations fell into variable
MTBC genomes. We also searched for orthologs of these
sites and nine, or ;2.6 times fewer, fell into conserved sites.
sequences in the ten more distantly related fully sequenced
Had selection been as strong in MTBC as it is in M. canettii,
mycobacterialstrainsmentionedearlier.Forthisanalysis,we
and assuming conservatively that the mutations at the
selectedgenesforwhichwecouldfindclearorthologsinallof
variable sites are entirely neutral,we would expect 2.6 times
the six MTBC strains and create multiple sequence align-
fewer mutations to be observed at constrained sites within
ments in these six strains that contained no gaps. We also
MTBC as well. Considering that there are 71 mutations at
requiredthatwefindorthologsforthesegenesinatleastfive
variable sites in MTBC we would expect ;27 mutations at
of the more distantly related mycobacterial strains. 1,970
conserved sites within MTBC. Instead, we find 97 such
genesfilledtheserequirements.WecouldfindnoSNPswithin
mutations. This indicates that close to 72% of mutations at
the six MTBC strains for 1,289 of these genes and thus
conserved sites (42% of the amino acid mutations overall)
removed them from the analysis. We further removed four
within MTBC would have been removed by selection in M.
genesbecausetheywereclearlyoutliersandhadwelloverten
canettii.Becausesuchaminoacidmutationsarelikelytoaffect
SNPs each, whereas the vast majority of genes had three or
fewer SNPs. These four genes include three-members of the proteinfunction,thissuggeststhatabout40%ofaminoacid
PE and PPE protein families, which have been implicated as changesinMTBC havefunctional consequences.
beinginvolvedinantigenicvariation[45],aswellasthepks12 To put this in context, we can consider the number of
gene,whichhasbeenshowntoproduceapolyketidethatisan nonsynonymous differences found between the fully se-
important T cell antigen [46]. After removing these four quencedM.tuberculosisH37Rvandothertwofullysequenced
genes,wewereleftwith677genesinwhichtherewasatleast strainsofM.tuberculosispresentinourdataset:CDC1551and
oneMTBCSNP.Withinthesegenes145,020ofthesiteswere F11.Weconsideredonlygenesthatareclearorthologs,could
conserved in the distantly related Mycobacteria (49%) and be aligned across their entire protein sequence, and had no
149,728 were variable (51%). 448 of the MTBC SNPs found gapsintheirDNAalignmentthatcouldrepresentframeshift
within these genes were at conserved sites (47%) and 515 of events.Again,weremovedgenesthatwereclearoutlierswith
theSNPsfellatvariablesites(53%).Aswasthecaseforthe89 respecttothenumberofmutationstheycontain(morethan
genes in our dataset, there is no significant difference tenmutations).Withintheremainingproteinpairs,wefound
between the proportion of SNPs falling in variable and 423 and 410 amino acid differences between H37Rv and
conservedsitesandtherandomproportionwewouldexpect CDC1551, and H37Rv and F11, respectively. Note that the
ifselectionwerenotdiscriminatingbetweenthetwotypesof genome-sequencedstrainsdonotrepresentthefulldiversity
sites (G-test, p¼0.1). Furthermore, when we calculated the within MTBC, because they all belong to the same strain
genome-wideP /P wefounditwasevenhigher(0.9)than lineage(redinFigure1).However,ifthesequencedproteins
cons var
for the 89 gene analysis (0.77). Taken together, our genome- are representative of the genome as a whole, these numbers
wide analyses show that the severe reduction in purifying andthediversityinourdatasetcanbeusedtoconservatively
selection we observed is not specific to the 89 genes in our estimate the range of functional amino acid differences
main dataset. Rather there appears to be a genome-wide betweenanaveragepairofMTBCstrains.Forthe89geneswe
reductionin purifyingselection in MTBC. sequenced, the number of differences between strain pairs
PLoSBiology|www.plosbiology.org 0007 December2008|Volume6|Issue12|e311
FunctionalDiversityinM.tuberculosis
Figure4.Out-Of-And-Back-To-AfricaScenariofortheEvolutionaryHistoryofHuman-AdaptedMTBC
(A)Globalhumanpopulationsizeduringthelast50,000y.Thelettersabovethegraphindicatethetimeperiodscorrespondingto(B),(C),and(D),
respectively (data source: http://en.wikipedia.org/wiki/Image:Population_curve.svg). (B) Hypothesized migration out of Africa of ancient lineages of
MTBC.Coloredarrowscorrespondtothesixmainhuman-adaptedMTBClineagesshowninFigure1.Thehypothesizedcommonancestorofthethree
modernlineages(inred,purple,andblue)isindicatedinblack.
(C)Recentincreaseofglobalhumanpopulation.Eachdarkgreydotcorrespondsto1millionpeople(datasource:http://www.pbs.org/wgbh/nova/
worldbalance/numb-nf.html).ThepopulationincreasewasstrongestinWesternEurope,India,andEastAsia(FigureS3).Thesethreegeographicregions
areeachassociatedwithoneofthethreemodernMTBClineages(red,purple,andblue).Recenthumanmigration,trade,andconquesthavepromoted
globalspreadofthesemodernMTBClineages.
(D)Thehumanpopulationhasreached6billion.Thedistributionofthesixmainhuman-adaptedMTBClineagesweobservetodayisshown(colors
correspondtoFigures1andS2;basedondatafrom[19,21]).
doi:10.1371/journal.pbio.0060311.g004
rangedfromzeroto46(Figure2),withanaverageof25.For of infected individuals. It has been shown that different
these genes we find 17 differences between H37Rv and lesions in lungs of tuberculosis patients can harbor genet-
CDC1551,and13differencesbetweenH37RvandF11.From ically distinct subclones of a particular infecting strain [49].
thesenumberswecanestimateconservativelythatanaverage Finally, clonal organisms are also prone to selective sweeps
pair of MTBC strains should have around 300 functional [6]. All of these features lead to small effective population
differences between them. At the same time, the least sizes, in which the effects of random genetic drift are
diverged strains may have close to no functional differences increased compared with those of natural selection
whereas the most diverged strains may differ at around 500 [6,26,27]. Random genetic drift allows mutations to reach
functionalsites(seeMethods).Furthermore,becauseweused high frequencies that in organisms with large effective
a highly conservative estimate of the number of amino acid population sizes would be eliminated by natural selection.
differences per genome and took only point mutations into Such mutations are likely to affect protein function delete-
account, this is likely to be a highly conservative estimate of riously (i.e., they are functional), which is why they are
thenumber offunctional differences. eliminatedinotherorganisms.Bycontrast,thelargeamount
Taken together, our results strongly suggest that the high of predicted functional variation we observed at conserved
valueofdN/dS inMTBC resultsfrom asignificantreduction sitesinMTBCshowsthatnatural selectionisnotasefficient
inpurifyingselection.Thisreducedselectiveconstraintmost atremovingdeleteriousmutationsinthisorganism.Thishigh
probably results from a numberof factors, all of which tend functional diversity in MTBC reported here further stresses
to reduce an organism’s effective population size. These the need to consider strain diversity in the development of
factors are high clonality (i.e., virtual absence of horizontal newtuberculosis diagnostics, drugs, and vaccines[21].
gene exchange) and the serial transmission bottlenecks
characteristic of this pathogen. These serial bottlenecks are An ‘‘Out-Of-And-Back-To-Africa’’ Scenario for Human
particularly tight in M. tuberculosis, because a single cell is MTBC
enough to establish a new infection. Furthermore, the Inadditiontonaturalselectionandgeneticdrift,migration
population structure of human MTBC is highly subdivided, and demographic changes can affect the generation and
both geographically (Figures 1 and 4D) and within the lungs distribution of genetic diversity in a particular organism.
PLoSBiology|www.plosbiology.org 0008 December2008|Volume6|Issue12|e311
FunctionalDiversityinM.tuberculosis
Thus, we decided to use our new sequence-based diversity humantraderoutesgobackseveralmillenniaandthuscould
dataset to further explore the complex interactions between havecontributedtothespreadofMTBC,itappearsthatthe
mycobacterial evolution and human population growth and large-scale migration events discussed above, which were
travel[50,51].Overall,ourdatasupportahypothesized‘‘Out- partiallydrivenbythemorerecentanddramaticincreasesin
of-and-back-to-Africa’’ scenario of the phylogeography of populationdensitiesinEurope,India,andEastAsia(Figures
human tuberculosis (Figures 1 and 4D) [19,21]. Specifically, 4A and S3), were key for the spread of the modern MTBC
mostevidenceindicatesthatMTBCoriginatedinAfrica.For lineages. Importantly, in contrast to the ancient migrations
example, Africa is the only continent in which all six major that followed mostly land routes, modern waves of spread
humanMTBClineagesoccur,andthedeeplyrootedlineages increasinglyfollowedroutesofocean-goingships(Figure4C)
are present exclusively in West Africa (brown and green in toarrive totheircontemporary distribution (Figure 4D).
Figures1and4).AnAfricanoriginforhumantuberculosisis Thequantitativenatureofournewdiversitydatapermitsa
alsoconsistentwithapreviousreportonM.canettiiandother more rigorous analysis of this out-of-and-back-to-Africa
so-called ‘‘smooth tubercle bacilli’’ that share a remote hypothesis.Theorypredictsthat foranorganismpopulating
common ancestor with the other MTBC and which are the world from a given point of origin, we will find a
primarily associated with countries at the Horn of Africa correlation between genetic variability and geographic
[28,52]. During hunter–gatherer times, human populations distance traveled [59,60]. Thus, if our hypothesis is correct,
remained small and geographically scattered (Figure 4A), we expect the genetic difference between any two strains to
which favored within-family transmission and the resulting increase with geographic distance between where they
stable host–pathogen associations [19,21,39]; it has been originated. For all isolates, we used the haversine method
postulated that the characteristic latency period in human with and without waypoints [59] to define pairwise distances
tuberculosis might be an adaptation of M. tuberculosis to low over land routes and over water routes, respectively.
hostdensities[53].Thus,wewouldexpecttheancienthuman Incorporation of waypoints into the calculation of geo-
migrationsoutofAfrica;50,000yago[54]tobereflectedin graphicdistancesallowsfortheadditionaldistancestraveled
thepopulationstructureofhuman-adaptedMTBC,similarto
via land routes during ancient human migrations. More
thestructurethathasbeenobservedforHelicobacterpylori[55].
direct geographic distances calculated without waypoints
In our dataset, this is seen in the distribution of the three
correspond to the water routes along which more recent
deeply rooted ancient strains. For example, the contempo-
human travel occurred once ocean-going ships became
rary distribution of the ancient pink lineage around the
available.Foreachtypeofgeographicdistance,wecomputed
Indian Ocean (Figures 1 and 4D) corresponds to the earliest
the coefficientof linearcorrelation between geographic and
spread of modern humans out of Africa [54]. Importantly,
genetic distances among ancient and modern strains, and
earlymigrationoutofAfricawouldhaveoccurredoverland
tested its statistical significance. For all human-adapted
asocean-goingtransporttechnologiesweredevelopedonlyat
MTBC strains we found a significant correlation between
amuch later timepoint inhumanhistory [51,54].
pairwise geographic and genetic distances, indicating that
Ancient MTBC Spread by Land, Modern MTBC by Sea eachgroupofstrainscarriesasignaloftheoverallpatternsof
human demography and migration (Figure 5A and 5B).
We hypothesize that further overland migration seeded
Furthermore, given the difference that we hypothesize
Western Europe, Northern India, and East Asia with what
betweenthemigratoryroutesofancientandmodernstrains
wouldbecomethethreemodernM.tuberculosislineages(red,
we would expect that the correlation would be stronger for
purple,andbluelineagesinFigures1and4B).Concurrently,
human populations started to increase dramatically and ancient strains spreading via land routes and for modern
disproportionately in these regions providing an ecological strainsspreadingviasearoutes;thisnotionwassupportedby
niche for the clonal expansion of these three strain lineages ourstatistical analysis(Figure 5Aand5B).
(Figures 4A, 4C, and S3). During the past few centuries with In summary, our quantitative analysis enabled by the
large-scalehumanmigration,tradeandconquestoutofthese sequence-based diversity data of MTBC strains supports the
threeareas[51],thesecladesfurtherspreadtootherareasof hypothesized relationship between the evolution of a very
the world (Figure 4C). Specifically, the presence of Euro- successfulhumanpathogenandpasthumanmigrations.Itis
American (red) strains on the American continent can be likelythatmorerecentairtravelmightalsobecontributingto
explained in terms of the exodus from the over-populated theglobalspreadofMTBCvariants,andpartofthesignalwe
EuropeancitiestoAmericaattheendofthe19thcentury—a detectedforthespreadofMTBCbywaterroutesmightinfact
‘‘vastmovementthatdwarfedallearliermigrations,’’accord- reflectmorerecentspreadbyairroutes.However,asdiscussed
ing to the historian William McNeill [51]. Furthermore, the above, it appears that changes in the global population
presence of strains from this cluster in Africa, Asia, and the structure of MTBC are linked to significant demographic
Middle East is in agreement with the cases of European changesandlarge-scalemovementsinhumanpopulations.On
colonizationofthepost-Columbianera.Inasimilarvein,the the long-term, increasing globalization and air travel could
presence of strains from the Indian-East-African (purple) leadtohomogenizationoftheglobalpopulationstructureof
lineageinEastAfricacanbeexplainedintermsoftherecent MTBC. Alternatively, it is conceivable that because of the
history of migration to this region from the Indian subcon- distinct evolutionary trajectories of ancient and modern
tinent [56]. Finally, the presence of the East Asian/Beijing MTBC lineages, fundamental differences in the optimal
(blue) lineage in South Africa is best explained by the virulencemightlimitthesuccessofancientlineagescompared
relatively recent import of slaves from Southeast Asia by withmodernstrains[61].Insupportofthispossibility,arecent
Dutch colonialists and later immigration of Chinese work- studyinTheGambiareportedthatmodernstrainswerethree
forces to South African gold mines [57,58]. Although some timesmorelikelytocauserapidprogressiontoactivedisease
PLoSBiology|www.plosbiology.org 0009 December2008|Volume6|Issue12|e311
FunctionalDiversityinM.tuberculosis
Figure5.CorrelationsbetweenGeneticandGeographicDistancesinAncientandModernMTBCStrains
To determine the route by which human-adapted MTBC strains were globally dispersed, we sought correlations between genetic and geographic
distances in ancient (A) and modern (B) MTBC strains. The shortest distances between geographic locations via water routes (pink diamonds) or
continentalwaypoints(i.e.,landroutes;bluesquares)areplottedagainstthecorrespondinggeneticdistancesexpressedasthenumberofSNPs.In
ancientstrains,thecorrelationusingeitherroutewashighlystatisticallysignificant(p,0.0001forboth),butthecorrelationusinglandrouteswas
slightlystronger.Bycontrastinmodernstrainsthecorrelationusingwaterrouteswasstrongerandmorestatisticallysignificantcomparedwithland
routes(p,0.0001andp¼0.05,respectively).Significantpvaluesandr2between0and1forallfouranalysessupporttheassociationbetweenMTBC
dispersalandhumanmigration,whichisconsistentwithearlyspreadofMTBCoutofAfrica.Therelativemagnitudeofr2forlandroutesversuswater
routessuggeststhatspreadofancientstrainsoccurredoverland(i.e.,throughancienthumanmigration)andmodernstrainsoverwater(i.e.,through
recenthistoricalmigration,trade,andconquest).
doi:10.1371/journal.pbio.0060311.g005
insecondarytuberculosispatientsthanstrainsthatbelonged Drug-resistant strains of M. tuberculosis can exhibit fitness
totheancientgreenlineage[62]. defects as a function of the specific drug-resistance-confer-
ring mutation and strain genetic background [63,64]. How-
Concluding Remarks
ever, reduced selection may allow strains harboring ‘‘costly’’
In conclusion, our analysis shows that human-adapted
drug-resistance-conferring mutations to persist, even in the
MTBCismoregeneticallydiversethantraditionallyassumed
and is under greatly reduced selective constraint. This absence of antibiotics. Furthermore, persistence in the
reduced selection pressure may have implications for the population of mutants resistant to a single drug, even if
emerging epidemic of multidrug-resistant tuberculosis [1]. relatively unfit, may contribute to the emergence of multi-
PLoSBiology|www.plosbiology.org 0010 December2008|Volume6|Issue12|e311