Efficient Inference of Recent and Ancestral Recombination within Bacterial Populations The

17 Pages·2017·1.3 MB·English

Efficient Inference of Recent and Ancestral
Recombination within Bacterial Populations
Citation
Mostowy, Rafal, Nicholas J. Croucher, Cheryl P. Andam, Jukka Corander, William P. Hanage,
and Pekka Marttinen. 2017. “Efficient Inference of Recent and Ancestral Recombination within
Bacterial Populations.” Molecular Biology and Evolution 34 (5): 1167-1182. doi:10.1093/molbev/
msx066. http://dx.doi.org/10.1093/molbev/msx066.
Published Version
doi:10.1093/molbev/msx066
Permanent link
http://nrs.harvard.edu/urn-3:HUL.InstRepos:32630649
Terms of Use
This article was downloaded from Harvard University’s DASH repository, and is made available
under the terms and conditions applicable to Other Posted Material, as set forth at http://
nrs.harvard.edu/urn-3:HUL.InstRepos:dash.current.terms-of-use#LAA
Share Your Story
The Harvard community has made this article openly available.
Please share how this access benefits you. Submit a story .
Accessibility
Efficient Inference of Recent and Ancestral Recombination
within Bacterial Populations
Rafal Mostowy,1 Nicholas J. Croucher,1 Cheryl P. Andam,2 Jukka Corander,3,4 William P. Hanage,2 and
Pekka Marttinen*,5
1DepartmentofInfectiousDiseaseEpidemiology,St.Mary’sCampus,ImperialCollegeLondon,London,UnitedKingdom
2DepartmentofEpidemiology,HarvardTHChanSchoolofPublicHealth,CenterforCommunicableDiseaseDynamics,Boston,MA
3DepartmentofMathematicsandStatistics,HelsinkiInstituteforInformationTechnologyHIIT,UniversityofHelsinki,Helsinki,Finland
4DepartmentofBiostatistics,UniversityofOslo,Oslo,Norway
5DepartmentofComputerScience,HelsinkiInstituteforInformationTechnologyHIIT,AaltoUniversity,Espoo,Finland
*Correspondingauthor:E-mail:pekka.marttinen@aalto.fi.
AssociateEditor:NaruyaSaitou
Abstract
Prokaryotic evolution is affected by horizontal transfer of genetic material through recombination. Inference of an
evolutionary tree of bacteria thus relies on accurate identification of the population genetic structure and
recombination-derived mosaicism. Rapidly growing databases represent a challenge for computational methods to
detectrecombinationsinbacterialgenomes.WeintroduceanovelalgorithmcalledfastGEARwhichidentifieslineages
in diverse microbial alignments, and recombinations between them and from external origins. The algorithm detects
bothrecentrecombinations(affectingafewisolates)andancestralrecombinationsbetweendetectedlineages(affecting
entire lineages), thus providing insight into recombinations affecting deep branches of the phylogenetic tree. In simu-
lations,fastGEARhadcomparablepowertodetectrecentrecombinationsandoutstandingpowertodetecttheancestral
ones,comparedwithstate-of-the-artmethods,oftenwithafractionofcomputationalcost.Wedemonstratetheutilityof
themethodbyanalyzingacollectionof616whole-genomesofarecombinogenicpathogenStreptococcuspneumoniae,
for which the method provided a high-resolution view of recombination across the genome. We examined in
detail the penicillin-binding genes across the Streptococcus genus, demonstrating previously undetected genetic
exchanges between different species at these three loci. Hence, fastGEAR can be readily applied to investigate
mosaicism in bacterial genes across multiple species. Finally, fastGEAR correctly identified many known recom-
binationhotspotsandpointedtopotentialnewones.MatlabcodeandLinux/Windowsexecutablesareavailableat
https://users.ics.aalto.fi/~pemartti/fastGEAR/ (last accessed February 6, 2017).
Key words: bacterial population genetics, recombination detection, population structure, hidden Markov models,
Streptococcuspneumoniae,antibioticresistance.
A
r
t
Introduction diversificationwasshowntobegreaterthanthatofmutation ic
alone(Vos,2009). l
e
Microbialgenomesareconstantlysubjectedtoanumberof Theprevalenceofrecombinationissuggestiveofitsimport-
evolutionary processes, including mutation, gene gain and anceformicrobialevolution,withpotentialadaptivebenefits.
loss, genetic rearrangement and recombination, the latter Geneticexchangebetweendifferentstrainshasbeenarguedto
here broadly defined as any form of horizontal transfer of play an important role in shaping of bacterial communities
DNA.Theimportanceofrecombinationinprokaryoticevo- (Polzetal.,2013;Marttinenetal.,2015;Shapiro,2016)andthe
lution has been recognized for some time (Feil and Spratt, emergence of new bacterial species (Fraseret al., 2007, 2009;
2001;DidelotandMaiden,2010;Hanage,2016)andgenomic Shapiroetal.,2012).Bacterialrecombinationhasalsoproveda
studieshavebecomeanimportantsourceofdatatomeasure powerful adaptive weapon against major forms of clinical
its contribution (Polz et al., 2013). Comparative studies of interventions: antibiotics and vaccines (Hanage et al., 2009;
prokaryotic genomes have found that the vast majority of Croucheretal.,2011;Perronetal.,2012).Asmostevolutionary
theirgeneshavebeenlaterallytransferredatleastonceinthe models(e.g.,phylogeneticorphylogenomicanalyses)assume
past (Dagan and Martin, 2007; Dagan et al., 2008), with no recombination, a good understanding of the impact of
around 20% of genes being acquired recently (Popa et al., recombinationonbacterialgenomesiscrucialforthecorrect
2011).Furthermore,whenmeasuredovershortertimescales, interpretationofanygenomicanalysis.
many bacterial species were found to recombine so fre- Currently,popularmethodsusedfordetectingrecombin-
quently that the impact of recombination on their genetic ationinbacterialgenomesincludeClonalFrame(Didelotand
(cid:2)TheAuthor2017.PublishedbyOxfordUniversityPressonbehalfoftheSocietyforMolecularBiologyandEvolution.
ThisisanOpenAccessarticledistributedunderthetermsoftheCreativeCommonsAttributionNon-CommercialLicense
(http://creativecommons.org/licenses/by-nc/4.0/),whichpermitsnon-commercialre-use,distribution,andreproductioninany Open Access
medium,providedtheoriginalworkisproperlycited.Forcommercialre-use,[email protected]
Mol.Biol.Evol.34(5):1167–1182 doi:10.1093/molbev/msx066 AdvanceAccesspublicationFebruary11,2017 1167
MBE
Mostowyetal. . doi:10.1093/molbev/msx066
Falush, 2007) and ClonalFrameML (Didelot and Wilson, whole-genomesofarecombinogenicpathogenStreptococcus
2015), Gubbins (Croucher et al., 2015), and BratNextGen pneumoniaesampledinMassachusetts,USA,inapaediatric
(Marttinenetal.,2012).Thethreeformerapproachesfollow carriagestudy(Croucheretal.,2013).ByapplyingfastGEARto
thelineofmethodsbasedonphylogenetictrees(Husmeier, theassembledpneumococcalpan-genome,weobtainahigh-
2005; Minin et al., 2005; Webb et al., 2009), and look for resolution view of the underlying mosaicism of all coding
clustersofpolymorphismsoneachbranchofaphylogenetic sequences. We demonstrate that fastGEAR can be used on
tree. On theother hand, BratNextGenuses HiddenMarkov onehandforin-depthunderstandingoftheextentofmosai-
Models(HMMs)tomodeltheoriginofchangesinthealign- cisminparticulargenesofinterest(e.g.,antibioticresistance
ment, where clonality (lack of recombinations) represents genes),andontheotherhandtocomparegenome-widethe
one origin and other origins represent foreign recombin- populationstructureandtheextentofmosaicismindifferent
ations. All these methods specialize in identifying imports bacterial genes. These analyses can be used in combination
originating in external sources, and are therefore appropri- with commonly used pan genome assembly pipelines like
atelyappliedtoasinglebacteriallineageatatime.Thus,they Roary(Pageetal.,2015),therebyimprovingourunderstand-
relyonanothermethodtoidentifytheunderlyingpopulation ingoftheimpactofrecombinationinbacterialgenomes.
structure, which limits their ability to provide insight into
New Approaches
species-wide or even inter-species patterns of exchange.
Withtherecentdevelopmentofhigh-throughputsequencing Here,wegiveageneralhigh-leveldescriptionofthemethod,
methods, which can process tens of thousands of bacterial and the details are presented in “Materials and Methods”
whole-genomes, such analyses have become increasingly section.Thealgorithmtakesasinputanalignmentofbacter-
interestingandnecessary. ialDNAsequencesandperformsthefollowingfourtasks:
Here, we present an approach to fulfill such a demand,
(1) Identifylineagesinthealignment.
whichidentifiesboththepopulationstructureofasequence
(2) Identify recent recombinations in the lineages, where
alignmentanddetectsrecombinationsbetweentheinferred
the“recentrecombinations”aredefinedasthosethat
lineages as well as from external origins. In particular, the
arepresentinasubsetofstrainsinalineage.
method first identifies clusters of sequences that represent
(3) Identify ancestral recombinations between the lin-
lineagesinthealignment.Thenitlocatesbothrecentrecom-
eages, where the “ancestral recombinations” are
binations, affecting some subset of strains in a lineage, and
defined as those that are present in all strains that
ancestraleventsbetweenlineagesthataffectallstrainsinthe
belongtothelineage.
lineage,asshowninfigure1.Ourapproachissimilartothe
(4) Testofsignificanceoftheputativerecombinations.
popular STRUCTURE software with the linkage model (Falush
et al., 2003) but with the following crucial differences: (1) it The distinction between recent and ancestral recombin-
iscomputationallyscalabletothousandsofbacterialgenomes ationsiswhethertherecombinationeventhappenedbefore
and(2)itprovidesinsightintothemosaicismofentirebac- orafterthemostrecentcommonancestorofthelineagein
terial populations by inferring also the ancestral recombin- whichitwasdetected(seefig.1).
ation events. As our method can quickly infer the genomic (1)Identifyinglineages:Toidentifylineagesinadataset,we
arrangementoflargebacterialdatasets,wecalleditfastGEAR. startbyrunningapreviouslypublishedclustering algorithm
We assess the accuracy of fastGEAR using extensive (Corander and Marttinen, 2006) included in the Bayesian
simulations, and compare it to three other state-of-the- AnalysisofPopulationStructure(BAPS)software(Corander
art methods. We then use it to analyze a dataset of 616 etal.,2003).ThisproducesCstrainclusterswhichrepresent
FIG.1.Simulationsofbacterialrecombinations.Thediagramshowstheunderlyingsimulationmethod,andherethecaseofP¼2populationsis
considered:blueandred.Populationsweresimulatedunderaclonalmodelofevolutionforagivensetofparameters(seeMethodssection).Three
typesofrecombinationswerethensimulatedusingtheclonalalignment.Ancestralrecombinations(case1)occurredbeforethemostrecent
common ancestor of bothpopulations, and thuswere present in all isolates of the recipient lineage. Intermediate recombinations (case2)
occurredsometimebetweenthetimewhenpopulationsemergedandpresenttime(t¼0),andthusweretypicallypresentinmultipleisolates.
Recentrecombinations(case3)occurredinthelastfewgenerations,andthusweretypicallypresentinfewisolates.Toclarify,ourmethod
identifieseventsoftype1asancestralrecombinations,whereasallotherrecombinations,affectinglessthananywholelineage(cases2and3),are
inferredasmultiplerecentrecombinationspresentinmultipleindividualstrains.
1168
MBE
RecombinationWithinBacterialPopulations . doi:10.1093/molbev/msx066
populationstructureamongthestrains.However,theclusters the direction of arecombination cannot beidentified using
as such are not optimal for recombination analysis for two thisapproach.Toresolvethisissue,wealwaysmarkthelin-
reasons.First,thealgorithmmayassigntwootherwiseiden- eagewithfewerstrainsastherecombinationrecipientinour
ticalsetsofsequencesintodistinctclustersduetooneofthe results. The convention may be justified by the principle of
sets experiencing a recombination event, because after the maximumparsimony,asitresultsinfewerstrainsinthedata
eventtherecombinantstrainsaremoresimilartoeachother set carrying a recombinant segment (but see “Discussion”
than to strains in the other set, resulting in a suboptimal section).
representation of the overall population structure across (4) Test of significance: The HMMs produce probabilities
the sequence (see supplementary fig. S1). Second, the algo- forsequencepositionsofhavingtheiroriginsinthedifferent
rithm may detect clusters that have diverged very recently, lineages,whichcanbeusedasameasureofstatisticalstrength
andtheclosenessofsuchclustersmayresultinaddednoisein ofthefindings.However,inourexperimentsweencountered
recombination detection. Due to these reasons, given the two kinds of false positive findings: first, recent recombin-
clustering pattern, we infer lineages using a Hidden Markov ations in strains that were outliers in the data set; second,
Model(HMM)approach(fig.2A).Inmoredetail,wecompare ancestral recombinations between lineages that were
allele frequencies for each cluster pair using a HMM, where diverged to the verge of not being considered the same by
the hidden states of the HMM represent equality of the the HMM, but not completely different either (see
frequencies at the polymorphic sites. All pairwise compari- “Discussion”onthelimitationsoftheHMMs).Toprunethese
sons are summarized as a distance matrix, which tells the falsepositivefindings,wemonitorthelocationsofSNPsbe-
proportionofthelengthofthesequenceswheretwoclusters tween the target strain and its ancestral lineage (for recent
are considered different. We apply the standard complete
recombinations) or between the two lineages (for ancestral
linkageclusteringwithcutoff0.5tothisdistancematrix, re-
recombinations)withinandbetweentheclaimedrecombin-
sultinginagroupingofclustersintoLgroups,whicharetaken
antsegments.Weapplyasimplebinomialtesttocomputea
as lineages in the data. This means that two clusters will
Bayes factor (BF; see, e.g., Bernardo and Smith, 2001), that
usually be considered as part of the same lineage if their
measureshowstronglythechangesinSNPdensitysupporta
sequences are considered similar for at least 50% of the se-
recombination, and we use a threshold BF¼1 for recent
quence length, although this is not strictly enforced by the
recombinations and BF¼10 for ancestral recombinations
completelinkagealgorithm.
for additional pruning of recombinations proposed by the
(2)Detectingrecentrecombinations:Toidentifyrecentre-
HMM analyses. These thresholds represent a compromise
combinations, we analyze each lineage by applying a HMM
between false positive rate and power to detect recombin-
approach for strains assigned to the lineage, this time with
ations.RecombinationswithBFslessthanthethresholdare
hidden states representing the origins of the different poly-
notreportedatall,andtheestimatedBFsfortheremaining
morphic sites in the strain (fig. 2B). Possible origins are the
recombinationsareincludedintheoutput.
other lineages detected in the data, as well as an unknown
origin,notrepresentedbyanystraininthedata.Thepositions Results
which are assigned to a different origin than the identified
PerformanceonSimulatedData
lineage of the strain are considered recombinations. After
TogiveanexampleoffastGEARperformance,weperformed
analyzing all strains in the lineage, the hyperparameters of
theHMMareupdated.Furtheriterationsofdetectingrecom- coalescent simulations involvingP¼3 lineages with a given
effective population size for each lineage N and the most
binationsandupdatinghyperparametersarecarriedoutuntil e
approximateconvergence.Thefinalreportedrecentrecom- recentcommonancestor(MRCA)timeT.Wethensimulated
binationsarethosesequencepositionswheretheprobability recombination events between them in three different
oftheassignedlineageislessthansomethreshold,wherewe modes: recent, intermediate,and ancestral. Finally, wecom-
have used a conservative threshold value equal to 0.05. If a pared the resulting, true population structure with the one
sequencepositionisconsideredrecombinant,thentheorigin inferred by fastGEAR in figure 3. We see that fastGEAR not
is set to be the lineage with the highest probability at this onlycorrectlyidentifiedthelineages,butalsofoundalmostall
position.Wenotethatalsothefullprobabilitydistributions recombinations.Theinferenceofrecentrecombinationswas
areavailablefromourimplementation. generally better than the inference of ancestral recombin-
(3) Detecting ancestral recombinations: To identify ances- ations. Of particular importance is that fastGEAR does
tral recombinations, we analyze all lineage pairs using the much better at predicting the direction of recombination
same approach as in step (1), such that the latent variables events for recent recombinations. Such direction is difficult,
forthedifferentsequencepositionshavetwopossiblestates, ifnotimpossible,todetermineforolder,ancestralrecombin-
eitherthelineagesarethesameordifferent,withtherecent ations.However,weseethatthepopulationgeneticstructure
recombination fragments treated as missing data. Putative wascorrectlyinferredinallthreeexamples,eveninthediffi-
ancestralrecombinationsbetweenlineagescorrespondtore- cult case of multiple, overlapping recombinations occurring
gions of the alignment where the inferred lineages are the atdifferenttimescales.
same—hence a portion of the genome in isolates that are TosystematicallyassesstheperformanceoffastGEARwe
overallassignedtodifferentlineages,maybeconsideredtobe performed three different sets of in silico experiments. First,
partofthesamelineage.However,itisimportanttonotethat weexaminedhowwellfastGEARdetectsrecombinationsfor
1169
MBE
Mostowyetal. . doi:10.1093/molbev/msx066
FIG.2.HiddenMarkovmodelstodetectrecombination.(A)HiddenMarkovmodelusedforidentifyinglineagesandinferringancestralrecom-
binations.Eachcolumnrepresentsapolymorphicsiteinthealignmentandrowsrepresentstrains.Theobservedstatesofthechainaretheallele
frequencieswithineachcluster(inthecaseofidentifyinglineages)orlineage(inthecaseofidentifyingancestralrecombinations).Thelatentstates
ofthechainrepresentidentityofallelefrequenciesinthetwolineagesatthepolymorphicsites.(B)HiddenMarkovmodelusedforidentifying
recentrecombinations.Theobservedstatesarenucleotidevaluesobservedinthetargetstrainandthelatentstatesarepossibleoriginsofthe
nucleotides.Thepossibleoriginsincludeallobservedlineagesplusanunknownorigin.
different population parameters. To this end, we varied the dependsonthestrengthoftheunderlyingpopulationgenetic
within-population distance (achieved by changing N) and structure.
e
the between-population distance (achieved by changing T). Inthesecondsetofexperiments,wecomparedfastGEAR
TheresultsareshowninsupplementaryfigureS2.Weseethat to other recombination-detection methods: STRUCTURE (link-
fastGEARgenerallydetectedrecombinationswell,particularly age model), Gubbins, and ClonalFrameML. To facilitate the
therecentonesastheysharehigherresemblancetotheorigin comparisonwiththelattertwophylogeny-basedapproaches,
andarethusbydefinitioneasiertodetect.Thefalse-positive designedtodetectrecombinationsfromexternaloriginsina
ratewaslowforalltypesofrecombinationsdetectedanddid lineage-by-lineage manner, and to explore the ability of
notvarywiththebetween-orthewithin-lineagedistance.By fastGEAR to detect recombinations from unknown origins,
contrast,weobservedthattheproportionofdetectedrecom- we ran fastGEAR both using the entire alignment and with
binationswashighlydependentonthebetween-lineagedis- each lineage individually. Figure 4 shows results in terms of
tance. This is because in the absence of clear population falsepositiverate(thenumberof incorrectrecombinations,
genetic structure, populations are relatively closely related left column), sensitivity (the proportion of true recombin-
andtherearetoofewpolymorphismstosignalthepresence ationsdetected,middle),andtheproportionoftotallength
ofarecombination.Furthermore,ahigherwithin-lineagedis- oftruerecombinationscoveredbydetectedrecombinations
tance often affected the inference of ancestral recombin- (right).
ations as it generated the intra-lineage population genetic Weseethatthefalsepositiveratestayedverylowwithall
structure. Thus, as expected, performance of fastGEAR methods. Furthermore, overall, fastGEAR had similar
1170
MBE
RecombinationWithinBacterialPopulations . doi:10.1093/molbev/msx066
FIG.3.Visualassessmentoftheinferredpopulationgeneticstructure.Thefigureshowsthepopulationgeneticstructureofthesimulateddata.In
eachpanel,therowscorrespondtosequences,columnscorrespondtopositionsinthealignmentandcolorsshowdifferentpopulations.Theleft
columnshowsthesimulated,truestructurealthoughtherightcolumnshowsthepopulationgeneticstructureinferredbyfastGEAR.Theorderof
thesequencesinbothcolumnsisidentical,andthecolorsareassignedrandomly,thuspopulationsareinthesameorder(1,2,3)butcanbeof
differentcolorontheleftandontheright.Threefigurerowscorrespondtothreedifferentsimulationscenarios:onlyrecentrecombinations(top),
onlyancestralrecombinations(middle),andallthreetypesofrecombinations(bottom).Thefollowingparameterswereusedinthesimulations:
P¼3,n¼20,N ¼50,T¼2e4; l¼2e(cid:2)6,L¼10kb(allrows);C ¼800andR ¼5,R ¼R ¼0(toppanel);C ¼800andR ¼3,R ¼R ¼0
e r r i a a a r i
(middlepanel);C ¼500; C ¼500andR ¼3,R¼4,andR ¼6(bottompanel).
a r a i r
sensitivitytodetectrecentandintermediaterecombinations generallyquitesimilar,exceptthatSTRUCTUREexplorestheen-
totheothermethods,andnomethodwassystematicallythe tireparameterspaceusingMCMCwhereasfastGEARusesthe
best. Particularly interesting is the comparison between mostprobableclusteringandpointestimatesforhyperpara-
STRUCTURE and fastGEAR, because the two methods are meters.Ingeneral,STRUCTUREhadabithighersensitivity,butat
1171
MBE
Mostowyetal. . doi:10.1093/molbev/msx066
FIG.4. ComparisonoffastGEARandotherrecombinationdetectionmethods.ThefigureshowstheperformanceoffastGEARcomparedwith
othermethods:STRUCTURE,GubbinsandClonalFrameML.Toprowshowsresultsforrecent,middlerowforintermediate,andbottomrowfor
ancestralsimulatedrecombinations.BothrecentandintermediatesimulatedrecombinationsweredetectedbyfastGEARinthesamewayas
“recent”recombinations.Theleftcolumnshowsthefalsedetectionrate,namely,themeannumberoffalse-positiverecentrecombinationsper
strain(top/middle)andancestralrecombinationsperalignment(bottom).Themiddlecolumnshowstheproportionoftruerecombinations
detected,andtherightcolumnshowstheproportionofthetotalrecombinationlengthdetected.Horizontalaxisshowsthebetween-population
distanceper100bp(simulatedbyvaryingTbetween1.0e3and2.0e4).Differentlinesshowperformanceofdifferentapproaches.STRUCTUREwasrun
for 400,000 generations (200,000 burn-in), with true populations set as prior and with three independent chains to test for convergence.
ClonalFrameMLwasconditionedonthetruetreetopology.MagentalineshowsresultsoffastGEARrunonthefullalignment,andredlinefor
fastGEARrunlineage-by-lineage.Eachpointrepresentstheaverageoftensimulations.Thefollowingparameterswereusedinthesimulations:
P¼3,n¼30,l¼2:0e–6,L¼20kb,R ¼300; R ¼600andR ¼5,R ¼R ¼5.
r a r i a
amuchhighercomputationalcost(between1.5and2orders between-species bacterial data. An additional advantage of
of magnitude; see supplementary fig. S3). Furthermore, fastGEARisthat itcanhandlemissing datainastraightfor-
STRUCTURE was here run conditioning on the true value of ward manner, although, as expected, a lot of missing data
thenumberofpopulationsKandthetruelineagemember- may reduce the accuracy of recombination detection (see
shipsofdifferentstrainsaspriors,whereasfastGEARhadno supplementaryfig.S5).
suchknowledge. Tofacilitatethecomparisonofdifferentmethods,thepre-
When run in a lineage-by-lineage manner, fastGEAR had viousbenchmarksimulationshadpredefinedrecombination
clearlylowersensitivitythanwiththefullalignment,andthe probabilitiesthatwerenotuniformintime,thusnotallowing
lengths of the external recombinations were often overesti- to investigate how the number of detected recombinations
mated(seesupplementaryfig.S4).ThisisbecausefastGEAR, depends on some assumed fixed underlying recombination
similarly to STRUCTURE, gains its statistical power from having rate.Toinvestigatethis,wecreatedathirdsetofexperiments,
theactualoriginsofrecombinations.Importantly,throughout wherewesimulatedgenealignmentsfromasinglebacterial
simulations fastGEAR detected ancestral recombinations lineagewithpredefinedinternalandexternalrecombination
equally well to recent and intermediate recombinations. ratesusingancestralrecombinationgraph(ARG)simulations
Thisisparticularlyencouragingasnoneoftheothermethods (Brown et al., 2016). Supplementary figures S6 and S7 show
could detect ancestral recombinations. These results there- theresults for different underlyingmutationrates andgene
foreshowthatalthoughfastGEARisabletodetectexternal sizes. Again, we see that the diversity levels affected recom-
recombinations in a lineage-by-lineage manner, it is at its binationdetectionsensitivity,withacleardetectionlimitfor
strongest when applied to investigate between-lineage or extremely high recombination rates which reduced the
1172
MBE
RecombinationWithinBacterialPopulations . doi:10.1093/molbev/msx066
linkagebetweenSNPstozero.Ontheotherhand,wedidnot S9 shows results separately for recent and ancestral recom-
observe a significant effect of gene length on the detection binations,demonstratingthatmuchofthestructureremains
rate: given a recombination rate, the number of detected hiddenifancestralrecombinationsarenotconsidered.
recombinationsperkbremainedapproximatelyconstant,ex- PreviouslyasoftwarecalledBratNextGen(Marttinenetal.,
cept for very short genes, where the length of a recombin- 2012),developedbythesameauthors,hasbeenusedtoin-
ation may exceedthelengthof thegene, inwhich casethe vestigatethepopulationstructureacrossdiversestrainsatthe
strains affected may be classified into a different lineage, ra- sameloci(Croucheretal.,2013;Jensenetal.,2015).Asmen-
therthanbeidentifiedasrecombinant. tioned previously, BratNextGen was designed for detecting
importsfromexternaloriginsinwhole-genomealignmentsof
AnalysisofStreptococcuspneumoniaeData closelyrelatedstrains.Forthisreason,itdoesnotautomatic-
WenextdemonstratetheutilityoffastGEARinthreediffer- allyidentifytheoriginoftheimports,therebyrequiringpost-
ent data analyses using a whole-genome collection of 616 processinglikemanualcurationusedbyJensenetal.(2015).
isolates of S. pneumoniae. First, we analyzed alignments of To better understand how the difference between the two
threeb-lactamgenesresponsibleforpenicillinresistance,aug- methodsaffectstheinferenceoftheunderlyinggeneticmo-
menting thealignmentswithsequences fromotherspecies. saicism, we compared the results using both approaches
This demonstrates the ability of fastGEAR to produce a de- using simulated data (see supplementary fig. S10) and the
tailedpopulationstructureofgenomeregionsofinterestand pbp’s (see supplementary fig. S11). The results show that
detectrecombinationbetweenspecies.Second,weanalyzed the lineages detected by BratNextGen as well as the strains
the pneumococcal pan-genome by running fastGEAR inde- identified as recombinant are, broadly speaking, similar to
pendentlyonthe2,113codonalignmentsforCOGswithat thosefromfastGEAR.However,asexpected,theBratNextGen
least 50 sequences from the previously-described pneumo- outputmaybearbitraryifmultiplediverselineagesareincluded
coccal population, to detect recombination hotspots in the intotheanalysis(seesupplementaryfig.S10).Thiscanbebest
pneumococcal genome. Third, we estimated a high- seen when comparing the results of the pbp analysis by
resolution genome-wide population structure of the fastGEAR and BratNextGen side-by-side, zooming in on
pneumococcus population, by summarizing the analyses of theresultsacrossmultiplespecies(seefig.6;supplementary
individualCOGsintermsoftheproportionofgenome-wide fig. S12 shows results separately for recent and ancestral
sharedancestrybetweenisolatepairs.Thethirdanalysisalso recombinations).OnecanseethatfastGEARidentifiedmo-
demonstrates how the pneumococcal species-wide phyl- saicism between species much better than BratNextGen.
ogeny,estimatedfromacorealignment,emergesasanaver- Overall, we conclude that fastGEAR is much more suited
ageofhighlyvariablepopulationstructuresatindividualloci, fortheanalysisofmosaicismofdiversebacterialgenesthan
rather than representing clonal descent at any genome other,single-lineagemethodslikeBratNextGen.
region. Wealsoinvestigatedtherelationbetweenrecombination
andthelevelsofantibioticresistanceforS.pneumoniaeiso-
Inter-SpeciesRecombinationatPenicillin-Binding lates analyzed here (see supplementary fig. S13). The
Proteins Spearmanrankcorrelationbetweenthedetectedrecombin-
The population structure of b-lactam genes, which encode ationsandtheminimuminhibitoryconcentration(MIC)was
for penicillin binding proteins, has been previously analyzed found significant in all genes (pbp1a, q¼0:66; p(cid:3)0;
to understand the genetic basis of penicillin resistance and pbp2b, q¼0:72; p(cid:3)0; pbp2x, q¼0:75; p(cid:3)0). Here
investigate the flow of antibiotic resistance determinants thenumberofbreakpointsduetobothrecentandancestral
among streptococci (Croucher et al., 2013; Jensen et al., recombinationswasusedasthevariabletestedagainstMIC,
2015).Here,weincludedthewholepopulationof616S.pneu- but all these correlations remained highly significant even if
moniae strains, together with sequences from other species usingeitherrecentorancestralrecombinationalone,orwhen
available in databases (list of accession numbers available in representing the amount of recombination by the total
supplementarytableS1),resultinginthemostcomprehensive lengthofrecombinantsequence.
collectionanalyzedforthispurpose.Weanalyzedeachofthe
b-lactamgenes,pbp1a,pbp2b,andpbp2x,usingfastGEAR.We ComparisonofRecombinationLevelsacrossDifferent
alsoinvestigatedthesensitivityoftheresultsonthepresence Proteins
ofdiversestrainsbyremovingtheoutlierstrains(i.e.,thosethat Wenextcomparedthelevelsofrecentandancestralrecom-
belongedtolineageswithfiveorfewerisolatesintheoriginal bination between different proteins. Consistent with a con-
analysis). stant rate of recombination over the history of this
Results with outliers removed are presented in figure 5, population, both measures were significantly correlated
andtheyshowahigh-resolutionviewofthepbpgenesacross (R2 ¼0:46; p(cid:3)0)withthemeannumberofrecentrecom-
streptococci, with multiple recombination events between binations almost twice the number of ancestral recombin-
different species. Supplementary figure S8 shows results be- ations (1.4 vs.2.7;seesupplementary fig.S14).Basedonthe
foreremovingoutliers,andweseethattheinclusionofhighly simulations, the overall diversity in a gene may affect the
divergent outlier strains slightly lowers the sensitivity to de- number of recombinations detected (see fig. 4 and supple
tect recombinations, however the overall picture remains mentary fig.S6). To investigatethepotential bias causedby
qualitatively identical. By comparison, supplementary figure this, we computed four common diversity measures and
1173
MBE
Mostowyetal. . doi:10.1093/molbev/msx066
FIG.5.Resultsshowinginter-speciesrecombinationatpenicillin-bindingproteins.ThefigureshowsfastGEARresultsforcombineddatasetswith
the616S.pneumoniaestrainsandsomenumberofadditionalsequencesfromotherspecies(104inpbp1a,129inpbp2b,127inpbp2x).The
phylogenyandthesequenceclusters(SCs)ontheleftshowthecore-genome-basedtreewith15majormonophyleticclustersfortheS.pneumoniae
strains.StrainsfromotherspeciesareshownontopoftheS.pneumoniaestrains.Thespeciesannotationisrepresentedbycolorsontheleftsideofthe
additionalstrains,abovethephylogeny.NotethatthecolorsusedtoannotatespeciesareindependentofthecolorsinthefastGEARoutputplots,
wherecolorsrepresentlineagesdetectedbyfastGEAR,exceptwhitewhichdenotesgapsinthealignment.
compared them with the estimated recombination counts hit for both recent and ancestral recombinations (Kadioglu
(see supplementary fig. S15). No obvious relationship was et al., 2008). Another top hit was the first of the rhamnose
seen, suggesting that results mainly represent true variation genes(rmlA),whichoftenservesasabreakpointinserotype
inrecombinationintensities. switching events (see Supplementary text and fig. S16 for a
Among the genes with the highest number of ancestral detaileddiscussion).Wealsofoundahighnumberofrecent
andrecentrecombinationsshownintables1and2,respect- recombinations in the zinc metalloprotease zmpA, which
ively,wefoundmanylocipreviouslyidentifiedasrecombin- cleaves human immunoglobulin A1 (Weiser et al., 2003).
ation hotspots (Supplementary tables S2 and S3 show top Thethirdgrouparegenesinvolvedindeterminingresistance
hits for recombination intensities normalized by alignment to antibiotics, including sulphamethoxazole resistance (folC,
length). Theseproteinscanbeclassifiedintoseveralgroups. see supplementary fig. S17), as well as b-lactams (pbp1a,
The first group consists of mobile genetic elements, which pbp2b and pbp2x) discussed above. With the exception of
include integrative and conjugative elements, prophages, themosaiczmpAsequences,theseproteinswerepreviously
phage-relatedchromosomalislands,andinsertionsequences. identifiedasrecombinationhotspotsingloballydisseminated
Thisisnotsurprisingasfrequentbetween-andwithin-lineage lineages(Croucheretal.,2011,2014b).
recombinationofmobilegeneticelementshasbeenreported Wealsofoundhighlyrecombinogenicproteinswhichhave
previously (Croucher et al., 2014a). The second group are not been previously identified as recombination hotspots.
proteins which are engaged in the interactions with the OneexampleisthechromosomepartitioningSMCprotein,
host.PneumococcalsurfaceproteinC(pspC),whichplaysa whichfunctionsinchromosomalsegregationduringcelldiv-
centralroleinpathogenesisofthepneumococcus,wasatop ision(Brittonetal.,1998)andisoneofthetophitsinboth
1174
MBE
RecombinationWithinBacterialPopulations . doi:10.1093/molbev/msx066
FIG.6.ComparisonoffastGEARandBratNextGenThefigurepresentsadetailedcomparisonoffastGEARandBratNextGen.Theresultsareshown
forthepbp2xgene,zoomingintothesequencesfrommultipledifferentspeciesappearingontopinfigure5.WeseethatfastGEARisabletodetect
mosaicstructurebetweenspecies.
Table1. SummaryofCOGswiththeHighestNumberofAncestralRecombinationEvents(>20)IdentifiedbyfastGEAR.
SNPs Sequences Clusters Lineages Ancestral GeneName
788 288 17 10 58 PneumococcalsurfaceproteinC
818 2149 24 17 51 InsertionsequenceIS1167
543 616 21 8 47 ChromosomepartitionproteinSMC
618 2320 25 19 39 InsertionsequenceIS1381
297 325 14 11 34 BacteriophageDNAbindingprotein
1177 378 10 7 32 Tn5252relaxase
720 616 12 7 31 Penicillin-bindingprotein2x
191 616 20 10 30 RibosomalproteinL11methyltranferase
332 612 14 7 23 Twocomponentsystemhistidinekinase
258 336 11 7 23 Capsulelocusglucose-1-phosphatethymidyltransferaseRmlA
320 616 22 9 22 Phenylalanyl-tRNAsynthetase(PheS)
310 229 11 7 22 Proteinfoundinphage-relatedchromosomalislands
recentandancestralrecombinations.Othergenesthatwere ancestry (PSA) matrix, which summarizes the fastGEAR re-
also inferred to undergo high levels of recombination are sultsforall2,113COGalignments.Specifically,foreachpairof
phenylalanyl- and valyl-tRNA synthetases, enzymes that at- isolates we analyzed the population structure at each COG
tachtheaminoacidsphenylalanineandvalinetotheircog- with putative lineages and the list of recent and ancestral
natetRNAmoleculesduringthetranslationprocess.Previous recombinationsdetectedintheCOG.IftheCOGwaspresent
reportsshowthatrecombinationandhorizontalgenetrans- inbothisolates,wecomputedtheproportionofthelengthof
fer frequently occur in aminoacyl-tRNA synthetases (aaRS) theCOGsequenceinwhichthetwoisolateswereassignedto
(Woeseetal.,2000).ThehorizontalacquisitionofaaRSvari- thesamelineage.IftheCOGwasabsentineitherorbothof
ants may be implicated in resistance to antibiotics (Woese theisolates,theywerenotcomparedatthisCOG;ifmultiple
etal.,2000),andatleastoneatypicaladditionalaaRShasbeen copiesofthegenewerefound,thenallpossiblecomparisons
found on Pneumococcal Pathogenicity Island 1 (Croucher betweenthetwoisolateswereincluded,andcorrespondingly
et al., 2009). Although within-species recombination of takenintoaccountinthetotallengthofsequencecompared.
aaRS has not been widely investigated, our results suggest TheresultingPSAmatrixtogetherwithapreviouslypub-
that this process plays an important role in the evolution lished core-gene-based phylogeny and 15 monophyletic se-
ofpneumococci. quence clusters (SCs), which can be taken as lineages, is
showninfigure7.Overall,thePSAresultsarehighlyconcord-
High-ResolutionViewofPopulationStructure antwiththetreeandtheSCs.First,strainswithinSCsshare
Toinvestigatethepopulation structureoftheentirecollec- almostalloftheirancestry,suchthattheaveragePSAwithin
tion of isolates, we calculated the proportion of shared differentSCsrangesfrom85%upto98%,whichisvisibleas
1175