Using ancient DNA to study the origins and dispersal of ancestral Polynesian chickens across the Pacific VickiA.Thomsona,Ophélie Lebrasseurb,JeremyJ.Austina,c,1,TerryL.Huntd,e,DavidA.Burneyf,TimDenhamg, NicolasJ.Rawlencea,h,JamieR.Woodi,JaimeGongoraj,LinusGirdlandFlinkb,k,AnnaLinderholmb,KeithDobneyl, GregerLarsonb,andAlanCoopera,1 aAustralianCentreforAncientDNA,SchoolofEarthandEnvironmentalSciences,UniversityofAdelaide,Adelaide,SA5005,Australia;bDurhamEvolutionand AncientDNA,DepartmentofArchaeology,andkSchoolofBiologicalandBiomedicalSciences,DurhamUniversity,DurhamDH13LE,UnitedKingdom; cSciencesDepartment,MuseumVictoria,Melbourne,VIC3001,Australia;dClarkHonorsCollegeandeDepartmentofAnthropology,UniversityofOregon, Eugene,OR97403;fNationalTropicalBotanicalGarden,Kalaheo,HI96741;gSchoolofArchaeologyandAnthropology,ANUCollegeofArtsandSocial Sciences,TheAustralianNationalUniversity,Canberra,ACT0200,Australia;hAllanWilsonCentreforMolecularEcologyandEvolution,Departmentof Zoology,UniversityofOtago,Dunedin9016,NewZealand;iLandcareResearch,Lincoln7640,NewZealand;jFacultyofVeterinaryScience,Universityof Sydney,Sydney,NSW2006,Australia;andlDepartmentofArchaeology,UniversityofAberdeen,AberdeenAB243UF,Scotland EditedbyDavidJ.Meltzer,SouthernMethodistUniversity,Dallas,TX,andapprovedFebruary20,2014(receivedforreviewOctober31,2013) The human colonization of Remote Oceania remains one of the resolved, the precise details of this intensive migratory episode great feats of exploration in history, proceeding east from Asia remainunclear(6). acrossthe vast expanseof the PacificOcean.Human commensal Humancommensalandearlydomesticatedspecieswerewide- anddomesticatedspecieswerewidelytransportedaspartofthis ly,butnotubiquitously,dispersedaspeoplecolonizedthePacific. diaspora, possibly as far as South America. We sequenced mito- As a result, they provide an opportunity to study colonization chondrial control region DNA from 122 modern and 22 ancient eventsandsubsequentmovementsforislandsandregionswhere they were successfully introduced, especially through the use of chicken specimens from Polynesia and Island Southeast Asia and biomolecular techniques, including ancient DNA. In the Asia– usedthesetogetherwithBayesianmodelingmethodstoexamine Pacific region, the complex histories of Pacific island colo- the human dispersal of chickens across this area. We show that nizationshavebeeninvestigatedusingthebiologicalelements specifictechniquesareessentialtoremovecontaminatingmodern associated with these cultures, such as bottle gourds (7, 8), DNA from experiments, which appear to have impacted previous sweet potatoes (9), pigs (10, 11), dogs (12), Pacific rats (13), studiesofPacificchickens.Incontrasttopreviousreports,wefind and chickens (14–17). However, studies of commensals and that all ancient specimens and a high proportion of the modern domesticates in the Pacific to date have provided limited res- chickenspossessagroupofunique,closelyrelatedhaplotypes olution of dispersal routes, due to low amounts of genetic di- found only in the Pacific. This group of haplotypes appears to versity in many groups and overwriting of genetic signals by represent the authentic founding mitochondrial DNA chicken subsequent introductions, especially for cotransported species lineages transported across the Pacific, and allows the early like rats (10, 13, 18). dispersal of chickens across Micronesia and Polynesia to be Ancient and modern DNA from chickens provide an oppor- modeled. Importantly, chickens carrying this genetic signature tunitytostudyhuman-mediateddispersalacrossthePacificdue persistonseveralPacificislandsathighfrequencies,suggesting totheextentofgeneticandphenotypicdiversityandtherangeof that the original Polynesian chicken lineages may still survive. archaeological material available. Although recent studies of No early South American chicken samples have been detected with the diagnostic Polynesian mtDNA haplotypes, arguing Significance against reports that chickens provide evidence of Polynesian contact with pre-European South America. Two modern speci- AncientDNAsequencesfromchickensprovideanopportunity mensfromthePhilippinescarryhaplotypessimilartotheancient tostudytheirhuman-mediateddispersalacrossthePacificdue Pacificsamples,providingcluesaboutapotentialhomelandfor tothesignificantgeneticdiversityandrangeofarchaeological thePolynesianchicken. material available. We analyze ancient and modern material | | | | andrevealthatpreviousstudieshavebeenimpactedby con- Lapita Pacificcolonization phylogeography archaeology migration tamination with modern chicken DNA and, that as a result, thereisnoevidenceforPolynesiandispersalofchickenstopre- ThecolonizationoftheremotePacificwasoneofthelastgreat Columbian South America. We identify genetic markers of humanmigrations,butdespitetherecentnatureoftheevents, authenticancientPolynesianchickensandusethemtomodel thetimingandroutesremainanareaofconsiderabledebate.The earlychickendispersalsacrossthePacific.Wefindconnections first colonization of Western Polynesia occurred around 3,250– betweenchickensintheMicronesianandBismarckIslands,but 3,100 calendar years before present (cal B.P.) as part of the noevidencethesewereinvolvedindispersalsfurthereast.We eastward migration of Lapita pottery-bearing peoples (1). This also find clues about the origins of Polynesian chickens in migrationoccurredonlyafewhundredyearsaftertheemergence thePhilippines. of this distinctive pottery tradition in the Bismarck Archipelago around 3,470–3,250 cal B.P., although its antecedents can be Authorcontributions:V.A.T.,J.J.A.,J.G.,andA.C.designedresearch;V.A.T.,O.L.,N.J.R., tracedtoIslandSoutheastAsia(ISEA)(2–5).Followingtheinitial J.R.W.,L.G.F.,andA.L.performedresearch;T.L.H.,D.A.B.,andK.D.contributednew reagents/analytictools;V.A.T.,J.J.A.,andA.C.analyzeddata;andV.A.T.,J.J.A.,T.D.,G.L., movementintoWesternPolynesia,aprolonged1,800-yhiatus,or andA.C.wrotethepaper. “pause,” is apparent before further colonization (6), potentially Theauthorsdeclarenoconflictofinterest. relating to the need to develop sailing technology essential for ThisarticleisaPNASDirectSubmission. crossing the vast ocean barrier to the east (between Samoa and Datadeposition:Thesequencesreportedinthispaperhavebeendepositedinthe the Society Islands, 2,400 km; Fig. 1). The huge navigational GenBankdatabase(accessionnos.KJ000585–KJ000642). achievementofcolonizingtheremoteEastPolynesiantriangle(an 1Towhomcorrespondencemaybeaddressed.E-mail:[email protected] oceanicregionroughlythesizeofNorthAmerica)thenoccurred [email protected]. rapidly (<300 y) (6). Although the overall chronology of the Thisarticlecontainssupportinginformationonlineatwww.pnas.org/lookup/suppl/doi:10. eastern Pacific island colonization has recently been further 1073/pnas.1320412111/-/DCSupplemental. 4826–4831 | PNAS | April1,2014 | vol.111 | no.13 www.pnas.org/cgi/doi/10.1073/pnas.1320412111 140o E 180o E 140o W 100o W 40o N Legend 40o N Haplogroup D Fig.1. Mapshowingsamplesandlocalitiesmentionedinthis Haplogroup E OcNeeaanr iaORceemaontiea Other haplogroups 7 17 sintuDdya.nScaamuspeleestfarol.m(2V6)anaureatuunadnedrliGnuedam, wpitrhevhioaupslloygpruobuplishfreed- 20o NPhilippines 4030500B0P-(9) n=5 Hawaii 20o N quenciesofthechickenspecimensindicatedbypiecharts(thick 23 Guam 2B outlinesindicateancientsamples).Colorsrefertohaplotype/ 2A 20o 0 S o MNaolur1tkhLuCCP 3a4p590u-3a550BGPPNua(18eip)2nwueaaS13olom3PW1oolenysn tIeessr.nia SaI4sn3ltaa2n8 CdsruVza nuatuNi4ue2 5FSroecniceht yP oIslTy6.nueasGmiaaomMtobsaierrq I5su.esa7s R1a3pa Nui 200 oo S homaawinorratahserhrpt)ooene.elwrowtr1eshgs,sahyriisnpaonrienfposuteoprdltpoorhlri:rietdcdeeDyasusdpeaitchsneertfaristproopomranwrlomooropseoguvofaarteioserrncmeceushhdpebictacnolebikeannstyeoctssrknel(teo.3sdtadge8nGw,is)itcr;Entiaatea3hhlydcB,iadrnaiaopnaiyrnslNttreosaroeSgodtwSahr(di4Corfesu9fuOce)rdapt,ercniiefoepiwanfnanrelheynbtresoseileaiursfnsenc;etceat2(,hdnss,amciaancabeskrodlnlhtieavoeaeancerrsdkls---,l 40o S PEoalysnteersnia 40o S fromNewGuineaintoMicronesia;4–7,spreadofchickensfrom WesternPolynesiainto,andwithin,EasternPolynesia.Dashed lineindicatesdemarcationbetweenNearandRemoteOceania, 0 1,250 2,500 5,000 Kilometers WesternPolynesiaisdefinedbyadashedcircle,andEastern 140o E 180o E 140o W 100o W Polynesiaisindicatedbyagrayshadedtriangle. domesticchickenbreedshavehighlightedhowthedomestication Results process and subsequent breeding have resulted in a 70% loss The 61 WMG dataset (25) contained 363 single-nucleotide of nuclear genetic diversity (19), substantial phylogeographic polymorphisms (SNPs), of which 154 were potentially phyloge- Y structureremainswithinthemitochondrial(mtDNA)sequences netically informative, with 62 (17%) located in the rapidly G O of chickens worldwide (20–23). Furthermore, an extensive ref- evolving CR (32). Bayesian and maximum likelihood inference OL erencedatasetof>3,000mtDNAcontrolregion(CR)sequences analyses of the WMG dataset supported the haplogroup frame- OP R and >60 whole mtDNA genomes is available from across the work defined by Liu et al. (24) and Miao et al. (25) and, im- TH natural range of wild and semiwild birds, as well as domestic portantly, produced robust support forhaplogroups A–G and Z AN breeds of chicken, permitting the reconstruction of phylogeo- (i.e.,haplogroupswheremultipleindividualsweresequenced),as graphic patterns of domestic chickens and associated human shown in Fig. 2. Robust support values were also obtained for cultures. Despite these intensive surveys, a resolved worldwide phylogenetic trees based on the WMG data without the CR sequences(SIAppendix,Fig.S1),butwerelessrobustwhenbased chicken mtDNA phylogeny has not been developed, and this is only on the highly variable CR sequences alone, likely due to an essential prerequisite to interpreting short ancient DNA issueswithsubstitutionrateheterogeneity(32)(SIAppendix,Fig. sequences. The current phylogenetic framework for chickens is S2).However,theshort(201bp),hypervariableregionoftheCR based on Liu et al. (24), who identified nine highly divergent usedinpreviousstudiescontains>12×theaveragediversityper haplogroups(designatedA–I)usingmtDNACRsequences,with basecomparedwiththerestoftheWMG,andhastheadvantage an additional four recently described on the basis of whole of being available for a worldwide dataset of >1,000 chicken mtDNA genomes (W–Z) (25). However, there is little infor- sequences. The comparative phylogenetic dataset constructed mation about the support for these topologies, and only neigh- fromthesesequencesidentified274uniquehaplotypes,whichwe bor-joining treeshave beenreportedtodate. termedH001–H274(SIAppendix,Dataset S6). Phylogeographic studies have identified that one particular Ofthe37Polynesianarchaeologicalchickenbonesanalyzedto mtDNA lineage (CR haplogroup D) is largely limited to the study the temporal and spatial patterns within Polynesia, 22 Asia–Pacificregion(24),whereasmanyoftheotherhaplogroups (59%) yielded positive and repeatable PCR amplification and areubiquitousworldwide,potentiallyasaresultofhistoricaldis- DNAsequencingresultsfora330bpregion(whichincludedthe persalwithEuropeancolonialists(e.g.,haplogroupsA,B,andE), hypervariable201bp;Niue,n=2/8;Hawai’i,n=7/11;RapaNui, and are therefore generally phylogeographically uninformative. n=13/18;SIAppendix,TableS1).Allofthe22positiveancient Previousstudiesofmodernandancientchickenshaveidentified samples produced mtDNA CR sequences belonging to haplo- both haplogroup D and E in the Pacific (14–17, 26), making group D. Two samples that could not be reliably reproduced interpretation of colonization history difficult due to poten- (from Niue and Rapa Nui) each generated a single PCR prod- tially contrasting origins and dispersal histories (24). Indeed, uct with different non-D haplotypes (from haplogroup A and the presence of haplogroup E in the Pacific has been used E, respectively; SI Appendix, Table S1 and Dataset S1). How- to infer a link between Polynesia and pre-Columbian South ever, when DNase (double-strand–specific Shrimp DNase) pre- treatment was used to remove potential contaminating DNA America, although both the phylogenetic signal and radiocarbon dating of the samples have been questioned (27–29). This issue fromreagents(33),thesesequenceswerenolongerdetected(SI Appendix). Two of the 124 modern feather samples examined hasrecentlytakenonmoresignificanceasotherstudiesofancient couldnotbesuccessfullyamplified(onefromtheMarquesasand genetic diversity in South America emphasize the importance of another from Hawai’i). The large majority of the resulting 122 evidenceforpre-EuropeanPolynesiancontact(17,30). modernsequencesbelongedtohaplogroupD(n=90/122,74%; In this study, we first quantify the support for previously de- SIAppendix,Figs.S3–S7),withhaplogroupEsequencespresent finedchickenmtDNACRhaplogroupsusingrecentlypublished at a lower frequency (n = 27/122, 22%). The remaining five whole mitochondrial genomes (WMGs) (25). We then use the samples fell within haplogroups A, B, and I (n = 1, 3, and 1, resulting robust evolutionary framework to analyze the spatial respectively, each<2.5%). andtemporalpatternsofmtDNACRhaplotypesinancientand Previous studies of Pacific chickens have reported elevated modern Pacific chickens to examine their origins in ISEA (31), levels of haplogroup E among ancient specimens (up to 48%) the dispersal of chickens into Near Oceania and Western Poly- (14–16), in direct contrast to our results. However, the con- nesia,potentialconnectionsbetweentheNewGuinearegionand tamination of laboratory consumables with DNA from modern Micronesia,andtheclaimedintroductionofPolynesianchickens domesticspecies,includingchickens,isawell-knownproblemin toSouthAmerica(14). ancient DNA research (34), and this would also likely generate Thomsonetal. PNAS | April1,2014 | vol.111 | no.13 | 4827 Polynesian D founding C lineage (4 SNPs) Y GGWUU226611668811GG__GGXUU11GUU5522UGG__22666666CC2UU11111611666622GG77G66177UU0011755U11227770__66__GG111261111118888UU666_444333____22177665___1111__996116666CCCCCGG9CCCCCC__779____UU11_111111111CCCC2227766W66____1111_11G11CC776U99001112__0776CC__1_22633800X0_____CC1383_C2[[1100[0111]1000000][51]3NN4CC__NN00CC00__77002200773322[6633__177__2222211220____1DDDD111100GG]UU22066116688GG7777UU__22662222211333336688___22__DDDD22[449111__DD211][80111]GU000261683_25_D[200091060]1999888033[[0888GGGUUU000222666]]]11[669977__92277__9DD]33[91]GG1UU2266[1166880GG55__UU2288822__66DD033116667777]__2266__DD8[33587][110100[001000]]GU2616193_[012_10YGU0G2U6012]667146_[911601__0Z101_0Z0]NZCG_U0GG20U6U7212G2676U13112754067A5__10_B60G04098570_9U___0_8BB2B0669_6G1_B1B0U6229_6501G_610U8_324A_6_A1072B_0A0_04_A GU26A1I698_61_I G HGUG26U1GG27U1G6UGF021U2U_66625219146660_7117_G671857157_69__5__5G35D58_5Q6G_G_6_UHGG4G2U8G627U1627716G11677GU17_0U213_2566__51125461068_8_98_G__FF8U4F2F_96_147F702G_U_246F616_9F1_45_F 1[10010000[]100[399]1][196[93]0HQ80572019_4]1_E02[55]0GGUU226611770088__4422__EE33351000[[HHHH11HHHQQQQ0088880055QQQ[77]]]221188811__44554499__77EE332211122__4488833__EE033]]]00GGGUUUU2222AA6666PP1111777700AAAAAAA1111002222PPPPPP____333333330000004444GG11__0000EE88UU333311__3333332233111166__7777EE11____AA11AAA333377YY1111PPP__3322EE__AAAAA00033113355PPPP000055GG__00000333377EE000000GG55555500UUUU11GG333333__8888UU33333322222UU00001116622666622__999__AAA661166EE333__11111133YYY7766996677__882222HH__9900EE66EE333344QQ9911__11__5555__22338855553388995577755____77__EE11EE22EE11__111111330022____44EE0011__EE11 E F(bA2iolg5lo.).th2sa.BtrpaalyPopehg(syirinoalounpgpaepsrneownestttiheitchreitsoremerse)opsbrrueoapsbpteahodbartoinlivntayoWlnuaeeMnsdGianrddmeiavsathixadoiumfwraounlmmohnaMlvibkeiareaolrniohecbohtuoeassdlt.. 0.0010 supportandconcurwiththedesignationsofLiuetal.(24). haplogroup E sequences, due to the ubiquity of the latter our results suggest that if haplogroup E was present at all in worldwide (SI Appendix, Dataset S2) (24). To examine this po- ancientPacificchickens,itmusthavebeeninlessthan13%(at tentialexplanation,wereexaminedkeysamplesfromaprevious the95%probabilitylevel;SIAppendix,Fig.S9).Itispossiblethat studythatlinkedancientPolynesianchickenstoSouthAmerican if haplogroup E was present in very low frequencies among an- archaeological specimens (14). Four of the six bone samples cient Pacific chickens (e.g., <10%), we did not detect it within from Rapa Nui used in the previous study were available for the 22 ancient samples we examined simply due to stochastic reexamination, but only three gave replicable results (SI Ap- sampling effects (P value = 0.098). However, if E was actually pendix,DatasetS1).However,theseincludedtheindividualbone present at only 10% in the ancient Pacific chickens, then it is reported to have generated the critical single haplogroup E se- also highly unlikely that haplogroup E sequences would have quence(H268ofouruniquehaploptypes)usedtolinkRapaNui been detected in 15/31 (48%) of the specimens in previous andSouthAmerica(samplePAQANA011;SIAppendix,Fig.S8) studies (P value= 6.9 ×10−9). (14).Indirect contrast to the previous results, ourreanalysis of A median-joining network of the haplogroup D chicken se- an independent sample of PAQANA011 using Shrimp DNase quences revealed that all of the ancient Pacific sequences gen- PCRpretreatmentyieldedahaplogroupDsequence(haplotype eratedinthisstudy(n=22)andthosefrompreviousstudies(n= H239; SI Appendix, Dataset S6) identical to those of the other 16) (14, 16) together comprise only five different haplotypes two Rapa Nui specimens we reexamined. This result was sub- (Fig. 3), none of which have been found outside the Pacific re- sequently confirmed through independent replication of a gion.FourofthesefivearefromPolynesiaandclustertogether, subsample of the same specimen at Durham University (SI Ap- possessingadiagnosticmotifoffourSNPs(A→Gatbase281, pendix,Dataset S3). C→Tatbase296,T→Catbase306,A→ Gatbase342com- Our results further revealed that the PAQANA011 specimen pared with NC_007235; SI Appendix, Dataset S5). The four di- contained low amounts of DNA, with elevated levels of DNA agnosticSNPswerealsodetectedinfouradditionalhaplotypes template damage (SI Appendix, Dataset S4), and strongly sug- within the diversity of sequences from modern chickens sam- gests the previously reported haplogroup E sequence was the pled across the western Pacific and the Philippines, but only resultofcontaminationwithmodernchickenDNA.Afurther10 from Vanuatu, Santa Cruz, Philippines, and Guam (Fig. 3) samples excavated from the same site on Rapa Nui (Anakena) (26).Indeedoneof thepreviouslypublished WMGs,fromthe were also examined, and all yielded replicable haplogroup D Philippines (NC_007236; 25), contains all four of these di- sequences(haplotypeH239;SIAppendix,DatasetS1).Together agnosticSNPs(SIAppendix,Figs.S1andS2).Fig.3showsthat with the haplogroup D results of the previous study (14), this the most common ancient haplotype, H239, forms the central meansthatall15differentbonesexaminedattheAnakenasite node from which the other three ancient Polynesian D hap- haveyieldedH239sequences. lotypes radiate, consistent with a recent rapid expansion. The To investigate the conflict between the results obtained here central haplotype was also the most common sequence in and those previously reported from ancient Pacific specimens modernPacificchickenpopulations,beingpresentonalmostall (14–16), we calculated the probability of detecting the reported Pacificislands sampled. proportionsofDandEhaplogroupsgiventhedifferentdatasets. IfhaplogroupEwasauthenticallypresentwithinancientPacific SouthAmerica.Giventhatatleastsomeofthe previouslyrepor- chickens at the levels previously reported (48%) (14–16), then ted ancient Pacific chicken data appear to be due to contami- theprobabilitythatall22ofourancientsampleswouldbelongto nation, and the fact that all of the authenticated or reliable haplogroup D is negligible (P value = 1.3 × 10−7). In contrast, ancient Pacific chicken sequences are restricted to the unique 4828 | www.pnas.org/cgi/doi/10.1073/pnas.1320412111 Thomsonetal. H1*46 Node size proportional to number of samples H126** 150 50 25 10 51 H038* H1*5*4 H268 H12H91*71H**085**H224** H053*H***189H188H**1**9**0 HH01*9*7*5*3* H*1*2*3 HH*01*95*40**H**0*7*7 H*1*5*1 (hapHlo0g0r2oup E) H1*8*5* H186**H*187*** H153*** H264** H035** H132* H*0*32 H*0*3*4 H*1*7*4 H*0*1*0 H076*** H144* H*0*1*6 H*0*7*9 H266*** 9H*0*0*9 LegenRdapa Nui China H*1*24H*2*65 H*0*3*3 H**2*3*9 H*2*6**3 H131*** H*1*7*2 Ancient VGHNPFMaaeuaiauinrwsarueuqma a(uihtieuisstaosric) MVMKMThieeyaanaatdlaiynnlawaaamgnmiadasrcar Fustisugind.gy3,.2i0nU1cnlubrdopiontogefd1mh4ia4tpoclsohegoqrnuodeurnpicaeDlsCnRgeetawnneoarralkytzegedednieinnrattthheiidss H*2**60H*2H**7*2*17*4 H269****H273****H**2*5*8 H**2*6*1H**2*6H**2*2*5*9 P‘aonlycneestsriaaln’ Modern ISEA SSJPPInaaoahdpnplioloatiumpanna pCe oNisrnniueae Iwzss. GuineaSZPrui*em*ds*e*abn 4cna ePb oSNfw SPNesPs s2nCt6ouu,lmdo2yrb7sea,rrn4es3df–lfe1r4oc,62tm)2.s6aLtmahsbeipsqelliussntegounndcloyencsoa(SftdrIioeoAsmnpa,pGrweeeintunhdnBiioaxqnu,uktDeli(anh1te4aa,psae2lroto1t–uyS2np64d)e,. *** 3 PSNPs Line le1nbgpth proportional to sequence differenHc*2e*7*s2 H**2*2*5 Chicken ISnrdi Liaanka *** 12 PPSSNNPPss pmieodecrhnarPtascifriecp(rwehseitnet)i,nagndamncoiednetrnsIaSmEApl(egsra(yb).lack), GY O L O P O Pacific group of haplogroup D sequences, we performed ap- around 3,850 years ago (ya), without further onward trans- R proximate Bayesian coalescent simulations to evaluate the evi- portation of chickens into Western and Eastern Polynesia (SI NTH dence for the pre-Columbian introduction of chickens to South Appendix, Fig. S12). In contrast, the origins of the chickens A America. The coalescent simulations provided no evidence to currently found in Polynesia appear to be via the standard supportprehistoricdispersalofchickensfromPolynesiatoSouth southern route from New Guinea to the Solomon Islands, the Americaeitherwhenthedatasetsincluded(i)ancientsequences Santa Cruz Islands, Vanuatu, and further eastward (Fig. 1, only from haplogroup D or (ii) all sequences reported from arrows1 and4–7,andSIAppendix,Fig.S12). ancient specimens (Haplogroups B, E, and D) (14–16) (SI Ap- pendix,Fig.S11andTablesS2–S5).Theanalysesrevealthatthe Discussion more likely route and explanation for South American chicken OurresultsindicatethatasmallclusterofmtDNAhaplogroupD diversity appears to be via Europe and early historical intro- sequences,definedbyadiagnosticcombinationoffourCRSNPs ductions, or as modern DNAcontamination of experiments (SI (whichwetermthe“ancestralPolynesianmotif”),representthe Appendix,Fig.S11).AsingleDhaplotypesequence(H033)has founding lineages of chickens transported as prehistoric domes- been reported from post-European contact Peru (16), but this ticates across the Pacific and ultimately ending up in Polynesia sequence is common within ISEA populations, and could have (i.e., “Polynesian chickens”). We suggest that the most common beenassociatedwithearly-colonialSpanishtrade.Importantly,it haplotypeinancientsamples(H239)representsthecoremtDNA hasnotbeenfoundamongtheancientPacificchickensequences. lineage of Polynesian chickens, and that the one- or two-step derivativesinancientPacificislandspecimens(Fig.3)representin Micronesia and Western Polynesia. To investigate early human- situevolutionfollowingcolonization.Thishypothesisissupported mediated dispersal patterns within the Asia–Pacific region, we by the geographic distribution of the ancient daughter lineages, examined modern chickens from islands across ISEA, Micro- which are unique to each Pacific island group, and the elevated nesia,andWesternPolynesia,becausefewspecimensofancient frequencyoflineageswiththefourdiagnosticSNPsintheeastern chickens were available from this area [however, see Fais D Pacific (SI Appendix, Fig. S10). Although mtDNA is maternally haplotype sample from Storey et al. (16)]. Although the ISEA inherited as a single genetic locus, limiting the ability to recover sequences are scattered across the haplogroup D network, the complexcolonizationhistories,ourdataestablishclearhypotheses majority of haplotypes from modern Pacific chickens are ge- thatcanbetestedwithgenomicdatafrombothmodernandan- netically clustered together (H032–35, H085, H224–225, H260, cient chickens, and other groups such as humans, commensals, H262,H271–274;Fig.3).WithinMicronesia,haplogroupDhas andother domesticates. Itisimportant to note thatinsituations been reported from modern chickens in Guam (n = 3/5; H032, like the Pacific, phylogeographic signals in domestic species are H224,and H225; 26), although interestingly, these particular D likely to represent processes of initial human dispersal and later haplotypesarenotsharedwithanyotherPacificislandgroup.In tradepatterns. fact,twoofthesehaplotypeshaveonlyeverbeenfoundinGuam Ourfindingscontrastsubstantiallywithpreviousstudies(14–16), (H224 and H225),whereas the third Guam haplotype is shared whichwesuggeststemsfromourstrictadherencetocontamination with the Philippines, Japan, Indonesia, and Papua New Guinea reduction measures—for example, the use of Shrimp DNase. By (H032). The modern haplogroup D chickens in Guam do not removing a key source of potential contamination with domestic appeartobesignificantlygeneticallydifferentiatedfromthosein chicken DNA (PCR reagents), the use of Shrimp DNase has thePhilippines, Japan,andIndonesia(SIAppendix,TableS6). allowed us to recharacterize the crucial ancient Rapa Nui sample An investigation of the discordant haplogroup D lineages in from a prior study (PAQANA011) as haplogroup D and not, as MicronesiaandPolynesiausingcoalescentsimulationsidentified previouslyreported,haplogroupE.Consequently,wecastdoubton an early movement of chickens between New Guinea and the authenticity of other haplogroup E sequences reported from Micronesiaasthemostlikelyoffivemodelstested(SIAppendix, ancientPacificchickenspecimens,wheresuchprocedureswerenot Fig. S12 and Tables S7 and S8). The simulations suggest that used.Perhapsmoreimportantly,wesuggestitwillbeverydifficult chickensweretransportedbetweenMicronesiaandislandsinthe tocategoricallyruleoutcontaminationasthesourceofhaplogroup Bismarck Sea off the coast of New Guinea and New Britain E sequences in ancient samples, due to the sporadic presence of Thomsonetal. PNAS | April1,2014 | vol.111 | no.13 | 4829 domesticate DNA in laboratory consumables (34) and the likeli- chicken populations, from the Santa Cruz Islands, Solomon Is- hood that any such contamination would result in haplogroup E lands,andVanuatu(26).Therefore,Polynesianchickensmaybe sequences. Importantly, sequencing longer stretches of such con- oneofthefewexampleswhereancestralgeneticpatternscanstill taminatingtemplates(17)doesnotprovideanyadditionalsupport beobserved inadomesticated species.Chickens onremote Pa- forauthenticity. cific islands may also contain Polynesian nuclear genomic line- Our recharacterization of the Rapa Nui PAQANA011 speci- ages,andifso,wouldrepresentoneofthefewsurvivingexamples menashaplogroupDhasimplicationsfortheotherEsequences ofprecolonialdomesticchickens. reported by Storey et al. (14), including the putative ancient Conclusion Chilean chicken sequence from El Arenal-1 used to propose a prehistoric link between Polynesia and South America. Co- AlthoughmtDNAlacksthepowerofgenomiclocitoreconstruct alescent simulations using “all ancient haplogroups” and the complex evolutionary histories,we showthat an informative re- modern data found that a European–South America route was gion ofthe chicken mitochondrial genome can beused totrace more likely than a direct link between haplogroup E chicken theirhumandispersalinthePacific.Theanalysisofancientand sequences in Polynesia and South America, due to the phylo- modernspecimensrevealsauniquePolynesiangeneticsignature, geographic signals within the worldwide dataset showing more whichcanbetracedbacktoISEA,andpromisestoallowfurther similarities between chickens from Europe and South America. resolution of migration and trading routes in the area. Impor- Perhapsmoregenerally,thesefindingshighlighthowhaplogroup tantly,we revealthat apreviously reported connection between E sequences are uninformative in nature and lack phylogeo- pre-European South America and Polynesian chickens most graphic signal worldwide. A clear understanding of the nature likelyresultedfromcontaminationwithmodernDNA,andthat andextentofPolynesiancontactwithSouthAmericawillrequire this issue is likely to confound ancient DNA studies involving genomic analyses of both ancient and modern populations of haplogroup E chicken sequences. These observations reaffirm humans,commensals, anddomesticates. the potential of coalescent simulations of genetic data to eval- ThedistributionofthenineDhaplotypescurrentlyknownto uate new hypotheses regarding the dispersal of humans, com- share the ancestral motif provides a unique genetic signature mensals, and domesticates derived from archaeology. These that can be used to trace the human dispersal of chickens hypothesescanbefurthergroundedusinggenomic-scalestudies throughISEAandthePacificislands.Ourreconstructionofthe in combination with direct dating and genetic investigation of chicken colonization history of Micronesia highlights how sim- newarchaeologicalsamples. ulations with CR data can provide sufficient phylogeographic signaltogeneratenewhypothesesregardingtradeandmigration MaterialsandMethods scenarios.Althoughithasbeenproposedthatmanycommensals Samples. Thirty-seven ancient chicken bones were collected for analysis, and domesticates are late arrivals to the Micronesian islands comprisingeightfromNiue,11fromHawai’i,and18fromRapaNuiexca- compared with humans (35), we have reconstructed a link be- vated from deposits at Anakena by T.L.H. [including the six samples pre- tweenchickensfromislandsintheBismarckSeaandMicronesia viously analyzed by Storey et al. (14); SI Appendix, Dataset S1]. Modern thatdatesto∼3,850B.P.Suchanearlydateisbroadlyconsistent feathersamplesfromISEAandthePacific(n=124)werealsoexaminedto with archaeological evidence for human settlement of Saipan investigatecurrentphylogeographicpatterns(forlocationdetails,seeFig.1 at 3,300–3,500 B.P. (36) and Palau at almost 4,000 B.P. (35), andSIAppendix,Figs.S3–S6,TableS1,andDatasetS6).Theancientsamples however few comparably early zoo-archaeological remains have wereextracted,amplified(usingprimersinSIAppendix,Fig.S13),andse- beenfoundinMicronesiatodate(10,13,37).Theinferredlink quenced at the Australian Centre for Ancient DNA (ACAD) in Adelaide, between chickens from the Bismarcks and Micronesia without SouthAustralia,accordingtoarangeofstrictprotocols(39),includingnu- subsequent eastward movement does not support a two-wave merouscontrols.Importantly,weincludedShrimpDNasepretreatmentinall PCRreactions,beforeaddingtemplateDNA,toremoveanycontaminating modelofPolynesianorigins(14,15,38)whereanearlierLapita migrationwave(2,800–3,500ya)wasmixedwithasecond, later double-strandedDNAintroducedviaPCRreagentsandplastic-ware(SIAp- wave moving through Micronesia to Western Polynesia (1,500– pendix)(33).Independentexternalreplicationwithdirectsequencingofthe PAQANA011 ancient sample was performed in a dedicated ancient DNA 2,000 ya). Our simulations suggest that there was little in- laboratoryintheArchaeologyDepartmentatDurhamUniversityfollowing teraction between chickens from Micronesia and the islands strict laboratory procedures (39). Theinitial and independently replicated further eastward. One caveat concerning the power of the sim- PCR fragments from bone sample PAQANA011 were also cloned and se- ulationanalysisisthesmallnumberofMicronesiansamples[one quencedattheACADlaboratories(SIAppendix,DatasetS4).Modernsam- ancientFais(16)andfivemodernGuam(26)specimens]andthe pleswereextracted,withthehighlyvariable201bpoftheCRamplifiedand expected historical and recent turnover of chicken populations sequencedinaphysicallyseparatepre-PCRcleanlaboratoryattheUniversity in the region. Reassuringly, the ancient Fais haplotype H260 is of Adelaide and in the Archaeology Department at Durham University, present in modern chickens from the Santa Cruz (n = 2) and followingstandardprotocols(39). Solomon Islands (n = 5), apparently surviving any later in- trogression. Our reconstruction of the colonization history of WMGAnalysis.Todeterminetherobustnessofthecurrentstandardchicken Micronesian chickens demonstrates the potential power of co- phylogeneticframeworkfortheanalysisoftheshortancientsequences,all61 alescent simulations to test hypothesized migration and trade WMG sequences (25) were downloaded and aligned; PartitionFinder (40) routesinarchaeologyandanthropology. was used to identify the number of preferred partitions and their sub- The only ISEA location where the ancestral SNP motif has stitutionmodel;andphylogenetictreeswereproducedusingbothBayesian beendetectedareCamiguinandManilainthePhilippines, and (MrBayesv3.2;41)andmaximumlikelihoodestimation(RaxMLv7.0.4;42). a link with this area is consistent with other lines of evidence SeeSIAppendixformoredetails. about early Polynesian origins (3, 4, 31). The other Philippine chicken haplotypes are spread throughout the haplogroup D CRSequenceAnalysis.Inadditiontothe144CRsequencesgeneratedinthis study,wedownloaded1,226worldwidemtDNACRchickensequencesfrom network (Fig. 3), reflecting relatively high genetic diversity (haplotype diversity=0.89;SIAppendix,TableS9). GenBank to establish the geographic distribution for each chicken hap- logroup(14,21–24,26,27,43–46).ToallowdirectcomparisonsoftheCR Despite extensive European settlement in the Pacific region haplotypes,the1,370chickensequenceswerealignedandtrimmedtothe overthelast fewcenturies,manynativechickenpopulationsap- highlyvariable201bpcommontoallofour144newlygeneratedsequences peartocontainrelativelyhighfrequenciesoffoundingmitochondrial (referredtoas“201bpCRdataset”).The201bpCRdatasetwascollapsedto lineages—forexample,theMarquesas,SolomonIslands,Vanuatu uniquehaplotypesusingCollapsev1.2,resultingin274uniquehaplotypes (26), and the Santa Cruz Islands—suggesting a high level of ge- (H001–H274;SIAppendix,DatasetS6;referredtoas“uniqueCRhaplotype netic continuity on these islands since prehistoric times. In addi- dataset”). ModelGenerator (47) was used to establish the best model to tion to the two ancient haplotypes detected in modern samples, fittheuniqueCRhaplotypedataset(GTR+I+G).Thehaplogroupofeachof many other D haplotypes are also present in modern Pacific our 144 newly generated sequences was established by comparison with 4830 | www.pnas.org/cgi/doi/10.1073/pnas.1320412111 Thomsonetal. sequences of known haplogroup designation from Liu et al. (24) (SIAp- tested whether coalescent simulations and approximate Bayesian compu- pendix, Dataset S6). As the majority of the new 144 CR sequences were tationofthe201bpCRdatasetcouldreconstructaprehistoriclinkbetween identifiedashaplogroupD,aMedianJoiningNetwork(usingNetworkv4.6; thePacificandSouthAmerica(SIAppendix).Toexplorelikelydemographic 48)wasgeneratedforjusttheDhaplogroup(SIAppendix).Allnewsequences historiesforchickensinMicronesiaandPolynesia,wealsousedBayeSSCto wereuploadedtoGenBank(KJ000585–KJ000642;SIAppendix,DatasetS6). simulatealternatehypothesesofmigrationroutesforcomparisonwiththe observedphylogeographicpatternswithinthePacific. StatisticalAnalysis.Toexaminethediscrepanciesbetweenthecompositionand phylogeographicdistributionofhaplogroupsreportedbyStoreyetal.(14,16) ACKNOWLEDGMENTS. We thank Jessica Metcalf, Peggy Macqueen, and andthosegeneratedinthisstudy,wetestedthelikelihoodofdetectingthe other members of the Australian Centre for Ancient DNA for assistance; reportedproportionsunderdifferentscenarios.Alinearregressionplotwasalso John Terrell for manuscript discussions; Richard Walter (Department generatedtovisualizethecorrelationbetweenoccurrenceofthefourcharac- of Anthropology, University of Otago) and Atholl Anderson (Australian teristic CR SNPs of the Polynesian chicken and longitude using the standard National University) for providing access to the Niue samples; and Will MillardforcollectingmodernfeathersamplesfromNewGuinea.G.L.and plottingfunctioninR. K.D. also thank Atholl Anderson, Hanneke Boon, James Wharram, Klaus Hymphendahl, Matt Fletcher, Ingo Isensee, and the Lapita Expedition for BayesianCoalescentSimulations.Given the importance of pre- and post- collectingfeathersamplesfromISEAandWesternPolynesia.Thisworkwas ColumbianmtDNAsequencesfromChileandPeru,respectively(14,16),we fundedbytheAustralianResearchCouncilandtheUniversityofAdelaide. 1. DenhamT,etal.(2012)DatingtheappearanceofLapitapotteryintheBismarck 26.Dancause KN, Vilar MG, Steffy R,Lum JK (2011) Characterizing genetic diversity Archipelago and its dispersal to Remote Oceania. Archaeology in Oceania 47(1): ofcontemporarypacificchickensusingmitochondrialDNAanalyses.PLoSONE 39–46. 6(2):e16843. 2. GreenRC(1991)TheLapitaCulturalComplex:Currentevidenceandproposedmodels. 27.GongoraJ,etal.(2008)Indo-EuropeanandAsianoriginsforChileanandPacific Indo-PacificPrehistoryAssociationBulletin11(1991):295–305. chickensrevealedbymtDNA.ProcNatlAcadSciUSA105(30):10308–10313. 3. KirchP(1997)TheLapitaPeoples:AncestorsoftheOceanicWorld(Blackwell,Oxford). 28.GongoraJ,etal.(2008)ReplytoStoreyetal.:MoreDNAanddatingstudiesneeded 4. SpriggsM(1997)TheIslandMelanesians(Blackwell,Oxford). forancientElArenal-1chickens.ProcNatlAcadSciUSA105(48):E100. 5. SummerhayesG(2000)LapitaInteraction(AustralianNationalUnivPress,Canberra, 29.StoreyAA,etal.(2008)Pre-Columbianchickens,dates,isotopes,andmtDNA.Proc Australia). NatlAcadSciUSA105(48):E99,authorreplyE100. Y 6. WilmshurstJM,HuntTL,LipoCP,AndersonAJ(2011)High-precisionradiocarbon 30.Gonçalves VF, et al. (2013) Identification of Polynesian mtDNA haplogroups in G O dAacatidngScsihUowSAsr1e0c8e(n5t):1a8n1d5r–a1p8i2d0i.nitialhumancolonizationofEastPolynesia.ProcNatl r6e4m65a–in6s46o9f. Botocudo Amerindians from Brazil. Proc Natl Acad Sci USA 110(16): OPOL 7. ClarkeAC,BurtenshawMK,McLenachanPA,EricksonDL,PennyD;SMBETri-National 31.BellwoodP(2007)PrehistoryoftheIndo-MalaysianArchipelago(AustralianNational HR Young Investigators (2006) Proceedings of the SMBE Tri-National Young Inves- UnivPress,Canberra,Australia),3rdEd,p385. NT tigators’Workshop2005.ReconstructingtheoriginsanddispersalofthePolynesian 32.BarkerFK,BeneshMK,VandergonAJ,LanyonSM(2012)Contrastingevolutionary A bottlegourd(Lagenariasiceraria).MolBiolEvol23(5):893–900. dynamicsandinformationcontentoftheavianmitochondrialcontrolregionand 8. WhitakerT,CarterG(1954)Oceanicdriftofgourds—Experimentalobservations.AmJ ND2gene.PLoSONE7(10):e46403. Bot41(9):697–700. 33.ChamplotS,etal.(2010)AnefficientmultistrategyDNAdecontaminationprocedure 9. RoullierC,BenoitL,McKeyDB,LebotV(2013)Historicalcollectionsrevealpatternsof ofPCRreagentsforhypersensitivePCRapplications.PLoSONE5(9):e13042. diffusionofsweetpotatoinOceaniaobscuredbymodernplantmovementsandre- 34.LeonardJ,etal.(2007)AnimalDNAinPCRreagentsplaguesancientDNAresearch. combination.ProcNatlAcadSciUSA110(6):2205–2210. JArchaeolSci34(9):1361–1366. 10.LarsonG,etal.(2007)PhylogenyandancientDNAofSusprovidesinsightsinto 35.WicklerS(2004)ModellingcolonisationandmigrationinMicronesiafromazooarch- neolithicexpansioninIslandSoutheastAsiaandOceania.ProcNatlAcadSciUSA aeologicalperspective.Colonisation,Migration,andMarginalAreas:AZooarchaeo- 104(12):4834–4839. logicalApproach,edsMondiniM,MunozS,WicklerS(OxbowBooks,Oxford),pp 11.Larson G, et al. (2010) Patterns of East Asian pig domestication, migration, and 28–40. turnoverrevealedbymodernandancientDNA.ProcNatlAcadSciUSA107(17): 36.CarsonM(2008)RefiningearliestsettlementinRemoteOceania:Renewedarchae- 7686–7691. ologicalinvestigationatUnaiBapot,Saipan.JournalofIsland&CoastalArchaeology 12.SavolainenP,LeitnerT,WiltonAN,Matisoo-SmithE,LundebergJ(2004)Adetailed 3(1):115–139. pictureoftheoriginoftheAustraliandingo,obtainedfromthestudyofmitochon- drialDNA.ProcNatlAcadSciUSA101(33):12387–12390. 37.Anderson A(2009) Theratandthe octopus:Initialhumancolonizationandthe prehistoricintroductionofdomesticanimalstoRemoteOceania.BiolInvasions11(7): 13.Matisoo-SmithE,RobinsJH(2004)OriginsanddispersalsofPacificpeoples:Evidence 1503–1519. from mtDNA phylogenies of the Pacific rat. Proc Natl Acad Sci USA 101(24): 9167–9172. 38.AddisonDJ,Matisoo-SmithE(2010)RethinkingPolynesiansorigins:AWest-Polynesia Triple-Imodel.ArchaeologyinOceania45(1):1–12. 14.Storey AA, et al. (2007) Radiocarbon and DNA evidence for a pre-Columbian 39.Cooper A, Poinar HN (2000) Ancient DNA: Do it right or not at all. Science introduction of Polynesian chickens to Chile. Proc Natl Acad Sci USA 104(25): 10335–10339. 289(5482):1139. 40.LanfearR,CalcottB,HoSY,GuindonS(2012)Partitionfinder:Combinedselectionof 15.Storey AA, et al. (2010) Mitochondrial DNA from 3000-year old chickens at the Teoumasite,Vanuatu.JArchaeolSci37(10):2459–2468. partitioningschemesandsubstitutionmodelsforphylogeneticanalyses.MolBiolEvol 16.StoreyAA,etal.(2012)Investigatingtheglobaldispersalofchickensinprehistory 29(6):1695–1701. usingancientmitochondrialDNAsignatures.PLoSONE7(7):e39171. 41.RonquistF,etal.(2012)MrBayes3.2:EfficientBayesianphylogeneticinferenceand 17.StoreyAA,etal.(2013)PolynesianchickensintheNewWorld:Adetailedapplication modelchoiceacrossalargemodelspace.SystBiol61(3):539–542. ofacommensalapproach.ArchaeologyinOceania48(2):101–119. 42.StamatakisA(2006)RAxML-VI-HPC:Maximumlikelihood-basedphylogeneticanaly- 18.Oskarsson MCR, et al. (2012) Mitochondrial DNA data indicate an introduction seswiththousandsoftaxaandmixedmodels.Bioinformatics22(21):2688–2690. throughMainlandSoutheastAsiaforAustraliandingoesandPolynesiandomestic 43.AdebamboA,etal.(2010)LackofphylogeographicstructureinNigerianvillage dogs.ProcBiolSci279(1730):967–974. chickensrevealedbymitochondrialDNAD-loopsequenceanalysis.IntJPoultSci9(5): 19.MuirWM,etal.(2008)Reviewoftheinitialvalidationandcharacterizationofa3K 503–507. chickenSNParray.WorldsPoultSciJ64(2):219–225. 44.Berthouly-SalazarC,etal.(2010)Vietnamesechickens:AgatetowardsAsiangenetic 20.FumihitoA,etal.(1996)Monophyleticoriginanduniquedispersalpatternsofdo- diversity.BMCGenet11:53. mesticfowls.ProcNatlAcadSciUSA93(13):6792–6795. 45.MuchadeyiFC,etal.(2008)MitochondrialDNAD-loopsequencessuggestaSoutheast 21.KanginakudruS,MettaM,JakatiRD,NagarajuJ(2008)GeneticevidencefromIndian AsianandIndianoriginofZimbabweanvillagechickens.AnimGenet39(6):615–622. redjunglefowlcorroboratesmultipledomesticationofmoderndaychicken.BMC 46.MwacharoJM,etal.(2011)MitochondrialDNArevealsmultipleintroductionsof EvolBiol8:174. domesticchickeninEastAfrica.MolPhylogenetEvol58(2):374–382. 22.OkaT,etal.(2007)AnalysisofmtDNAsequencesshowsJapanesenativechickens 47.KeaneTM,CreeveyCJ,PentonyMM,NaughtonTJ,MclnerneyJO(2006)Assessment havemultipleorigins.AnimGenet38(3):287–293. ofmethodsforaminoacidmatrixselectionandtheiruseonempiricaldatashowsthat 23.SilvaP,etal.(2009)MitochondrialDNA-basedanalysisofgeneticvariationandre- adhocassumptionsforchoiceofmatrixarenotjustified.BMCEvolBiol6(29):29. latednessamongSriLankanindigenouschickensandtheCeylonjunglefowl(Gallus 48.Bandelt HJ, Forster P, Röhl A (1999) Median-joining networks for inferring in- lafayetti).AnimGenet40(1):1–9. traspecificphylogenies.MolBiolEvol16(1):37–48. 24.LiuYP,etal.(2006)Multiplematernaloriginsofchickens:OutoftheAsianjungles. 49.SpechtJ(2007)Smallislandsinthebigpicture:TheformativeperiodofLapitainthe MolPhylogenetEvol38(1):12–19. BismarckArchipelago.OceanicExplorations:LapitaandWesternPacificSettlement, 25.MiaoYW,etal.(2013)Chickendomestication:Anupdatedperspectivebasedon TerraAustralis,edsBedfordS,SandC,ConnaughtonSP(AustralianNationalUniv mitochondrialgenomes.Heredity(Edinb)110(3):277–282. Press,Canberra,Australia),Vol26,pp51–70. Thomsonetal. PNAS | April1,2014 | vol.111 | no.13 | 4831 Supplementary Materials and Methods Sample collection. Thirty-seven ancient chicken bones were collected for analysis, comprising: eight ancient chicken bones from archaeological sites at Paluki and Anatoloa in Niue; 11 ancient Hawaiian chicken bones from an excavation at Makauwahi Cave on Kauai, Hawaii collected by DB; and 18 Rapa Nui chicken bones excavated from deposits at Anakena collected by TH. The 18 Rapa Nui bones include the six samples previously analyzed by Storey et al. (1) (Table S1). One hundred and twenty four modern feather samples were also examined to investigate recent phylogeographic patterns. These included 107 modern feathers from ISEA and Remote Oceania collected by GL and KD in 2008 and 2009: 28 from the Santa Cruz Islands, 31 from the Solomon Islands, 13 from Papua New Guinea, 10 from Indonesia, 23 from the Philippines and two from Vietnam. An additional 17 naturally shed modern feather samples were collected from the Marquesas (French Polynesia, n=6) by TH, and Kokee, Kauai (Hawaii, n=11) by TH/DB. Details on the locations of these modern samples can be found in Table S1 and are shown in Fig. 1. Ancient DNA Extraction, Amplification, and Sequencing. The samples were extracted, amplified, and sequenced in specialist ancient DNA (aDNA) laboratories at the Australian Centre for Ancient DNA (ACAD) in Adelaide, South Australia, according to a range of strict protocols and including controls (2). Ancient bone samples (n=37) were extracted and PCR experiments set up in the physically remote ACAD ancient laboratory, whereas the feathers (n=17) were extracted and PCR experiments set up in the physically remote ACAD pre-PCR clean-room laboratory. Independent external replication of the ACAD9068 (PAQANA011) ancient sample was performed in a dedicated aDNA lab in the Archaeology Department at Durham University following strict laboratory procedures (2). ACAD ancient bone extractions. Each chicken bone was ground to fine powder in a Mikrodismembrator (5000 rpm, for 10 seconds). Approximately 70 mg of bone powder was decalcified concurrently with protein digestion by incubation at 55 °C overnight in 1mL of extraction buffer (consisting of 0.4725 M EDTA (pH=8.0), 0.2 % sodium dodecyl sulphate (SDS), and 0.7 mg.ml-1 Proteinase K). After digestion, samples were centrifuged at 10,000 rpm for 5 mins and the supernatant was transferred to an Amicon ultra-4 (Millipore), which was centrifuged at 4000 xg until only 100 µL supernatant remained. The supernatant was washed with 1 mL molecular grade water and centrifuged again (at 4000 xg until only 100 µl remained). An equal volume of ATL buffer (Qiagen DNeasy kit) was then added, mixed, and the supernatant removed to a 2 mL screw-cap tube. The supernatant was incubated for 10–60 mins at room temperature on a rotary mixer after the addition of an equal volume of AL buffer (Qiagen DNeasy kit) and 0.02 µg.µl-1 of carrier RNA. After the incubation, an equal volume of ethanol (100 %) was added, and then the total volume was transferred to a Qiagen DNeasy spin column where it was incubated at room temperature for 10–60 mins. The extraction then followed the Qiagen DNeasy kit instructions, with the following exceptions at the elution stage: 100–150 µL of warmed AE buffer was added and then incubated at room temperature for 10–30 mins, before being centrifuged at 8,000 rpm for 1 min, this step was repeated to finish with 200–300 µL of total volume. ACAD PCR amplification and sequencing of ancient samples. A 330 base pair (bp) segment of the mtDNA CR was amplified and sequenced from each specimen in short overlapping fragments (Table S10, Fig. S13), which is necessary to ensure amplification of the short damaged fragments of ancient DNA samples. PCRs were set up using 25 µL volumes containing a final concentration of 1 U Platinum Taq DNA Polymerase High Fidelity (Invitrogen), 1 x PCR Buffer (Platinum, Invitrogen), 3 mM MgSO4, 200 µM each dNTP, 2 mg.ml-1 rabbit serum albumin (Sigma), 1 µM forward and reverse primers and 2-3 µl of template DNA. PCR reactions were performed on a Corbett Research Palm Cycler using the following cycling conditions: 94 °C for 2 min, 55 cycles of 94 °C for 30 s, 55 °C for 30 s, 68 °C for 30 s, and a final extension of 10 min at 68 °C. Amplifications of extraction and PCR controls were performed in all experiments to monitor contamination. PCR products were separated by electrophoresis on a 3.5 % agarose gel. Successful PCR products (10 µl) were purified using 0.8 µl of EXOSAP (Fermentas) at final concentration of 0.38 U/µl Exonuclease I, and 0.05 U/µl Shrimp Alkaline Phosphotase, and thermal cycled at 37 °C for 30 mins, 80 °C for 15 mins, and 15 °C for 3 mins on a Corbett Research Palm Cycler. The forward and reverse complements of each fragment were sequenced from the same PCR reaction using the same primers as for the PCR, and Big Dye Terminator v3.1 cycle-sequencing chemistry, followed by vacuum clean up on a Multiscreen® SEQ plate (Millipore). The sequencing run was conducted on an ABI 384 3130XC capillary sequencer. Primers GG144F/GG387R and GG316F/GG586R (1) were used initially to amplify a portion of the mitochondrial (mtDNA) control region but as the PCR products amplified from these primers (fragment 1 and 2) are 250bp and 305bp respectively; additional primers were designed to cover the same range of mtDNA control region. Primer GG144F was paired with A1781 (187bp as fragment 3) and A1780 was paired with GG387R (151bp as fragment 4) to cover the equivalent DNA sequence as fragment 1 but in two overlapping fragments (Table S10). Primers A1958 and A1959 (192bp as fragment 5) were used to cover the balance of the mtDNA CR under study for the ancient samples. The use of this alternative primer set meant that a sequence gap was introduced in some ancient sequences equivalent to the primer binding region (Fig. S13). Further trimming to the sequence length shared across all chicken specimens resulted in a final sequence length of 201bp. Durham Bone extractions as replication for PAQANA011. DNA extraction of the replicate ancient chicken bone fragment PAQANA011 was performed in a dedicated aDNA lab in the Archaeology department at Durham University following strict laboratory procedures as per commonly used guidelines (2). All equipment and work surfaces were cleaned before and after each use with a dilute solution of bleach (10 %) followed by ethanol (99 %). The ancient chicken bone (~0.05 g) was pulverized in a Micro-dismembrator, digested in 0.425 M EDTA, 0.05 % SDS, 0.05 M Tris- HCI and 0.333 mg.ml-1 proteinase K and incubated overnight on a rotary mixer at 50 oC until fully dissolved. 2 ml of solution was then concentrated in a Millipore Amicon Ultra-4 30 KDa MWCO to a final volume of 100 µl. The concentrated extract was purified using the QIAquick PCR Purification Kit following manufacturers recommendations, except that the final elution step was performed twice to produce a final volume of 100 µl. A negative extraction control was performed alongside the ancient bone sample. Durham PCR amplification and sequencing of ancient samples. PCRs were setup in 25 µl reactions using 1.25 U Taq GOLD (Applied Biosystems), 1 x Gold buffer (Applied Biosystems), 2.5 mM MgCl , 0.5 µg.µl-1 bovine serum albumin (BSA), 200 µM of each dNTP, 0.8 2 µM of each forward and reverse primers, and 2-5 µl of aDNA extract. We used PCR primers (5’-3’) GG144F and GG387R; GG316F, and GG586R (1). One PCR negative control was included for every three aDNA template PCR tubes. We ran a total of 22 PCRs with aDNA template, eight PCR negative controls and two PCR negative extraction control. Neither the PCR negative controls nor the negative extraction control produced bands (PCR product) when analyzed by gel-electrophoresis. PCR cycling conditions were 95oC for 5 min, 50 cycles of 94 °C for 45 sec, 54 °C for 45 sec and 72 °C for 45 sec, followed by 72 °C for 10 min. PCR products were stored at -20 °C. Sanger sequencing on the Applied Biosystems 3730 DNA Analyser was performed at the DNA sequencing service in the School of Biological and Biomedical Sciences at Durham University. Modern DNA Extraction, Amplification, and Sequencing. ACAD modern feather extractions. Approximately 5 mm of each feather tip was rehydrated overnight with 1 ml phosphate buffered saline (PBS) on a rotary mixer at room temperature. On day 2, the supernatant was removed, the feather tip was macerated using a clean scalpel blade, and the sample was digested in 440 µl of digestion buffer (comprising ATL buffer (Qiagen DNeasy kit) with 1.8 mg.ml-1 Proteinase K, and 90 mM Dithiothreitol) overnight at 55 °C on a rotary mixer. After digestion, 400 µL of AL buffer (Qiagen DNeasy kit) and 0.02 µg.µl-1 of carrier RNA was added and incubated at room temperature on a rotary mixer for 10–30 mins, after which 400 µL of 100 % ethanol was added. The supernatant (650 µl) was incubated on a Qiagen DNeasy spin column for 10–30 mins before being centrifuged at 8000rpm for 1 min. This incubation was then repeated until all of the supernatant had been centrifuged through the column. The feather extraction protocol then followed that of the bone extraction procedure above. ACAD PCR amplification and sequencing of modern feather samples. PCR amplifications and sequencing of the 2 overlapping fragments were performed as per the ancient bone samples (see above). Durham modern feather extraction At Durham University, modern feathers from ISEA and Near Oceania were extracted in a pre-PCR clean room after Cooper & Poinar (2), using a protocol designed by Pfeiffer et al. (3) alongside the QIAquick PCR purification Kit (QIAGEN Ltd, UK). The tip of each feather was sampled (approximately 1cm cut into smaller fragments) and digested in 340µl extraction buffer containing 100mM Tris-HCl, pH8, 100mM NaCl, 3mM CaCl 2% SDS (w/v), 40mM DTT and 250µg/ml proteinase K following the protocol 2, by Pfeiffer et al. (3). The samples were incubated overnight at 56°C on a rotary mixer. Following digestion, the samples were purified using the QIAquick PCR purification Kit (QIAGEN Ltd, UK) following the manufacturers’ instructions. An extraction control was used for every run of seven samples. The quantity of DNA present within each extract was measured using the Quant-iT HS Assay Kit (Invitrogen) used with the Qubit fluorometer following the manufacturers’ instructions. Durham PCR amplification and sequencing of modern feather samples. The amplification of a 201bp fragment of the CR (a subset of the 330bp amplified from the ancient samples) was undertaken through PCR in a physically separated clean laboratory. The forward primer GG144F and the reverse primer GG387R (see Table S10) were used to amplify this 201bp fragment (excluding primers). The PCR amplifications were performed in a 25µl reaction mix containing 1µl of extract, 0.96x PCR Gold Buffer, 2.4mM MgCl , 1.2U Taq, 0.24mM dNTP and 0.96µM of each primer. The 2 PCR thermal cycling reactions consisted of 90s initial denaturation step at 94°C, followed by 35 cycles of 30s denaturation at 94°C, 30s annealing at 54°C, 30s extension at 72°C then a 10 minute final extension step at 72°C. The PCR products were visualized on a 0.5x agarose gel. Sequencing was performed on an ABI 3730 sequencer in the DNA-dedicated laboratory of the School of Biological and Biomedical Sciences. Cloning of PAQANA011 at ACAD. The PCR products generated from bone sample PAQANA011 were cloned using Stratagene and/or Topo cloning kits using manufacturers instructions (after an A-tailing reaction). The A-tailing reaction consisted of a 20 µl volume reaction containing 0.125 U HotMaster Taq, 2.5 µM dATP, 10x HotMaster buffer, 17 µl cleaned PCR products. The Buffer, dATP’s and Taq were activated at 94 °C for 2 mins prior to addition of the PCR products then a further incubation at 72 °C for 10 mins. The A-tailed PCR products were then cleaned up using an Isopropanol precipitation and resuspended in 10 µl of PCR grade water. Sanger sequencing of the cloned PCR products were performed according to the procedures outlined above. Phylogenetic inference WMG: To determine the robustness of the current phylogenetic framework used for chicken research, the 61 WMG sequences from that study were downloaded and aligned using the Muscle algorithm in Geneious v5.6 (4). PartitionFinder v1.0.1 (5) was used to identify the number of preferred partitions and their substitution model (CR with HKY plus Gamma; codon 1, codon 2 and tRNA with HKY; and codon 3 with GTR). MrBayes v3.2 was used to generate a phylogenetic tree using four runs of four independent chains of 100 million iterations, less 25% as burnin (6). Tests for convergence to stationarity were performed by analyzing the standard deviation of split frequencies (< 0.01). RaxML v7.0.4 was used to generate a maximum likelihood tree with the same partitions as above, with bootstrapping performed via 100 iterations followed by an optimized maximum likelihood search (7). To establish the level of phylogenetic concordance between topologies produced by WMGs versus the highly variable 201bp of the CR, the WMG data was split into two subsets, the 201bp fragment of the CR and the WMG excluding all of the CR. Each subset was rerun for the PartitionFinder and MrBayes analyses separately (i.e. the CR was run separately from the WMG data minus the CR), using the same parameters as above except only 2 million iterations were required to obtain a standard deviation of less than 0.01 for the four chains. mtDNA CR: In addition to the 144 CR sequences generated in this study, we downloaded 1226 worldwide mtDNA CR chicken sequences from Genbank (1, 8-17) to establish the geographic distribution for each chicken haplogroup (n=1370). Although additional CR sequences have since been uploaded to Genbank (total chicken CR sequences are currently >3000), overall haplogroup designations are not changed with the inclusion of additional sequences (18). To allow direct comparisons of the CR haplotypes, the 1370 chicken sequences were aligned and trimmed to the 201bp common to our 144 newly generated sequences (referred to as ‘full CR dataset’), with any indels removed. The 201bp hypervariable fragment is a useful region for reconstructing recent evolutionary events when DNA template length is a constraint (19, 20), such as in ancient DNA studies. For ease and clarity, the 1370 CR sequences were collapsed to unique haplotypes using Collapse v1.2 with manual adjustments where missing data caused short sequences to be considered different haplotypes, resulting in 274 unique haplotypes (H001-H274, see Dataset S6; referred to as ‘unique CR haplotype dataset’). The haplogroup of each of our 144 newly generated sequences was established by comparison to sequences of known haplogroup designation from Liu et al. (13) (see Dataset S6). The phylogenetic robustness of the full 330bp length (both fragment 1 & 2) was investigated using PhyML (21) to establish that inclusion of additional length sequences did not change the haplogroup designation of the new sequences (Fig. S14), with ModelGenerator (22) used to establish the model of best fit. We also explored the unique CR haplotype dataset in SplitsTree4 (23), using the NeighborNet algorithm, and found that the data appeared not to be tree-like, probably due to saturation and substitution rate heterogeneity (18). As the majority of the new 144 CR sequences were identified as haplogroup D, a Median Joining Network (using Network v4.6; 24) was also generated for just the D haplogroup. DNAsp was used to generate the input file for the Network program. As DNAsp does not allow ambiguous bases and as these ambiguous bases were assumed to reflect sequencing errors, each ambiguous base was modified to reflect the more common of the possible bases within the haplogroup. Default weights were used in Network. To examine the discrepancies between the composition and phylogeographic distribution of haplogroups reported by Storey et al. (1, 25, 26) and those generated in this study, we tested the likelihood of detecting the reported proportions. Tests of statistical significance were performed using the binom.test command and probability distribution graphs were created using the dbinom command (Fig. S9), in the R ‘stats’ package (27). A linear regression plot (Fig. S10) was also generated to visualize the correlation between occurrence of the characteristic 4 CR SNPs of the Polynesian chicken and longitude using the standard plotting function in R (27). Population genetic and differentiation statistics were estimated in Arlequin v3.5 (28) for each population. Bayesian Serial Simcoal (BayeSSC) simulations Bayesian coalescent simulations (using Bayesian Serial Simcoal – BayeSSC v1.0; 29) were used to model eight possible scenarios of chicken colonization of the New World via either 1) Polynesia or 2) Europe. Low level migration between populations was 1) permitted or 2) not permitted, and two separate datasets were examined: 1) only containing haplogroup D ancient samples (representing authenticated Polynesian chicken signals); and 2) containing all putative ancient haplotypes (ancient samples from haplogroups B, D, E; 1, 25, 26; this paper ). In order to test between the different migration routes in BayeSSC, we modeled the same uniform priors for modern population deme size and population growth for each of the migration scenarios to maintain similar demographic parameters. All eight of the South American migration simulations were performed using common uniform priors on modern effective population sizes (MSEA: 10,000-2,000,000; ISEA: 10,000-1,000,000; Europe: 10,000-1,000,000; South America: 1,000-1,000,000; and Pacific: 1,000-1,000,000), with the total panmixia model having a uniform prior with a slightly lower minimum and slightly higher maximum (10,000- 10,000,000). The uniform prior on the growth rate since the last migration event (which differs for each model – see Figure S11) was also common across all eight migration scenarios (growth rate of -0.00001, which equates to 0.001% per generation). Although the generation time of free-ranging domestic chickens is not known, we have estimated a generation time of a year. We considered this appropriate as we were attempting to model early historic chicken populations, which would have had relatively short life spans and low fecundity due to their value as a food source of both meat and eggs. The samples included in the BayeSSC simulations and the migration matrices used are provided in Tables S2-3 and S4-5, respectively. To explore likely demographic histories for chickens in western Polynesia, we also used BayeSSC to simulate alternate migration route hypotheses for comparison with the observed phylogeographic patterns within the Pacific. Sequences in the 201bp CR dataset from the Pacific and ISEA that had location details (n=177) were used to model five possible scenarios of migration routes through western Polynesia, Micronesia and eastern Polynesia (see Fig. S12): a total panmixia model; two models that describe the colonization of Micronesia but with no onward link to Polynesia (one from the Philippines-Micronesia [P-M; arrow 2A in Fig. 1] and the other from New Guinea-Micronesia [NG-M; arrow 3]); and two models that describe Micronesia as a stopping point in an onward route to Polynesia (one from the Philippines- Micronesia-West Polynesia [P-M-WP; arrows 2A and 2B] and the other from New Guinea-Micronesia-West Polynesia [NG-M-WP; arrows 3 and 2B]). Note that the alternate scenario of migration from Micronesia to New Guinea was not tested. The Pacific migration scenarios also had common uniform priors on modern effective population sizes (Philippines: 10,000-2,000,000; PNG: 10,000-2,000,000; Micronesia: 1,000- 1,000,000; Melanesia: 1,000-1,000,000; Western Polynesia: 1,000-1,000,000; Eastern Polynesia: 1,000- 1,000,000), and a common uniform prior on growth rate since the last migration event at 750 BP (growth rate of -0.00001, which equates to 0.001% per generation). The samples included in the Pacific BayeSSC simulations and the migration matrices used are provided in Tables S7 and S8, respectively Supplementary Information Ancient Pacific sample (PAQANA011) Repeated amplifications and Sanger sequencing of Storey et al.’s PAQANA011 sample (1) placed it within the D haplogroup (Dataset S3), however it also highlighted 10 type 2 transitions (C-to-T or G-to-A) across the 12 amplicons. This type of transition is commonly observed in aDNA because of post-mortem template damage, with the hydrolytic loss of amino-groups from cytosine converting the base to uracil, which DNA polymerases read as a thymine base (30). As these internal PCR replications confirmed the discrepancy between our extraction (ACAD9068) and Storey et al.’s (1) published sequence (EF535246) for