Early Divergent Strains of Yersinia pestis in Eurasia 5000 Years Ago

45 Pages·2015·7.71 MB·English

by  Simon Rasmussen

Article
Yersinia pestis
Early Divergent Strains of in Eurasia
5,000 Years Ago
Graphical Abstract Authors
SimonRasmussen,MortenErikAllentoft,
KasperNielsen,...,RasmusNielsen,
KristianKristiansen,EskeWillerslev
Correspondence
[email protected]
In Brief
Theplague-causingbacteriaYersinia
pestisinfectedhumansinBronzeAge
Eurasia,threemillenniaearlierthanany
historicalrecordsofplague,butonly
acquiredthegeneticchangesmakingita
highlyvirulent,flea-bornebubonicstrain
(cid:1)3,000yearsago.
Highlights
d YersiniapestiswascommonacrossEurasiaintheBronze
Age
d ThemostrecentcommonancestorofallY.pestiswas5,783
yearsago
d Theymtgenewasacquiredbefore951calBC,givingriseto
transmissionviafleas
d BronzeAgeY.pestiswasnotcapableofcausingbubonic
plague
Rasmussenetal.,2015,Cell163,571–582
October22,2015ª2015TheAuthors
http://dx.doi.org/10.1016/j.cell.2015.10.009
Article
Yersinia pestis
Early Divergent Strains of
in Eurasia 5,000 Years Ago
SimonRasmussen,1,18MortenErikAllentoft,2,18KasperNielsen,1LudovicOrlando,2MartinSikora,2Karl-Go¨ranSjo¨gren,3
AndersGormPedersen,1MikkelSchubert,2AlexVanDam,1ChristianMoliinOutzenKapel,4HenrikBjørnNielsen,1
SørenBrunak,1,5PavelAvetisyan,6AndreyEpimakhov,7MikhailViktorovichKhalyapin,8ArtakGnuni,9AivarKriiska,10
IrenaLasak,11MaitMetspalu,12VyacheslavMoiseyev,13AndreiGromov,13DaliaPokutta,3LehtiSaag,12LiiviVarul,10
LevonYepiskoposyan,14ThomasSicheritz-Ponte´n,1RobertA.Foley,15MartaMirazo´nLahr,15RasmusNielsen,16
KristianKristiansen,3andEskeWillerslev2,17,*
1CenterforBiologicalSequenceAnalysis,DepartmentofSystemsBiology,TechnicalUniversityofDenmark,Kemitorvet,Building208,
2800KongensLyngby,Denmark
2CentreforGeoGenetics,NaturalHistoryMuseumofDenmark,UniversityofCopenhagen,ØsterVoldgade5–7,1350Copenhagen,Denmark
3DepartmentofHistoricalStudies,UniversityofGothenburg,40530Gothenburg,Sweden
4SectionforOrganismalBiology,DepartmentofPlantandEnvironmentalSciences,UniversityofCopenhagen,Thorvaldsensvej40,
1871FrederiksbergC,Denmark
5NovoNordiskFoundationCenterforProteinResearch,UniversityofCopenhagen,2200Copenhagen,Denmark
6DivisionofArmenologyandSocialSciences,InstituteofArchaeologyandEthnography,NationalAcademyofSciences,0025Yerevan,
RepublicofArmenia
7InstituteofHistoryandArchaeologyRAS(SouthUralDepartment),SouthUralStateUniversity,454080Chelyabinsk,Russia
8OrenburgMuseumofFineArts,460000Orenburg,Russia
9DepartmentofArchaeologyandEthnography,YerevanStateUniversity,0025Yerevan,RepublicofArmenia
10DepartmentofArchaeology,UniversityofTartu,51003Tartu,Estonia
11InstituteofArchaeology,UniversityofWroc1aw,50-139Wroc1aw,Poland
12DepartmentofEvolutionaryBiology,EstonianBiocentreandUniversityofTartu,51010Tartu,Estonia
13PetertheGreatMuseumofAnthropologyandEthnography(Kunstkamera)RAS,199034St.Petersburg,Russia
14LaboratoryofEthnogenomics,InstituteofMolecularBiology,NationalAcademyofSciences,0014Yerevan,Armenia
15LeverhulmeCentreforHumanEvolutionaryStudies,DepartmentofArchaeologyandAnthropology,UniversityofCambridge,Cambridge
CB21QH,UK
16CenterforTheoreticalEvolutionaryGenetics,UniversityofCalifornia,Berkeley,California94720-3140,USA
17DepartmentofZoology,UniversityofCambridge,DowningStreet,CambridgeCB23EJ,UK
18Co-firstauthor
*Correspondence:[email protected]
http://dx.doi.org/10.1016/j.cell.2015.10.009
ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/).
SUMMARY INTRODUCTION
The bacteria Yersinia pestis is the etiological agent Plague is caused by the bacteria Yersinia pestis and is being
of plague and has caused human pandemics with directly transmitted through human-to-human contact (pneu-
millions of deaths in historic times. How and monicplague)orviafleasasacommonvector(bubonicorsepti-
when it originated remains contentious. Here, we cemic plague) (Treille and Yersin, 1894). Three historic human
report the oldest direct evidence of Yersinia pestis plaguepandemicshavebeendocumented:(1)theFirstPandemic,
which started with the Plague of Justinian (541–544 AD), but
identified by ancient DNA in human teeth from Asia
continuedintermittentlyuntil(cid:1)750AD;(2)theSecondPandemic,
and Europe dating from 2,800 to 5,000 years ago.
whichbeganwiththeBlackDeathinEurope(1347–1351AD)and
By sequencing the genomes, we find that these
includedsuccessivewaves,suchastheGreatPlague(1665–1666
ancient plague strains are basal to all known AD),untilthe18thcentury;(3)theThirdPandemic,whichemerged
Yersinia pestis. We find the origins of the Yersinia
inChinainthe1850sanderuptedthereinamajorepidemicin1894
pestis lineage to be at least two times older than beforespreadingacrosstheworldasaseriesofepidemicsuntil
previous estimates. We also identify a temporal themiddleofthe20thcentury(Bosetal.,2011;Cuietal.,2013;
sequenceofgeneticchangesthatleadtoincreased Drancourtetal.,1998;Harbecketal.,2013;Parkhilletal.,2001;
virulence and the emergence of the bubonic PerryandFetherston,1997;Wagneretal.,2014).Earlieroutbreaks
plague. Our results show that plague infection such as the Plague of Athens (430–427 BC) and the Antonine
was endemic in the human populations of Eurasia Plague (165–180 AD) may also have occurred, but there is no
directevidencethatallowsconfidentattributiontoY.pestis(Dran-
atleast3,000yearsbeforeanyhistoricalrecordings
courtandRaoult,2002;McNeill,1976).
of pandemics.
Cell163,571–582,October22,2015ª2015TheAuthors 571
Figure 1. Archaeological Sites of Bronze
A B
AgeYersiniapestis
(A)MapofEurasiaindicatingtheposition,radio-
carbondatedagesandassociatedculturesofthe
samplesinwhichY.pestiswereidentified.Dates
aregivenas95%confidenceintervalcalendarBC
years.IA:IronAge.
(B) Burial four from Bulanovo site. Picture by
MikhailV.Khalyapin.SeealsoTableS1.
assemble the Y. pestis genomes to an
averagedepthof0.14–29.5X,with12%–
95% of the positions in the genome
coveredatleastonce(Table1,TableS2,
S3, and S4). We also recovered the
sequences of the three plasmids pCD1,
pMT1, and pPCP1 (0.12 to 50.3X in
average depth) the latter two of which
Theconsequencesoftheplaguepandemicshavebeenwell- are crucial for distinguishing Y. pestis from its highly similar
documentedandthedemographicimpactsweredramatic(Little ancestor Y. pseudotuberculosis (Table 1, Figure 2, Table S3)
etal.,2007).TheBlackDeathaloneisestimatedtohavekilled (Bercovieretal.,1980;Chainetal.,2004;Parkhilletal.,2001).
30%–50%oftheEuropeanpopulation.Economicandpolitical ThehostindividualsfromwhichY.pestiswasrecoveredbelong
collapses have also been in part attributed to the devastating to Eurasian Late Neolithic and Bronze Age cultures (Allentoft
effects of the plague. The Plague of Justinian is thought to etal.,2015),representedbytheAfanasievocultureinAltai,Sibe-
have played a major role in weakening the Byzantine Empire, ria(2782calBC,2794calBC,n=2),theCordedWareculturein
andtheearlierputativeplagueshavebeenassociatedwiththe Estonia(2462calBC,n=1),theSintashtacultureinRussia(2163
declineofClassicalGreeceandlikelyunderminedthestrength calBC,n=1),theUneticecultureinPoland(2029calBC,n=1),
oftheRomanarmy. theAndronovocultureinAltai,Siberia(1686calBC,n=1),and
MolecularclockestimateshavesuggestedthatY.pestisdiver- an early Iron Age individual from Armenia (951 cal BC, n = 1)
sifiedfromthemoreprevalentandenvironmentalstress-tolerant, (TableS1).
butlesspathogenic,entericbacteriumY.pseudotuberculosisbe-
tween2,600and28,000yearsago(Achtmanetal.,1999,2004; AuthenticationofYersiniapestisAncientDNA
Cui et al., 2013; Wagner et al., 2014). However, humans may Besidesapplyingstandardprecautionsforworkingwithancient
potentiallyhavebeenexposedtoY.pestisformuchlongerthan DNA (Willerslev and Cooper, 2005), the authenticity of our
thehistoricalrecordsuggests,thoughdirectmolecularevidence findings are supported by the following observations: (1) The
forY.pestishasnotbeenobtainedfromskeletalmaterialolder Y. pestis sequences were identified in significant amounts in
than 1,500 years (Bos et al., 2011; Wagner et al., 2014). The shotgun data from eight of 101 samples, showing that this
most basal strains of Y. pestis (0.PE7 clade) recorded to date finding is not due to a ubiquitous contaminant in our lab or in
wereisolatedfromtheQinghai-TibetPlateauinChinain1961– thereagents.Indeed,furtheranalysisshowedthatoneofthese
1962(Cuietal.,2013). eightwasmostlikelynotY.pestis.Wealsosequencedallnega-
WeinvestigatedtheoriginofY.pestisbysequencingancient tiveDNAextractioncontrolsandfoundnosignsofY.pestisDNA
bacterialgenomesfromtheteethofBronzeAgehumansacross in these (Table S3). (2) Consistent with an ancient origin, the
Europe and Asia. Our findings suggest that the virulent, flea- Y. pestis reads were highly fragmented, with average read
borne Y. pestis strain that caused the historic bubonic plague lengthsof43–65bp(TableS3)andalsodisplayedclearsignsof
pandemicsevolvedfromalesspathogenicY.pestislineagein- C-T deamination damage at the 50 termini typical of ancient
fecting human populations long before recorded evidence of DNA(Figure3,FigureS1).Becausetheplasmidsarecentralfor
plagueoutbreaks. discriminating between Y. pestis and Y. pseudotuberculosis,
wetestedseparatelyforDNAdamagepatternsforthechromo-
RESULTS someandforeachoftheplasmids.Forthesevensamples,we
observesimilarpatternsofDNAdamageforchromosomeand
IdentificationofYersiniapestisinBronzeAgeEurasian plasmidsequences(Figure3,FigureS1).(3)Weobservecorre-
Individuals lated DNA degradation patterns when comparing DNA degra-
We screened c. 89 billion raw DNA sequence reads obtained dation in the Y. pestis sequences and the human sequences
fromteethof101BronzeAgeindividualsfromEuropeandAsia fromthehostindividual.GiventhatDNAdecaycanbedescribed
(Allentoftetal.,2015)andfoundthatsevenindividualscarriedse- asarateprocess(Allentoftetal.,2012),thissuggeststhatthe
quencesresemblingY.pestis(Figure1,TableS1,Supplemental DNA molecules of the pathogen and the human host have a
Experimental Procedures). Further sequencing allowed us to similar age (Figure 3, Figure S1, Table S3 and Supplemental
572 Cell163,571–582,October22,2015ª2015TheAuthors
Table1. OverviewoftheY.pestisContainingSamples
Sample Country Site Culture Date(calBC) CO92 pMT1 pPCP1 pCD1
RISE00 Estonia Sope CordedWare 2575–2349 0.39 0.36 1.40 0.66
RISE139 Poland Chociwel Unetice 2135–1923 0.14 0.24 0.76 0.28
RISE386 Russia Bulanovo Sintashta 2280–2047 0.82 0.96 1.12 1.60
RISE397 Armenia Kapan EIA 1048–885 0.25 0.40 6.88 0.50
RISE505 Russia Kytmanovo Andronovo 1746–1626 8.73 9.15 34.09 17.46
RISE509 Russia AfanasievoGora Afanasievo 2887–2677 29.45 16.96 31.22 50.32
RISE511 Russia AfanasievoGora Afanasievo 2909–2679 0.20 0.24 1.19 0.60
ThedatingisdirectAMSdatingofbonesandteethandisgivenas95%confidenceintervalcalendarBCyears(detailsaregiveninTableS1).The
columnsCO92,pMT1,pPCP1andpCD1correspondtosequencingdepth.Additionalinformationonthearchaeologicalsitesandmappingstatistics
canbefoundintheSupplementalExperimentalProceduresandTableS1,S2,andS3.EIA:EarlyIronAge,AMS:AcceleratorMassSpectrometry.
Experimental Procedures). (4) Because of the high sequence 2013;Morellietal.,2010;Wagneretal.,2014),whichwerebased
similarity between Y. pestis and Y. pseudotuberculosis, we onlessthan2,300nucleotidesthatwereascertainedtobevari-
mapped all reads both to the Y. pestis CO92 and to the ableinY.pestis,likelyleadingtolowerstatisticalaccuracythan
Y. pseudotuberculosis IP32953 reference genomes (Chain with whole-genome analyses. Furthermore, the use of SNPs
etal.,2004).ConsistentwithbeingY.pestis,theseveninvesti- ascertained to be variable in Y. pestis would downwardly bias
gated samples displayed more reads matching perfectly (edit estimatesofbranchlengthsinY.pseudotuberculosisandlead
distance=0)towardY.pestis(Figure3,FigureS2).Onesample tounderestimatesoftheY.pestisversusYpseudotuberculosis
(RISE392)wasmostlikelynotY.pestisbasedonthiscriterion. divergence time, as seen in the branch length of the Y. pestis
(5)AnaiveBayesianclassifiertrainedonknowngenomespre- cladetoY.pseudotuberculosis(FigureS3).Thetopologyofour
dictsthesevensamplestobeY.pestiswith100%posteriorprob- whole genome tree shows Y. pestis as a monophyletic group
ability, while RISE392 is predicted to have 0% probability of within Y. pseudotuberculosis with RISE505 and RISE509 (Fig-
being Y. pestis (Figure S2, Table S3). (6) If the DNA was from ure 4A, black arrow, Figure S4) clustered together within the
other organisms than Y. pestis, we would expect the reads to Y.pestisclade.TheY.pestissub-treetopology(Figure4B,Fig-
be more frequently associated with either highly conserved or ure S4) is similar to that reported previously (Bos et al., 2011;
low-complexityregions.However,wefindthereadstobedistrib- Cui et al., 2013; Morelli et al., 2010; Wagner et al., 2014), but
utedacrosstheentiregenome(FigureS2),andcomparisonof withthetwoancientstrains(RISE505andRISE509)fallingbasal
actual coverage versus the coverage that would be expected toallotherknownstrainsofY.pestis(100%bootstrapsupport).
fromreadlengthdistributionsandmappabilityofthereference
sequencesarealsoinagreementforthesevensamples(Figure3). DeterminationofYersiniapestisDivergenceDates
(7)Inamaximumlikelihoodphylogeny,therecoveredY.pestis To determine the dates for the mostrecent common ancestor
genomicsequencesofRISE505andRISE509areclearlywithin (MRCA)ofY.pestisandY.pseudotuberculosis,andforallknown
theY.pestiscladeandbasaltoallcontemporaryY.pestisstrains Y.pestisstrains,weusedaBayesianMarkovChainMonteCarlo
(Figure4)(seebelow). approachimplementedinBEAST2(Bouckaertetal.,2014)ona
subsetofthesupermatrix.WeestimatedtheMRCAofY.pestis
ThePhylogeneticPositionoftheBronzeAgeYersinia and Y. pseudotuberculosis to be 54,735 years ago (95% HPD
pestisStrains [highest posterior density] interval: 34,659–78,803 years ago)
To determine the phylogenetic positions of the two high (Figure 4C, Figure S5, Table S5), which is about twice as old
coverageancientY.pestisstrains,RISE505(Andronovoculture compared to previous estimates of 2,600–28,000 years ago
1686 cal BC, 8.7X) and RISE509 (Afanasievo culture, 2746 cal (Achtman et al., 1999, 2004; Cui et al., 2013; Wagner et al.,
BC, 29.7X), we mapped the reads, together with reads from 2014). Additionally, we estimated the age of the MRCA of all
strains of Yersinia similis (n = 5), Y. pseudotuberculosis (n = knownY.pestisto5,783yearsago(95%HPDinterval:5,021–
25),andY.pestis(n=139),totheY.pseudotuberculosisrefer- 7,022 years ago). This is also significantly older and with a
ence genome (IP32953). Only high confidence positions were much narrower confidence interval than previous findings of
extracted. To assess whether the individuals were infected 3,337yearsago(1,505–6,409yearsago)(Cuietal.,2013).
with multiple strains of Y. pestis we investigated the genotype
heterozygosity levels of the ancient genomes and found no BronzeAgeYersiniapestisStrainsLackingYersinia
indicationsofmixedinfection(FigureS3).Therewasnodecay MurineToxin
in Linkage Disequilibrium (LD) across the chromosome (Fig- Forthehigh-depthancientY.pestisgenomes,weinvestigated
ureS3),indicatingnodetectablerecombinationamongstrains. the presence of 55 genes that have been associated with the
We therefore used RAxML (Stamatakis, 2014) to construct a virulenceofY.pestis(Figure5A,TableS6).Wefoundallvirulence
Maximum Likelihood phylogeny from a supermatrix concate- genestobepresent,excepttheYersiniamurinetoxin(ymt)gene
natedfrom3,141genesandatotalof3.14Mbp(Figure4).This thatislocatedat74.4–76.2kbonthepMT1plasmid(Figure2C,
contrasts with earlier phylogenies (Bos et al., 2011; Cui et al., arrow1).TheymtgeneencodesaphospholipaseDthatprotects
Cell163,571–582,October22,2015ª2015TheAuthors 573
A B
D
C
Figure2. Y.pestisDepthofCoveragePlots
(A–D)Depthofcoverageplotsfor(A)CO92chromosome,(B)pCD1,(C)pMT1,(D)pPCP1.Outerring:Mappability(gray),genes(RNA:black,transposon:purple,
positivestrand:blue,negativestrand:red),RISE505(blue),RISE509(blue),Justinianplague(orange),BlackDeathplague(purple),modernY.pestisD1982001
(green),Y.pseudotuberculosisIP32881(red)sample.ThemodernY.pestisandY.pseudotuberculosissamplesareincludedforreference.Thehistogramsshow
sequencedepthin1kbwindowsforthechromosomeand100bpwindowsfortheplasmidswithamaxof20Xdepthforeachring.Arrow1:ymtgene,arrow2:
transposonatstartofmissingregiononpMT1,arrow3:transposonatendofmissingregiononpMT1,arrow4:plagene,arrow5:missingflagellinregionon
chromosome.TheplotsweregeneratedusingCircos(Krzywinskietal.,2009).SeealsoTablesS2,S3andS8.
Y.pestisinsidethefleagut,thusenablingthisentericbacteria (modern and ancient), except three strains (5761, 945, and
touseanarthropodasvector;itfurtherallowsforhighertiters CA88)thatarelackingthepMT1plasmidcompletely.
of Y. pestis and higher transmission rates (Hinnebusch, 2005; Althoughlargersamplesizesareneededforconfirmation,our
Hinnebuschetal.,2002).Wheninvestigatingallsevensamples data indicate that the ymt gene was not present in Y. pestis
forthepresenceofymt,weidentifieda19kbregion(59–78kb, before1686calBC(n=6),whileafter951calBC,itisfoundin
Figure2Carrow2–3,Figure5B)tobemissingexceptintheyoun- 97.8%ofthestrains(n=140),suggestingalateandveryrapid
gestsample(RISE397,951calBC)(Figure5B,TableS7).Wefind spread of ymt. This contrasts with previous studies arguing
thisregiontobepresentinallotherpublishedY.pestisstrains thattheymtgenewasacquiredearlyinY.pestisevolutiondue
574 Cell163,571–582,October22,2015ª2015TheAuthors
to its importance in its life cycle (Carniel, 2003; Hinnebusch, DifferentRegion4PresentintheAncestralYersinia
2005; Hinnebusch et al., 2002; Sun et al., 2014). Interestingly, pestis
we identified two transposase elements flanking the missing Besides the 55 pathogenicity genes, we also investigated the
19kbregion,confirmingthattheymtgenewasacquiredthrough presenceofdifferentregion4(DFR4)thatcontainsseveralgenes
horizontalgenetransfer,aspreviouslysuggested(Lindleretal., withpotentialroleinY.pestisvirulence(Radnedgeetal.,2002).
1998).Moreover,ithasrecentlybeenshownthatthetransmis- ThisregionwasreportedaspresentinthePlagueofJustinian
sionofY.pestisbyfleasisalsodependentonlossoffunction andBlackDeathstrains,havingbeenlostintheCO92reference
mutationsinthepde2,pde3,andrcsAgenes(Sunetal.,2014). genome(fromtheThirdPandemic)(Chainetal.,2004;Wagner
The RISE509 sample carries the promoter mutation of pde3 etal.,2014).Consistentwiththeancestralpositionofoursam-
andthefunctionalpde2andrcsAalleles(FigureS6).Incombina- ples,wefindevidencethattheregionispresentinallofourseven
tionwiththeabsenceofymt,theseresultsstronglysuggestthat samples(FigureS6).
theancestralY.pestisbacteriaintheseearlyBronzeAgeindivid-
ualswerenottransmittedbyfleas. YersiniapestisflagellarFrameshiftMutationAbsentin
BronzeAgeStrains
NativePlasminogenActivatorGenePresentinBronze AnotherimportantfeatureofY.pestisistheabilitytoevadethe
AgeYersiniapestis mammalianimmunesystem.Flagellinisapotentinitiatorofthe
Another hallmark gene of Y. pestis pathogenicity is the plas- mammalian innate immune system (Hayashi et al., 2001).
minogen activator gene pla (omptin protein family), located on Y. pseudotuberculosis is known to downregulate expression
thepPCP1plasmid(6.6–7.6kb).Thegenefacilitatesdeeptissue of flagellar systems in a temperature-dependent manner, and
invasionandisessentialfordevelopmentofbothbubonicand none of the known Y. pestis strains express flagellin due to a
pneumonic plague (Sebbane et al., 2006; Sodeinde et al., frameshift mutation in the flhD regulatory gene (Minnich and
1992; Zimbler et al., 2015). We identify the gene in six of the Rohde, 2007). However, we do not find this mutation in either
sevengenomes,butnotinRISE139,thesamplewiththelowest RISE505orRISE509,suggestingthattheyhavefullyfunctional
overalldepthofcoverage(0.75XonpPCP1)(Figure2D,arrow4, flhDgenesandthatthelossoffunctionoccurredafter2746cal
Table S6). Recently, it has been proposed that pPCP1 was BC.Interestingly,theyoungestofthesetwoY.pestisgenomes
acquiredafterthebranchingofthe0.PE2clade(Zimbleretal., (RISE505, 1686 cal BC) shows partial loss of one of the two
2015);however,weidentifiedpPCP1inoursamples,including flagellasystems(758–806kb),with39of49genesdeleted(Fig-
in the 0.PE7 clade (strains 620024 and CMCC05009), which ure2A,arrow5,TableS8).Thisdeletionwasnotfoundinanyof
divergedpriortothecommonancestorofthe0.PE2lineage(Fig- theotherY.pestissamples(n=147).Thismaypointtoselective
ure4B,Figure5A).ThisshowsthatpPCP1andplalikelywere pressureonancestralY.pestiswhenemergingasamammalian
present in the most basal Y. pestis (RISE509), suggesting that pathogen,yieldingvariablyadaptivestrains.
the 0.PE2 strains lost the pPCP1 plasmid. Interestingly, three
2.ANT3 strains (5761, CMCC64001, and 735)are also missing DISCUSSION
theplagene,indicatingthatthelossofpPCP1occurredmore
thanonceintheevolutionaryhistoryofY.pestis. Ourcalibratedmolecularclockpushesthedivergencedatesfor
Additionally,weinvestigatedwhetherRISE397,RISE505,and theearlybranchingofY.pestisbackto5,783yearsago,anaddi-
RISE509hadtheisoleucinetothreoninemutationataminoacid tional2,000yearscomparedtopreviousfindings(TableS5,Fig-
259 in the Pla protein. This mutation has been shown to be ureS5)(Cuietal.,2013;Morellietal.,2010).Furthermore,using
essential for developing bubonic, but not pneumonic, plague thetemporallystampedancientDNAdata,weareabletoderive
(Zimbleretal.,2015).Wefoundthatthesesamples,inagreement atimeseriesforthemolecularacquisitionofthepathogenicity
withtheirbasalphylogeneticposition,carrytheancestralisoleu- elements and immune avoidance systems that facilitated the
cineresidue.However,wealsoidentifiedavalinetoisoleucine evolution from a less virulent bacteria with zoonotic potential,
mutationatresidue31forRISE505(1686calBC)andRISE509 suchasY.pseudotuberculosis,tooneofthemostdeadlybacte-
(2746calBC).Thismutationwasnotfoundinanyoftheother riaeverencounteredbyhumans(Figure6).
140Y.pestisstrains,butwaspresentinotheromptinproteins, From our findings, we conclude that the ancestor of extant
suchasEscherichiacoliandCitrobacterkoseri,andverylikely Y. pestis strains was present by the end of the 4th millennium
represents the ancestral Y. pestis state. The youngest of the BCandwaswidelyspreadacrossEurasiafromatleasttheearly
samples, RISE397 (951 cal BC) carries the derived isoleucine 3rd millennium BC. The occurrence of plague in the Bronze
residue,showingthatthismutation,similartotheacquisitionof Age Eurasian individuals we sampled (7 of 101) indicates that
ymt,wasonlyobservedafter1686calBC. plague infections were common at least 3,000 years earlier
Analternativeexplanationtotheacquisitionofymtandthepla than recorded historically. However, based on the absence of
I259T mutation, given the disparate geographical locations of crucial virulence genes, unlike the later Y. pestis strains that
our samples, could be that the Armenian strain (RISE397, 951 wereresponsibleforthefirsttothirdpandemics,theseancient
calBC)containingymtandtheisoleucineresidueinplahada ancestralY.pestisstrainslikelydidnothavetheabilitytocause
longerhistoryintheMiddleEastandexperiencedanexpansion bubonicplague,onlypneumonicandsepticemicplague.These
duringthe1stmillenniumBC.Thiswouldhaveledtoitsexportto early plagues may have been responsible for the suggested
Eurasiaandpresumablytheextinctionoftheothermoreances- populationdeclinesinthelate4thmillenniumBCandtheearly
tralandlessvirulentY.pestisstrains. 3rdmillenniumBC(Hinzetal.,2012;Shennanetal.,2013).
Cell163,571–582,October22,2015ª2015TheAuthors 575
A B
C
D
E F
(legendonnextpage)
576 Cell163,571–582,October22,2015ª2015TheAuthors
It has recently been demonstrated by ancient genomics as well as bases with quality 2 or less. Hereafter, the trimmed reads
that the Bronze Age in Europe and Asia was characterized with a length of at least 30 nt were mapped using bwa mem (local
by large-scale population movements, admixture, and re- alignment) (Li and Durbin, 2009) to the database of Y. pestis,
Y. pseudotuberculosis, Y. similis, and Y. enterocolitica mentioned above.
placements (Allentoft et al., 2015; Haak et al., 2015), which
Readswithamatchtoanyofthesequencesinthisdatabasewerealigned
accompanied profound and archaeologically well-described
separately to three different reference genomes: Yersinia pestis CO92
social and economic changes (Anthony, 2007; Kristiansen genome includingtheassociatedplasmidspCD1,pMT1,pPCP1(Parkhill
and Larsson, 2005). In light of our findings, it is plausible et al., 2001); Yersinia pseudotuberculosis IP32953 including the associ-
that plague outbreaks could have facilitated—or have been ated plasmids (Chain et al., 2004); Yersinia pestis biovar Microtus 91001
facilitated by—these highly dynamic demographic events. andassociatedplasmids(Zhouetal.,2004).Thisalignmentwasperformed
using bwa aln (Li and Durbin, 2009) with the seed option disabled for
However, our data suggest that Y. pestis did not fully adapt
better sensitivity for ancient data, enforcing global alignment of the
as a flea-borne mammalian pathogen until the beginning of
readtothereferencegenome.Eachsequencingrunwasmergedtolibrary
the 1st millennium BC, which precipitated the historically re- level and duplicates removed using Picard-1.124 (http://broadinstitute.
corded plagues. github.io/picard/), followed by merging to per sample alignment files.
These files were filtered for a mapping quality of 30 to only retain high
quality alignments and the base qualities were re-scaled for DNA
EXPERIMENTALPROCEDURES
damage using MapDamage 2.0 (Jo´nsson et al., 2013). We defined
Y.pestisaspresentinasampleifthemappeddepthoftheCO92refer-
SamplesandArchaeologicalSites
ence sequences were higher or equal to 0.1X and if the reads covered
Weinitiallyre-analyzedthedatafromAllentoftetal.(Allentoftetal.,2015)and
atleast10%ofthechromosomeandeachoftheplasmids.Theassembly
identifiedY.pestisDNAsequencesin7ofthe101individuals.Descriptionsof
of Justinian, Black Death, and the modern samples were performed
thearchaeologicalsitesaregiveninSupplementalExperimentalProcedures
similarly and is described in detail in the Supplemental Experimental
andTableS1.
Procedures.
GenerationofAdditionalSequenceData
Coverage,DepthandMappabilityAnalyses
In order to increase the depth of coverage on the Y. pestis genomes we
We calculated the coverage of the individual sample alignments versus
sequencedmoreonthesesevenDNAextracts.Libraryconstructionwascon-
the Y. pestis CO92 reference genome using Bedtools (Quinlan and Hall,
ductedasin(Allentoftetal.,2015).Briefly,doublestrandedandblunt-ended
2010) and plotted this using Circos (Krzywinski et al., 2009). For the
DNAlibrarieswerepreparedusingtheNEBNextDNASamplePrepMaster
chromosome,thecoveragewascalculatedin1kbpwindowsandforthe
MixSet2(E6070)andIllumina-specificadapters(MeyerandKircher,2010).
plasmids in 100 bp windows. Mappability was calculated using GEM-
Thelibrarieswere‘‘shot-gun’’sequencedintwopoolsonIlluminaHiSeq2500
mappabilitylibraryusingak-mersizeof50,whichissimilartotheaverage
platformsusing100-bpsingle-readchemistry.Wesequenced32lanesgener-
length of the trimmed and mapped Y. pestis reads (average length
atingatotalof11.2billionnewDNAsequencesforthisstudy.Readsforthe
43–65 bp). Statistics of the coverage and depth are given in Tables S3
sevenY.pestissamplesareavailablefromENA:PRJEB10885.Individualsam-
andS4.
pleaccessionsnumbersareavailableinTableS2.
DNADecayRates
CreationofDatabaseforIdentificationofY.pestisReads We investigated the molecular degradation signals obtained from the
ToidentifyY.pestisreadsintheBronzeAgedataset(Allentoftetal.,2015)we sequencing data. Based onthe negative exponential relationship between
firstcreatedadatabaseofallpreviouslysequencedY.pestisstrains(n=140), frequency and sequence length, we estimated for each sample the DNA
Y.pseudotuberculosisstrains(n=30),Y.similisstrains(n=5),andaselection damagefraction(l,perbond),theaveragefragmentlength(1/l),theDNA
ofY.enterocoliticastrains(n=4)(SupplementalExperimentalProceduresand decayrate(k,perbondperyear),andthemolecularhalf-livesof100bpfrag-
TableS2).ThegenomeswereeitherdownloadedfromNCBIordownloadedas ments(Allentoftetal.,2012).WecomparedtheseDNAdecayestimatesfor
readsanddenovoassembledusingSPAdes-3.5.0(Bankevichetal.,2012) Y.pestistothedecayofendogenoushumanDNAfromthehostindividuals.
withthe–carefuland–cov-cutoffautooptions. IftheplagueDNAisauthenticandancient,acorrelationisexpectedbetween
the rate of DNA decay in the human host and in Y. pestis, because the
IdentificationandAssemblyofY.pestisFromAncientSamples DNA has been exposed to similar environmental conditions for the same
Raw reads were trimmed for adaptor sequences using AdapterRemoval- amountoftime.SeeSupplementalExperimentalProceduresforadditional
1.5.4(Lindgreen,2012).AdditionallyleadingandtrailingNswereremoved information.
Figure3. AuthenticityofY.pestisDNA
(A)DNAdamagepatternsforRISE505andRISE509.ThefrequenciesofallpossiblemismatchesobservedbetweentheY.pestisCO92chromosomeandthe
readsarereportedingrayasafunctionofdistancefrom50(leftpanel,first25nucleotidessequenced)anddistanceto30(rightpanel,last25nucleotides).The
typicalDNAdamagemutationsC>T(50)andG>A(30)arereportedinredandblue,respectively.
(B)AncientDNAdamagepatterns(n=7)ofthereadsalignedtotheCO92chromosomeandtheY.pestisassociatedplasmidspMT1,pCD1andpPCP1.The
boxplotsshowthedistributionofC-Tdamageinthe50ofthereads.Thelowerandupperhingesoftheboxescorrespondtothe25thand75thpercentiles,the
whiskersrepresentthe1.5inter-quartilerange(IQR)extendingfromthehinges,andthedotsrepresentoutliersfromthese.
(C)DNAfragmentlengthdistributionsfromRISE505andRISE509samplesrepresentingboththeY.pestisDNAandtheDNAofthehumanhost.Thedecliningpart
ofthedistributionsisfittedtoanexponentialmodel(red).
(D)Linearcorrelation(red)betweenthedecayconstantintheDNAofthehumanhostandtheassociatedY.pestisDNAextractedfromthesameindividual
(R2=0.55,p=0.055).Thedecayconstant(l)describesthedamagefraction(i.e.,thefractionofbrokenbondsontheDNAstrand).
(E)DistributionofeditdistanceofhighqualityreadsfromRISE505andRISE509samplesmappedtoeitherY.pestis(darkgray)orY.pseudotuberculosis(light
gray)referencegenomes.ThereadshaveahigheraffinitytoY.pestisthantoY.pseudotuberculosis.
(F)Plotsofactualcoverageversusexpectedcoverageforthe101screenedsamples.Expectedcoveragewascomputedtakingintoaccountreadlengthdis-
tributions,mappablefractionsofreferencesequences,andthedeletionsinpMT1forsomeofthesamples.SamplesassumedtocontainY.pestisareshownin
blueandRISE392thatisclassifiedasnotY.pestisappearsisshowninred.SeealsoFigureS1andS2,TableS3.
Cell163,571–582,October22,2015ª2015TheAuthors 577
A Figure4. PhylogeneticReconstructions
BC***60RR0IISS002EE0..0550000E998-O40((Knn5==581165*))000*0 *** 4000.0P0E7*(*n*=2*)3*0*0**0O*0K*6B00*8l.a8AR*c04*NkI..*3SAA*T.DA0ENN23*P*e.2N15 T1.1TTA*(aA*6100.1.n.T3 8N12AtII(12.N512.0h1= N N*nI20..I2J(TN.(.N.TI.A2N0=nIA.n.(1u1I(2.NI1*N2PATn)N nA1N0=N1=1s.N( .2=IE2(N=(12n)(5t4N2M22NTn(nTin(1.1=)4n)nT(.AnT02(=(==E32)()(Mn5inn=1na.=(n21N81D)=P(n==(1I2(En=1=))P)n(nTn12=E(23.)31Dn)1=nMp==3 )39))313=()1131=l1)2*n2Ea7)13(. ()3370)=(ngOnD10)n))5INu=Y=PR01=)79e311 10113IY(IY)I282 P0P6(n077(5)3n32)1n=4222=I69*P54=3143)1843)285II*PP)0332360730I08HIP.12P031821E2189564852II5I03PP1P4/333(0222n2894=8261132N)o5 *RC6IPPS194BA2CCA4C4RMC77C9MMCMM1CC4CMCCUACCCCCK7D37ZCD3CC32CND5G9CCCCHMDD3CCCCGCHACCCFCECIYEYYC2K5DCJ97AC7CJIIICCCCIISSCIKCACAMIMCE222FYFFCYYCCYe1eE121211n3P2o22n14IGM21M61M181M5407n411MGMM0MGM7M13GGGMG23MMMMGMMNM1N1MNM1M1eHH1MMMMMM1M1M1MM1I11M1M11M15SM56VMMMMn11Ms111OA11iMNNNNsMMMG6C5990sd990M3t002c099950t4001000972901639t9p95821999999CJ909igJJJJC9Et0058C99701CAAJ99411CCCiCCCC_0C8CCCCCCCCCCCCCC0CCCCCoCCCi6C70969C9CCC22q2r1C7o0CC896q720CC5C94ZZZ7Z30Z607150a297a7Z17789806685o7956659492o57694984206mNNC8CiCC02C6CCiC9C015u551C41CCCCCC3CCCCCCCu3CCCCCCCC8C5CC21CCCCCd1639C765106dC40410581l77198CC2C485097l428319604t66206148253011a587a4000003a2070ee8u244N1o1e90965001130235990991188308065909107111996456100310000000000053K8K003121_9181980_00100rs100223ss811314023015741025664400444000229270305000010000007000g0104108116_0B0_U103_1251612000000000000000070000004608001010051111170431e0061411104461201A0F000190314801940019010G21000034300002000000002d020031611729201404112432000117860030173341000010000077132222415000051113004 STEWgtpSsfY(Ysag(pcb1uusl(mgeY(pYs0YenBCArhieuarls9nlaiiilhrh.....oe5xfiansavvveoeae)ae))etmpp9cbyypemppetdeeetgpi.ptdtnaBmrniiopll9slersnnnoenooetpaEeeseiwcnMMseEh,aloesSdenaYgglsisclsxhdhianaeuSsnitatAura2naauttsttieeu.ptoieotpsscsdiihissd30xTeS.xamdosannpeswpasosa(F.oi0siardrnsoTyytdpmcbmyeetcapeeNtr(hna(ui4tbeest2megtuirlelidsultevevuae)nnuebau;eeodlega.butebocecmembcadmdismdisdiCd(irnatsreTnFfweorvAikgserstrSe)fie,gaceriuc.edshiersl)iret)giccltEf5scxtnetuYceioaebYneaihuysnuTShcuLuL.uixnttrupv.erheme.satashiariloi4aeplafwcaoekoakltsfooepwuelloetvs.ouneeAner)neeessepwsiclaqssaatShrmldrillshxdEsidMchessoaieliitiusrgihae.shmt5rpunoxse,tusvrdepweeTnatooc.suwbgeetw,pbutdehem(n2aeehajdelootbSiecsoawdeaayretoinorcn0trebtrdadholttaeerrsdlorhittuiin.1aiesilisnsgutairstnepemCet>pFenl3suidehgkbgyuetl.bcrao)osrr)ertaoereie9crb,Sc.eealenivetPpoaoruc.olahnvr5dt(aessernB2ccStoueOrhnafdoncdeteh%eomalreondooacrltdeacavcisua:seffprseanncYnn1uluoenhiFoeubialttbdnttbgpnegos.s5illeirrhislaroaglnwPgicodFaopteetsa(lrttn9(rrtistndgseuhuisigrnirorAutenueuehassgc5eoyeeesiroerstscbrcd,kccsrnuei%eescerainteeettimnhthtadoaoltritrescioiaeeonsaSieergootesDivnr()ledclamsdntfneo.b.sca3nencacaiboupYSgendvrYllpl,orYYalOuadiosrrou)tse.4ggs.eSrend..heddhoondobbnrrctisurs4isneg)eshofippwtlnffiarreitpoynssm,liaaaep.siweegamntoipratgaaatTnnntossrtiaatttiihanmolnnneibannhcchhhcattpialtrrliasiiirddddpgog.nhhneeeeeeeesssl--t,.ll
IP32953 (cid:2) (cid:3)
260 c=1(cid:3)YN 1(cid:3)li ri
IH111554 g
IP32463 i=1
IP32921
wherethereadshaveNdifferentlengths,l tol with
IP32881 1 N
countsr tor .Toaccountformappabilitywedeter-
RISE509 1 N
RISE505 mined the mappable fraction for each reference
0.PE7 sequence using kmers of length 40, 50, and 60,
0.PE7 andthenusedthemappabilityvaluewiththek-mer
0.PE2 lengthclosesttotheactualaveragereadlengthfor
7000.0 0.PE2 eachsample/referencecombination.Formoreinfor-
0.PE3
mationseeSupplementalExperimentalProcedures.
0.PE4A
0.PE4A
GenotypingForPhylogeneticAnalyses
Branch 0-4 AIcPoli3gn2ns9me5ne3snuwtssaossfeauqlsluesedtnracaiensssrvefoeferrsreaunlslcYse.apmfosprelugedseontuoustbyeepdricniunglottshhiees
phylogeny.Thesamplesweregenotypedindividually
usingsamtools-0.1.18andbcftools-0.1.17(Lietal.,
ComparisonofSamplestoY.pestisandY.pseudotuberculosis 2009)andhereafterfiltered(SupplementalExperimentalProcedures).Based
ReferenceGenomes on Y. pseudotuberculosis IP32953 gene annotations, the consensus se-
We used the alignments of several sets of reads (Y. pestis, quencesforeachgeneandsamplewereextracted.Becauseofthedivergence
Y. pseudotuberculosis, and Y. similis) to Y. pestis CO92 and the betweenY.pestisandY.pseudotuberculosis,anumberofgenesequences
Y.pseudotuberculosisIP32953genomes.Persamplewedeterminedthedis- displayedhighratesofmissingbasesandweremovedgeneswhere20or
tribution of edit-distances (mismatches) of the reads versus the particular moremodernY.pestissampleshad>10%missingness.Thiscorresponded
referencegenome.WeusedthesedistributionstobuildaNaiveBayesianclas- toatotalof985genes,leavingdatafrom3,141genesthatweremergedinto
578 Cell163,571–582,October22,2015ª2015TheAuthors
A Figure5. IdentificationofVirulenceGenes
(A)Genecoverageheatmapof55virulencegenes(rows)
in 140 Y. pestis strains (columns). Sample ordering is
basedonhierarchicalclustering(notshown)ofthegene
coverage distributions. RISE505 and RISE509 are
markedwitharedasterisk.Coloringgoesfrom0%gene
coverage(white)to100%genecoverage(blue).
(B) Depth of coverage of high quality reads mapping
across pMT1. Outer ring is mappability (gray), genes
(RNA:black,transposon:purple,positivestrand:blue,
negativestrand:red)andthentheRISEsamplesordered
afterdirectAMSdating.SampleorderingareRISE509,
RISE511, RISE00, RISE386, RISE139, RISE505 and
RISE397. See also Figure S6, Tables S2, S6, and S7.
AMS:AcceleratorMassSpectrometry.
B SNVs,theLDr2wascalculatedusingPLINK1.9(Chang
etal.,2015)andplottedagainstthephysicaldistancebe-
tweenthepairs.Wereconstructedthephylogenyfrom
thecodon-partitionedsupermatrixusingRAxML-8.1.15
(Stamatakis,2014)withtheGTR+G+Isubstitutionmodel.
Bootstrapswereperformedbygenerating100bootstrap
replicates and their corresponding parsimony starting
trees using RAxML. Hereafter, a standard Maximum
Likelihoodinferencewasrunoneachbootstrapreplicate,
andtheresultingbesttreesweremergedanddrawnon
thebestMLtree.InitialphylogeniesplacedtheY.pestis
Harbin strain with an unusual long branch inside the
1.ORIcladeanditwasexcludedfromfurtheranalysis.
Additionally Y. pseudotuberculosis SP93422 (serotype
O:15), Y. pseudotuberculosis WP-931201 (serotype
O:15) and Y. pseudotuberculosis Y248 (serotype un-
known) was in a clade with long branch lengths and
werethereforealsoomitted(seeFigureS4).
HeterozygosityEstimates
We determined heterozygosity by down-sampling the
Y.pestisbam-filestothesameaveragedepthasthecorre-
spondingRISEsamples,genotypedeachofthesamples
andextractedheterozygotecallswithadepthequaltoor
higherthan10.Alltransitionswereexcluded.SeeSupple-
mentalExperimentalProceduresfordetailedinformation.
DivergenceEstimations
TodatethedivergencetimeforY.pestisandnodeswithin
theY.pestiscladeweperformedBayesianMarkovChain
Monte Carlo simulations using BEAST-2.3.0 (Bouck-
aertetal.,2014)andtheBEAGLElibraryv2.1.2(Ayres
asupermatrix.Wecreatedtwodifferentsupermatrices,onewithY.similis, et al., 2012). We used the codon-partitioned supermatrix that included the
Y.pseudotuberculosis,andY.pestiscontaining173taxa33,141genesthat twoclosestY.pseudotuberculosisclades,withunlinkedsubstitutionmodels,
was used for the initial phylogeny (Figure 4A). The second supermatrix GTR+G+Iwitheightgammaratecategoriesandunlinkedclockmodels.Dates
consisted of all Y. pestis strains and the genomes from the two closest weresetasyearsagowiththeRISE509,RISE505,JustinianandBlackDeath
Y. pseudotuberculosis clades, which was used for the divergence time samplessetto4,761,3,701,1,474,and667yearsago,respectively.Allunknown
estimations. datesweresetto0yearsago.Wefollowedpreviouswork(Cuietal.,2013;Wag-
neretal.,2014)andappliedalognormalrelaxedclock,assumingaconstantpop-
Phylogenetics ulationsize.Were-rootedtheMLtreefromRAxMLsothattherootwasplaced
Thealignmentswerepartitionedbycodonpositionandanalyzedwithjmodelt- betweenthetwoY.pseudotuberculosisclades(IP32953,260,IH111554)and
est-2.1.7(Darribaetal.,2012)totestforthebestfittingsubstitutionmodel.All (IP32921, IP32881, IP32463) and used this as the starting tree. Based on
decisioncriteria(Akaike,Bayesian,andDecisiontheory)foundtheGeneralized the ML tree we defined the closets Y. pseudotuberculosis clade (IP32921,
Time Reversible substitution model with gamma distributed rates, using IP32881,IP32463)andtheY.pestiscladeasamonophyleticgroupanddefined
fourratecategories,andaproportionofinvariablesites(GTR+G+I)tobethe auniformpriorwith1,000and100,000yearsasminimumandmaximumbounds.
bestfitforeachofthethreecodonpartitions.Totestforrecombinationacross Weran20independentparallelBEASTchainssamplingevery2,000statesfor
thechromosomeweestimatedlinkagedisequilibrium(LD)using141Y.pestis between52and64millionstatesusingatotalof240,000corehours.Thechains
strains.Atotalof482bi-allelicsinglenucleotidevariations(SNVs),withaminor werecombinedusingLogCombinerdiscardingtheinitial10millionstatesas
allelefrequencyof5%orhigherwereextracted.Forallpairsoftheextracted burn-in.Thecombinedpostburn-indatarepresented961millionstatesand
Cell163,571–582,October22,2015ª2015TheAuthors 579