Abstract
The inhabitants of New Guinea and its outlying islands have played an important role in the human history of the Pacific region. Nevertheless, the genetic diversity, particularly of pre-colonial communities, is still understudied. Here we present the ancient genomes of 42 individuals from Papua New Guinea (PNG). The ancient genomic results of individuals from Watom Island (Bismarck Archipelago) and the south and northeastern coasts of PNG are contextualized with new (bio-) archaeological data. The individuals’ accelerator mass spectrometry (AMS) dates span 2,500 years of human habitation, and our results demonstrate the influences of different dispersal events on the genetic make-up of ancient PNG communities. The oldest individuals show an unadmixed Papuan-related genetic signature, whereas individuals dating from 2,100 years before present carry varying degrees of an East-Asian-related contribution. These results and the inferred admixture dates suggest a centuries-long delay in genetic mixture with local communities after the arrival of populations with Asian ancestry. Two geographically close communities on the South Coast, AMS dated to within the past 540 years, diverge in their genetic profiles, suggesting differences in their interaction spheres involving groups with distinct ancestries. The inferred split time of these communities around 650 years before present coincides with intensified settlement activity and the emergence of regional trade networks.
Similar content being viewed by others
Main
Humans occupied Near Oceania from at least 50,000 years before present (BP)1,2,3,4,5. The region comprises the Aru Islands, New Guinea, the Bismarck Archipelago and the Solomon Islands. Despite the critical role that its occupants have played in the genetic history of the adjacent regions of Indonesia6,7, northwestern Remote Oceania8 and western Remote Oceania9,10,11,12, its past genetic diversity remains unstudied through ancient genomes.
Apart from being the setting of some of the earliest maritime dispersals by modern humans, Near Oceania was the embarkation point for human dispersals into Remote Oceania, culminating with the settlement of some of the last islands on earth to be permanently inhabited. As part of a dispersal throughout island Southeast Asia, starting around 5,000 BP13,14,15, people, presumably Austronesian-speaking, arrived in the Bismarck Archipelago by 3,300–3,200 BP16,17. Recognized in the archaeological record by distinct ornately decorated pottery and a lifestyle that incorporated horticulture, domesticated animals, long-distance seafaring and maritime and terrestrial foraging18,19, this combination of traits is known as the Lapita cultural complex20. The linguistic diversity of modern Near Oceania includes closely related languages of the Austronesian family and various unrelated non-Austronesian (that is, Papuan) languages, thought to have derived from languages spoken by the regionʼs original inhabitants21. Archaeological evidence suggests this diversity is the result of high mobility, a deeply complex history of interactions before the arrival of Austronesian-speaking peoples, their repeated settlement of the coasts22,23,24 and interactions with local inhabitants. However, the nature and full extent of the interactions remain unknown25.
Until recently, occupation sites of Lapita-associated populations in Near Oceania seemed primarily restricted to smaller offshore islands near the larger island arcs26. Archaeological excavations have challenged this view for the middle and late Lapita period on mainland Papua New Guinea (PNG) with the discovery of sites featuring Lapita pottery and associated settlements in well-dated contexts on the south coast of PNG27 and the Massim region off of southeastern New Guinea28. A small number of Lapita sites are documented on the north coast of New Guinea and the Vitiaz Strait. The Bismarck Archipelago is home to a small number of Lapita sites, suggesting similar timing and process of settlement to the south coast, and implying that Austronesian languages appeared or proliferated in these regions during the middle to late Lapita period29,30,31,32,33.
From 2,200 BP onwards, intensive settlements on the south coast of PNG are associated with shell-impressed pottery and tools manufactured in an exchange sphere extending across ~700 km (refs. 25,34), providing evidence for an occupation by groups culturally and perhaps genetically descended from the Lapita Cultural Complex27,35,36. The restricted present distribution of Austronesian-speaking people along the south coast37 supports this idea, as do the local oral traditions38,39. In the period between 1,200 and 500 BP, a major change in settlement strategies is observed, and long-distance connections between coastal communities were disrupted. During this so called ‘Papuan Hiccup’40 or ‘Ceramic Hiccup’41, settlement along the eastern south coast continued, but locations further west, previously inhabited for centuries, were abandoned, suggesting the re-organization and relocation of populations locally and regionally42. Pottery with increasingly localized forms and decorative motifs emerged within this timeframe40,41,42,43. After 700 BP, many new settlements were established, and several former locals were reoccupied. However, it remains unclear whether this is the return of descendants of the previous occupants or the arrival of new populations44.
This study provides a genetic perspective on the diversity of ancient inhabitants of different regions in PNG and the Bismarck Archipelago. We investigate critical questions regarding the different ancestries involved in the formation of the genetic make-up in the region, specifically the genetic impact of the dispersal of Austronesian-speaking groups. We provide further insights on different human dispersal events or migration processes in the region, such as the settlement of Vanuatu and the Mariana Islands. The genetic data are contextualized with partly new (bio-) archaeological data in the form of isotopic information and microparticles from dental calculus to place these individuals into a cultural and environmental context in PNG by providing evidence for patterns of human diet and mobility.
Results
AMS dating, corrections, isotopic and microparticle evidence
To understand possible chronological changes in ancestry, we produced new accelerator mass spectrometry (AMS) radiocarbon dates for 28 individuals and included previously published dates for four individuals (Supplementary Tables 1 and 2). Stable isotope data were used to correct for the marine reservoir effect45,46,47,48,49, observed as a result of marine based diets, when necessary50,51,52,53,54.
The radiocarbon determinations (Supplementary Table 2) date the individuals at Nebira 2 on the southern coast from 540–320 cal BP (calibrated years before present) up to 420–0 cal BP. Largely overlapping in time, individuals from the nearby site of Eriama date to around 470–310 cal BP to sometime within the last 280 cal BP. An individual from the site of Tilu, on the northeastern coast, dates to around 690–500 cal BP. The five individuals excavated from Watom show wide time intervals, covering a period from 2,690–2,110 cal BP to 630–510 cal BP (Table 1 and Fig. 1b). From the stable isotope data, Eriama and Nebira 2 were found to be consuming a fully terrestrial diet. Tilu and Watom individuals consumed a partially marine diet, consistent with zooarchaeological analysis55,56, and required AMS dating corrections (Supplementary Tables 2 and 6).
a, Map of Near Oceania showing the ___location of the sites discussed in this study, the number of individuals analysed per site and other places mentioned in the text. b, Date ranges for each site/individual in calibrated years BP. Date ranges are based on directly dated skeletal remains and do not necessarily represent the entire occupation of the site. The colours indicate the different archaeological sites or complexes and are used to represent the individuals from these sites in other figures. The asterisk indicates that the date is based on layer information70.
Additionally, the isotope data provided information about past mobility and diet at the sites (Supplementary Table 6 and Extended Data Fig. 1). For Watom in the Bismarck Archipelago, human strontium isotope (87Sr/86Sr; Supplementary Table 4 and Extended Data Fig. 1c) data were consistent with a local signature for all individuals, suggesting an upbringing on the island or one with a similar underlying geology. On the south coast of PNG, wallaby (n = 13) tooth enamel 87Sr/86Sr ratios provide a local baseline, helping to interpret local and non-local individuals at the Nebira 2 site (Supplementary Tables 5 and 6 and Extended Data Fig. 1d). The majority of human 87Sr/86Sr ratios analysed from the Nebira 2 assemblage were consistent with the local signature, suggesting the childhood of these individuals was spent at the site or nearby. Five individuals displayed a non-local 87Sr/86Sr ratio, suggesting they spent their early years away from the site, possibly nearer to the coast, before moving to the site sometime before their death54. The results of a Ward’s hierarchical cluster analysis of the δ13Ccollagen, δ15Ncollagen, δ13Cdentine, δ15Ndentine, δ13Ccarbonate and 87Sr/86Sr of the humans from the Nebira site (Extended Data Fig. 1e) show that the non-local individuals grouped together, indicating their child diet, adult diet and childhood residency were similar to one another but different than the ‘local’ Nebira individuals. The tooth δ13Ccarbonate and δ13Cdentine values of these five non-local individuals were higher than the ‘local’ individuals from the site, indicating that they may have been eating more marine foods at a young age (Supplementary Table 6b). The δ13Cdentine, δ15Ndentine and δ13Cenamel values for the wallabies provided a reference for herbivores eating C3 and C4 plants (Extended Data Fig. 1f). The human dietary isotope results suggest that both the Nebira and Watom communities were eating a diet high in terrestrial resources such as horticultural foods and wild and domestic animals that included wallaby at the Nebira site on the south coast (Supplementary Table 2 and Extended Data Fig. 1f).
Additionally, microparticle analysis recovered from the dental calculus at Watom showed that arboriculture was an important part of the diet. The majority of phytoliths recovered originated from palm or other trees (Supplementary Table 7). Some of the tree phytoliths may also indicate medicinal use as they are likely to come from the bark or leaves of trees57.
Ancient genetic variation
To explore the genetic variation in ancient Near Oceania and archaeological evidence of mobility in the region, we generated genome-wide data for 41 ancient individuals from the sites on the PNG mainland and the Bismarck Archipelago (Supplementary Table 1 and Supplementary Information). To overcome the poor ancient DNA (aDNA) preservation, we used a targeted capture approach, enriching for 1.2 million single nucleotide polymorphisms (SNPs) across the genome. Contamination estimates (Supplementary Table 9) were low with an average of 2% contamination.
Genetic variability between present day, published ancient and newly produced ancient genomes was investigated through a principal component analysis (PCA) (Supplementary Table 10 and Fig. 2) and an ADMIXTURE analysis (Extended Data Fig. 2). Ancient individuals from the Papuan coast included in this study do not cluster with present-day individuals from the PNG highlands but instead with present-day populations inhabiting the southern and northern coastlines of PNG and the Massim Archipelago. Illustrating a cline, extending from the present-day individuals from the PNG highlands to present-day East Asian and ancient Early Remote Oceanians (ERO), the individuals from Eriama on the southern coast cluster on the Papuan end of the cline together with the two individuals from the site of Tilu on the northern coast. However, two individuals (ERI004/ACV-4 and ERI006/ACV-6) are removed from the group and shifted towards East Asians. Individuals from the nearby site of Nebira 2 cluster close together around the midpoint of the same cline. The observations from the PCA are supported by a genetic grouping analysis using f4-statistics and qpWave (Supplementary Table 11), confirming a genetically homogeneous community at Nebira 2, whereas individuals from Eriama form two groups. The placement of the Eriama and Nebira 2 individuals suggests a genetic composition derived from a mixture of Papuan-related ancestry most similar to present-day highland populations of PNG (Extended Data Fig. 3) and of East-Asian-related ancestry. The three older individuals excavated from Watom were analysed separately (WAT002/B10, WAT005/B14, WAT006/B15), considering low coverage in WAT006/B15. WAT005/B14 and WAT006/B15 were consistent with exclusively Papuan-related ancestry (Supplementary Table 12) and have a considerable time interval of >600 years to WAT002/B10 (1900BP and WAT005/B14 ~ 2600BP), which also showed a notable shift towards Asian populations in the PCA (Figs. 1b and 2 and Supplementary Tables 1 and 2). WAT001/B1 and WAT003/B12 were grouped based on more recent dates and genetic similarity.
PCA of present-day individuals (upwards triangles) from Asia, Island Southeast Asia (yellow colours, left side of the plot), Near Oceania (purple colours, right side of the plot) and Remote Oceania (green and turquoise colours) with ancient individuals (circles) projected. Outlined individuals are newly reported in this study. Individuals with insufficient number of SNPs ( < 20 000) are shown transparent (WAT006, NBR008 and TIL004).
Ancestry modelling
Having established that, except for WAT005/B14 and WAT006/B15, none of the individuals derived from only one stream of ancestry using qpWave (Supplementary Table 12), we modelled the East-Asian-related ancestry using present-day Indigenous Taiwanese (Amis) and the Papuan-related ancestry using New Guinea Highlanders. qpAdm tests (Supplementary Table 12) show that 34 of 36 individuals can indeed be modelled as a mixture of these two ancestries (Fig. 3a). However, f4-statistics add detail, showing ancient individuals from the mainland of Papua New Guinea and the Bismarck Archipelago differ in their affinities to present-day Near Oceanic populations. Whereas the individuals from Eriama, Nebira 2 and Tilu show higher affinity to highland New Guinean populations, all individuals from Watom show higher affinities to Baining-speakers (a non-Austronesian language) from New Britain (Supplementary Table 11 and Extended Data Fig. 3a). The Nebira 2 and Eriama individuals with higher East-Asian-related ancestry proportions show more affinity to the ERO from Vanuatu and Tonga compared to present-day Indigenous Amis from Taiwan, ~2,650-year-old individuals from Guam or ~4,500-year-old individuals from the site of Suogang on Penghu Island, west of Taiwan (Supplementary Table 11 and Extended Data Fig. 3b).
Ancestry modelling with qpWave using Ami (yellow) and New Guinea (purple). a,b, Representing the ancestral components for all individuals from Papua New Guinea (a) and for Watom Island (b). c, Ancestry modelling for WAT002 using Lapita-related individuals from Vanuatu and Tonga (green) and an individual from Vanuatu with Papuan ancestry closely related to present-day Baining Marabu from the Bismarck Archipelago (TAN002, teal green). White lines indicate the standard error for each component, p values are indicated on the right-hand side of each bar, values above 0.05 are considered working models, non-working models are shown in grey font. d, Inferred dates of admixture (entire range indicated by filled blocks) and the midpoint of the calibrated AMS date range, based on the 95% probability ranges, for each individual (circles). Grey lines indicate the full AMS range.
The patterns observed in the PCA and f4-statistics are supported by the qpAdm analysis. East-Asian-related ancestry in individuals from Nebira 2 ranges from 45–60% and is higher compared to Eriama and Tilu, where the proportion is around 20%. The two individuals from Eriama (ERI004/ACV-4 and ERI006/ACV-6) not included in the main group are the exception, showing higher East-Asian-related ancestry proportions of around 35%.
The individual WAT002/B10, dated to around 1,900 BP, shows admixture between Papuan-related and East-Asian-related ancestry, at 40% to 60%, respectively. However, the East-Asian-related proportion is reduced to 20% when ancient individuals are used as a proxy for this ancestry (Fig. 3c).
Timing of the admixture events
To estimate the time of admixture between the two ancestry components present in the individuals, we performed an admixture dating analysis using DATES (Fig. 3d and Extended Data Fig. 4). For the ~1,900-year-old individual from Watom Island (WAT002/B10), we inferred an admixture around ten generations before, resulting in a date of ~2.100 BP for the admixture event and for the two individuals dating to ~500 BP from the same site.
We again observe slight differences between the two sites on the southern PNG coast. Inferred from the dated individuals only, individuals from Eriama have an average admixture date of 1,100 BP (ranging from 1,600 to 600 BP), whereas the average date for Nebira 2 is slightly older (around 1,500 BP, ranging from 1,800 to 900 BP) (Supplementary Table 13). ERI007/ACV-7 and NBR020/ACJ-34 are the exception with younger admixture dates of 370 BP and 650 BP, respectively. Coincidently, NBR020/ACJ-34 displayed a non-local 87Sr/86Sr ratio and childhood dietary values (Extended Data Fig. 1d). Because of low coverage, individuals from Tilu were grouped for this analysis, resulting in an inferred admixture event of 900 BP.
Sex-biased mixture
We determined the Y-chromosomal haplogroups for 21 of the 25 genetically male individuals and mitochondrial DNA (mtDNA) haplogroups for 38 of the 41 individuals in the dataset (Table 1 and Extended Data Fig. 5).
The majority of individuals carry mtDNA haplogroups associated with present-day and ancient East Asian populations58,59,60, dominated by the mtDNA haplogroup B4a1a1. Only four individuals carry Papuan-related mtDNA haplogroups. Conversely, the Y-chromosome haplogroups indicate the opposite pattern. Apart from one Y-chromosome haplogroup with East Asian origin (O2a2b2), all other 18 identified haplogroups are most common in Near Oceania (Table 1 and Extended Data Fig. 5). The pattern of sex-biased mixture is also detectable at a genome-wide level, shown by the analysis of mixture proportions on the X chromosome, comparing to those inferred from the autosomes. The analysis shows an excess of Austronesian-related ancestry on the X chromosome ranging from 10 to 60 percentage points, suggesting a sex-biased mixture as previously observed in similar to ancient individuals from Vanuatu and Wallacea7,9 (Supplementary Table 12).
Settlement sequence of the broader region
To investigate the relationships with other ancient groups within the region, we modelled an admixture graph using qpGraph. The source region for the initial settlement of the Mariana Islands is under debate. Radiocarbon dates indicate settlement ~ 3,200 BP61, but some archaeological evidence suggests even earlier occupation several centuries earlier62. Potential source regions include the Philippines63,64 and, based on highest probabilities of landfall in the Marianas during a sea voyage, northern Near Oceania65. In the absence of ancient genomes of reasonable quality from the Philippines, we evaluate this indirectly. Using ancient genomes from Guam and Saipan, dating to between 2,800 and 200 BP8,66, we tested two competing hypotheses. First, the Marianas could have been settled from the Philippines, in which case the group is modelled to have split before the formation of the Lapita-related genetic component, including the East-Asian-related component in WAT002/B10 (Extended Data Fig. 6b). Alternatively, the Marianas could have been settled from Near Oceania, in which case their group should split after they arrived in the Bismarck Archipelago (Extended Data Fig. 6c). The data provide a better fit for the model of a split of the genomes from Guam and Saipan before WAT002/B10 (|Z| = 2.883 vs |Z| = 4.196 for the alternative settlement from Near or Remote Oceania), suggesting the origin of the first settlers of the Mariana Islands lies in the Philippines.
Genetic relationships and burial patterns
Finally, we investigated the genetic relationships between the individuals for Nebira 2, the only site where familial relationships could be addressed (Extended Data Fig. 7 and Supplementary Table 14). The Nebira 2 cemetery contained primary burials, some of which were used as multiple or simultaneous burials67. People were interred in a supine extended position, mostly in western or northwestern orientation. The resulting relatedness network shows that the first- and second-degree relations at the site include both genetically male and female individuals, showing no clear signs of patri- or matrilocality. Genetically related individuals are not interred together but are buried close by, whereas unrelated individuals can be found in multiple or simultaneous burials. These results suggest that the reason for multiple burials or reusing graves was not necessarily family related or that perhaps family was not restricted to genetic relatedness as is often the case in PNG social contexts, that is, local Motu and Koita groups who inhabit the region today68.
Analysis of runs of homozygosity (RoH) can potentially reveal higher levels of consanguinity and approximate the effective population size, which in a simplified scenario reflects the number of individuals reproducing, not the total number of individuals. The RoH show that the population at Nebira had lower levels of consanguinity and an estimated effective population size of 800–3,200 (Extended Data Fig. 8). For Eriama, a large number of short RoH is suggestive of higher background relatedness and a smaller effective population size of ~400–1,600.
The analysis of segments in the genome that are identical by descent (IBD) reveals distant genetic relations between various individuals from Eriama and Nebira (Supplementary Table 14, Extended Data Fig. 7 and Supplementary Fig. 7). The number and size of segments shared between multiple individuals from the two sites and a two-population split model suggests the two groups ceased contact around 13 generations prior, resulting in an estimated split of ~640 years ago (Supplementary Information).
Discussion
Time transect in Watom
The Reber–Rakival (SAC) site on Watom Island has a well-established sequence, providing evidence for occupation by people of the Lapita cultural complex from 2,800 to 2,350 BP, the middle and late Lapita phases in the Bismarck Archipelago46,69. The lowermost layer (Layer D) presents evidence for an earlier occupation70. Ceramics and obsidian are absent in this earliest layer, possibly pointing to a different population occupying the island before the arrival of the Lapita cultural complex or, potentially, co-occupation with Lapita communities during the middle and late Lapita periods on the island. At other sites on the island (SAD, SDI, SAB) both dentate-stamped Lapita pottery and sherds with opposed pinching nail impression decoration are found in secure contexts that date to before 3,000 BP (middle Lapita period).
The first analysis of ancient genomes from the Bismarck Archipelago show diverse ancestries on Watom, changing through time. WAT005/B14 and WAT006/B15 show an occupation by people with only Papuan (Baining-like)-related ancestry during the middle to late Lapita period (about 2,700–2,200 BP), at least 600 years after the archaeologically attested arrival of the Lapita cultural complex in the Bismarck Archipelago about 3,300–3200 BP (Fig. 2a and Extended Data Fig. 3a). The 87Sr/86Sr ratio of WAT006/Burial 15 indicated this individual was probably local to the island, as were the other individuals recovered from the site (Extended Data Fig. 1e). The archaeological and anthropological analysis of WAT006/B15 provides clues that support the different ancestry compared to the later individuals. WAT006/B15 exhibited a rare case of cultural cranial deformation to elongate the cranium (details in the Supplementary Information and Supplementary Fig. 2)—a practice not currently known from the Lapita Cultural Complex but observed in the Aware region of southwest New Britain during the historic period71,72. WAT006/B15 was also the only definitely seated burial, with their head fallen into their lap. Although WAT005/B14 was disturbed with their lower half extending into the baulk, it could be determined that the individual was interred with their upper body in a supine extended position and their head oriented west to northwest. Supine extended and supine with limbs semi-flexed was, with a few exceptions, the general pattern of the later burials at the site (Petchey et al.70).
By 1,900 BP, we find a genetically female individual (WAT002/B10) who fits the model of mixture between Papuan-related and East-Asian-related ancestry, similar in the proportions to the individuals from the late Lapita period in Vanuatu (Extended Data Fig. 3b). Additionally, WAT002/B10 shows higher affinity to present-day non-Austronesian-speaking Baining from New Britain, similar to contemporary individuals from Vanuatu (Extended Data Fig. 3a). The connection to Vanuatu is mirrored in the incised and applied relief pottery from Watom, which shares similarities with ‘Mangaasi’-style pottery from Vanuatu73 and was found alongside dentate Lapita pottery at the Watom SDI ___location in a layer dating to 2,210–1,755 BP. At the Watom SAB site, no Lapita pottery but two types of locally made incised and applied relief style pottery were found in stratigraphic contexts dating to 2,455–2,310 cal BP. There may be some connection between the later arrival of Lapita individuals to Watom and incised and applied relief pottery in the late Lapita sequences on the island. Finally, the two most recent individuals from Watom (about 670–510 BP) carry a genetic make-up very similar to Austronesian-speaking Tolai people living in the region today (Fig. 2). However, they are not consistent with entirely deriving from that population (Supplementary Table 12). This is an expected observation when considering the island of Watom was highly affected by natural disasters74,75 and colonialism that incited population mobility76. Evidence for population movements within the last ~500 years also comes from the oral traditions of the Tolai, which records that they came from southern New Ireland77. This origin is supported by the belief system and other cultural practises the Tolai share with groups in southern New Ireland but not with other groups in the Willaumez Peninsula72.
Proxy information regarding the settlement of the Mariana Islands
The branching pattern of the East-Asian-related ancestry of the admixed individual from Watom in a tree-like model (WAT002/B10; Extended Data Fig. 6) favours a settlement of the Mariana Islands from island Southeast Asia. The favoured route requires sailing against prevailing winds and currents and excludes a less challenging passage from the islands northeast of PNG. The lack of adequate ancient samples from the Philippines and Taiwan leaves some uncertainty to the exact geographic ___location from where the initial settlement of the Mariana Islands started, but regardless, expert navigation skills have to be assumed for the populations settling the islands in northwestern Remote Oceania.
Population history of coastal PNG
The genetic affinities of the inhabitants of the southern coast of PNG show that the occurrence of pottery associated with the Lapita Cultural Complex was probably linked to the arrival of populations with East-Asian-related ancestry. All individuals from the coast of PNG harbour East Asian ancestry not observed in present-day individuals from the highlands of New Guinea but present in coastal groups in PNG today (Fig. 2a and Extended Data Fig. 1)78. Whereas both the Eriama and Nebira 2 sites on the southern coast are geographically close, their genetic composition suggests different population histories. On average, the individuals from Eriama show higher Papuan-related ancestry of ~80%, albeit with higher variation and younger admixture dates ranging from 1,600 to 600 BP (Fig. 3d) compared to the individuals from Nebira 2. The latter are more homogeneous, with higher genetic variation, higher proportions of East-Asian-related ancestry of > 50% and, on average, older admixture dates ranging from 1,800 to 900 BP (Fig. 3a,d,e).
The two individuals from the site of Tilu, situated on the northeastern coast of New Guinea, show a genetic profile and admixture dates similar to Eriama (Fig. 3a,d and Supplementary Table 13). First evidence for sustained settlement of the site securely dates to 650 BP, with possible earlier traces of habitation 900–800 BP79, by people possibly speaking Austronesian languages of the Ngero–Vitiaz network37,80. The admixture event inferred from the individuals from Tilu falls around the time of first occupation, suggesting admixture occurred during initial settlement of the site or establishment by an already admixed population. Archaeological research places the site of Tilu in a local trading network extending to the Willaumez Peninsula on the northwestern coast of New Britain in the Bismarck Archipelago81. Linguistic evidence and oral traditions suggest an origin in the Vitiaz Strait in the Bismarck Sea for the Ngero–Vitiaz languages of the Bel subgroup spoken today in the area of the Tilu site, Gedaged and Bilbil80,82,83. However, the material exchange seems not to have extended to genetic exchange, as there is no higher affinity of the individuals of Tilu to the populations of the Bismarck Archipelago, compared to Eriama (Supplementary Table 11 and Extended Data Fig. 3a).
Timing and genetic impact of regional dispersals
The date of ~2,100 BP for admixture between the two ancestries for all individuals from Watom with East-Asian-related ancestry suggests that admixture with local Papuan people either occurred repeatedly starting with the arrival of the Lapita Cultural Complex or only started ~1,000 years after the arrival of the Lapita cultural complex in the Bismarck Archipelago (Fig. 3d). We acknowledge the limitations of interpretations deriving from only five individuals covering a wide time transect. As previously observed7,9, repeated admixture events may result in seemingly recent admixture dates. However, we would cautiously interpret the results as further support that the early settlers in Remote Oceania largely derived from a Lapita-associated population of East Asian ancestry that did not mix with local groups in Near Oceania. A local admixture upon or shortly after arrival in the Bismarck Archipelago would have resulted in a much earlier admixture date of ~3,000 BP. Additionally, the presence of individuals with entirely Papuan-related ancestry in the late Lapita period on Watom is consistent with the observation that the oldest individuals with exclusively Papuan-related ancestry reached Vanuatu ~2,600 BP, around 300–400 years after the initial settlement of the archipelagos9,10,84.
The admixture dates inferred from the dated individuals on the south coast show this admixture event for people in Nebira 2 occurred at around 1,100–1,500 BP, much later than in other places such as Vanuatu (~2,600 BP) and Watom (~2,100 BP). This could be interpreted as the result of a later arrival of the descendants of people associated with the Lapita Cultural Complex along the southern coast of PNG. However, evidence for the first occurrence of Lapita settlements in the south coast of PNG is dated to starting from ~2,900–2,800 BP27,28, 400 years after the initial Lapita formation in the Bismarck Archipelago16,26 and contemporaneous with archaeological evidence from Vanuatu85. From our analysis, it is not clear whether this late admixture date is a result of long genetic isolation of the first settlers related to Lapita Cultural Complex or the result of repeated admixture with local populations carrying varying proportions of Lapita-related ancestry. Several sites demonstrating cultural contact between the Lapita Cultural groups and local groups show that interactions started shortly after the appearance of the Lapita-associated settlers28, rendering the isolation hypothesis less likely.
The archaeological record of the south coast of PNG shows ‘pulses’ of settlement intensification and ceramic abundance86. This period, coined the ‘Ceramic Hiccup’41 but also known as ‘Papuan Hiccup’40 was possibly the result of an abandonment or relocation of sites, with a disruption in ceramic production86,87 and styles40,88, or the contraction of interactions to the eastern parts of the southern PNG coast42,89, where networks were maintained, although perhaps less intensive41. This period is believed to have lasted between ~1,200 and 500 BP, however the exact timing remains uncertain as a result of insufficient radiocarbon dating42. The assumed period coincides with the medieval climate anomaly (1,250–700 BP), which was possibly delayed in Oceania, however a lack of proxy sites from PNG and Western Oceania leaves the magnitude of the effect open for the specific region90. Still, an increase in El Niño/Southern Oscillation events91 may have led to seasonal droughts, suggesting the changes in settlement patterns were at least partly influenced by the availability of fresh water and the viability of cultivation22,92. Despite whether climate had a direct effect on the coastal communities, the interruption of long-distance trade links might have led to abandonment of many long-term settlement locales, relocation further inland40 and perhaps been a catalyst for conflict between groups. Certainly, settlement of Nebira 4, at the base of the Nebira hill, was abandoned by the time settlement of Nebira 2, on the saddle of the double-peaked hill commenced93,94,95. A shift to defensive hilltop settlement suggests that the relationships between groups might have changed, and defensive sites on hilltops became the preferred occupational area42.
The majority of individuals from Nebira 2 show a mixture event between 1,600 and 1,100 BP, overlapping with this period. The individual NBR020/ACJ-34 shows an admixture date of only ~650 years, much younger than the majority of the group. Additionally,this individual displayed a non-local 87Sr/86Sr ratio54 and a different childhood diet, suggesting they may have originated from a population with a different ancestral history. The admixture date of NBR020/ACJ-34 coincides with a period when coastal and inland settlements had become archaeologically visible, and ceramic wares show stylistic innovations with regional variation indicating the establishment or reinforcement of social boundaries96,97,98,99. Trade supposedly resumed, perhaps with reconfigured networks, in ~800–500 BP and has been associated with different Austronesian-speaking groups arriving on the coasts78.
Overlap of maritime trade networks such as the Kula and hiri trade networks by at least 500 BP42,86,100,101, suggests a (re-)emergence of these networks in the period following the ‘Papuan Hiccup’42.
Admixture dates inferred from the individuals from Eriama (Fig. 3d), together with the higher Papuan ancestry (Fig. 3a and Extended Data Fig. 3) indicate a continuation of genetic exchange beyond that identified for Nebira 2. The strong shift from 40% to 15% Asian-related ancestry suggests they were part of an interaction sphere with people carrying higher Papuan-related ancestry in the model represented by Papuan Highlanders. However, we lack the resolution to identify which populations contributed and where their geographical source was. During the ‘Papuan Hiccup’, communities might have retreated into the interior areas, where interactions with highland populations might have intensified40. One individual from Eriama (ERI007/ACV-7) shows more recent admixture inferred to around 600–400 years ago, also implying that the population in Eriama continued genetic exchange with other populations.
Despite the many differences between the two sites, analyses of IBD blocks reveal distant connections between Nebira 2 and Eriama, suggesting that they constituted one reproductive unit ~13 generations before the analysed population, estimating the split to ~640 BP. The different genetic compositions of the two nearby sites show that the southern coast was a genetic, and possibly also a cultural and linguistic mosaic of people, matching the situation in Motu and Koita speaker communities today. The Motu language is a Central Papuan Tip language of the Western Oceanic branch of the Oceanic group of Austronesian languages102,103 and is spoken mostly by people living in the coastal region of the Central Province. The Papuan language dominant in the region, Koita104, is spoken in settlements more inland38,39,105,106. The two major cultural groups occupying the coastal and interior areas in the vicinity of Nebira 2 and Eriama have historically close social connections that regularly included intermarriage, with both groups involved in the maritime hiri trade expeditions107. Differences in burial practices, specifically primary burial at Nebira 2 and secondary cave burial at Eriama, also point to distinct cultural practices regarding the treatment of the dead.
The oral traditions of present-day groups remembering their ancestors residing at Nebira 2 and Eriama include tales of a relocation from the mountains67,105. The reason for the separation of the two communities, despite no evident geographical barriers, remains enigmatic and might point to the introduction of a cultural barrier connected to different interactions spheres of the two groups after relocation.
Methods
Sample processing
All samples were processed in dedicated ancient DNA laboratories at the Max Planck Institute for the Science of Human History (now MPI for Geoanthropology) in Jena, Germany.
Bone powder from the petrous part of the temporal bone was obtained through cutting along the margo superior partis petrosae (crista pyramidis) and drilling 50–150 mg bone powder from the densest part around the cochlea108. Teeth were sampled by cutting along the junction of the root and the crown and drilling ~50 mg from the pulp chamber. In total, 46 samples were destructively sampled of which 41 could be included in the analysis.
DNA extraction was carried out following established protocols109. Negative and positive controls were included. To release DNA from 50–100 mg of bone powder, a solution of 900 μl ethylenediaminetetraacetic acid (EDTA), 75 μl H2O and 25 μl Proteinase K was added. In a rotator, samples were digested for at least 16 h at 37 °C, followed by an additional hour at 56 °C (ref. 110). The suspension was then centrifuged and transferred into a binding buffer as previously described109. Samples were purified over silica columns for high volumes (High Pure Viral Nucleic Acid Large Volume Kit; Roche) with two washing steps using the manufacturer’s wash buffer. DNA was eluted in TET (10 mM Tris, 1 mM EDTA and 0.05% Tween) in two steps for a final volume of 100 μl.
Double-stranded DNA libraries were built from 25 μl of DNA extract in the presence of uracil DNA glycosylase (UDGhalf libraries), following a double-stranded (ds)‘UDG-half’ library preparation111. Negative and positive controls were carried alongside each experiment. Libraries were quantified using the IS7 and IS8 primers112 in a quantification assay using a DyNAmo SYBP Green qPCR Kit (Thermo Fisher Scientific) on the LightCycler 480 (Roche). Each aDNA library was double indexed113 in 1–4 parallel 100 μl reactions using PfuTurbo DNA Polymerase (Agilent). The indexed products for each library were pooled, purified over MinElute columns (Qiagen), eluted in 50 μl TET and again quantified using the IS5 and IS6 primers112 with the quantification method described above. Four μl of the purified product were amplified in multiple 100 μl reactions using Herculase II Fusion DNA Polymerase (Agilent) following the manufacturer’s specifications with 0.3 μM of the IS5/IS6 primers. After another MinElute purification, the product was quantified using the Agilent 2100 Bioanalyzer DNA 1000 chip. An equimolar pool of all libraries was then prepared for shotgun sequencing on Illumina platforms in 75 base pair single-end-run cycles using the manufacturer’s protocol. For nine individuals with low DNA content in ds libraries, we produced a second, single-stranded (ss) library114 with an automated protocol as detailed in ref. 115.
Libraries were further amplified with IS5/IS6 primers to reach a concentration of 200–400 ng μl−1 as measured on a NanoDrop spectrophotometer (Thermo Fisher Scientific). Mitochondrial DNA capture116 was performed on screened libraries which, after shotgun sequencing, showed the presence of aDNA, highlighted by the typical CtoT and GtoA substitution pattern towards 5′ and 3′ molecule ends, respectively. Furthermore, samples with a percentage of human DNA in shotgun data 0.1% or greater were enriched for a set of 1,237,207 targeted SNPs across the human genome117. The enriched DNA product was sequenced on an Illumina HiSeq 4000 instrument with 75 single-end-run cycles using the manufacturer’s protocol. The output was de-multiplexed using bcl2fastq version 2.17.1.14 (Illumina conversion Software) and dnaclust version 3.0.0118.
Pre-processing of the sequenced reads was performed using EAGER version 1.92.55119. The resulting reads were clipped using Clip&Merge119 and AdapterRemoval version 2120. Clipped sequences were mapped against the human reference genome hg19 using the Burrows–Wheeler Aligner (BWA) version 0.7.12121 disabling seeding (−l 16500, –n 0.01). Duplicates were removed with DeDup version 0.12.2119. A mapping quality filter of 30 was applied using SAMtools version 1.3122. In ds libraries, reads were trimmed for two base pairs. Different sequencing runs and libraries from the same individuals were merged, duplicates removed and sorted again using SAMtools122. Trimmed and untrimmed reads were genotyped separately using pileupCaller version 8.6.5 (https://github.com/stschiff/sequenceTools). We combined the genotypes keeping all transversions from the untrimmed genotypes and transitions only from the trimmed genotypes to eliminate problematic, damage-related transitions on the ends. ss-libraries were genotyped based on the untrimmed reads using the –singleStrandMode. The generated pseudo-haploid calls from both ss and ds libraries were merged using a custom python script, which keeps all identical positions across the two genotypes and the sites covered only in one of the two. For sites covered in both libraries, but with different base calls, the state of the genotype was randomly picked from one of the libraries. The final genotypes of all ancient individuals were merged to a pulldown of the 1,240 K SNPs from the Simons Genome Diversity Project123, a set of individuals from Asia and the Pacific, genotyped on the Human Origins array124, on the Illumina Infinium Multi-Ethnic Global Array100 and ancient Asian and Oceanian individuals8,10,66,124,125,126,127.
The typical features of ancient DNA were inspected with DamageProfiler version 0.3.1 (http://bintray.com/apeltzer/EAGER/DamageProfiler)119 (Supplementary Table 8). Sex determination was performed by comparing the coverage on the targeted X-chromosome SNPs ( ~ 50 K positions within the 1,240 K capture) normalized by the coverage on the targeted autosomal SNPs to the coverage on the Y-chromosome SNPs ( ~ 30 K), again normalized by the coverage on the autosomal SNPs128 (Supplementary Table 1). For male individuals, ANGSD version 0.919129 was run to provide an estimate of nuclear contamination in males. All male samples with at least 100 X-chromosome SNPs covered twice exhibited X-chromosome contamination levels below 7%, hence all reads were retained for further analyses (Supplementary Tables 1 and 4). For both genetically male and female individuals, mtDNA-captured data were used to jointly reconstruct the mtDNA consensus sequence and estimate contamination levels with Schmutzi130 (Supplementary Table 1) and used as reliable predictors for nuclear contamination131,132. The software ADMIXTURE version 1.3.0 (69) was used in unsupervised mode to allow for free genetic clustering with a worldwide set of individuals123,133,134,135 (Extended Data Fig. 1).
AMS and isotope analysis
The 14 new radiocarbon dates for this study and the stable isotope analysis was produced at the Curt–Engelhorn–Zentrum Archäometrie gGmbH in Mannheim, Germany. Details on the protocols can be found in the Supplementary Information. Marine reservoir correction was applied to correct for the offset in 14 C between the atmosphere and the oceans (∆R)45,46,47,48,49. ∆R was applied to the most recent marine calibration curve136 after determining the percentage of marine food in the diet from the bone stable isotope values46,49 (Supplementary Tables 2–4). For individuals from Eriama and Nebira, fully terrestrial calibrations were produced. For Burials 1 and 12 (WAT001/003) from Watom, a marine deltaR of 172 ± 72 for shell/mixed diet animals was applied using the marine20 calibration curve and for the individual from Tilu, a ∆R of −111 ±16 and 28 ± 10 marineC137. For burials 15, 10 and 14 (WAT002, 005 and 006) from Watom, a marine carbon contribution of 28 ± 10% marineC was used in the calibration.
Approximately 0.5–0.8 g of cortical bone was sampled for the stable isotope analysis50,51. For the dentine stable isotope analyses, the distal third of the root of the permanent first molar (~0.2–0.3 g, measured from the cemento–enamel junction, formation time between the ages of 5–9 years)138,139 was sampled with a Dremel drill fitted with a diamond edge saw. The dentine was sampled horizontally rather than along growth increment lines for the purposes of attaining enough sample to analyse140. Secondary dentine was removed from inside pulp cavity of the tooth root with a Dremel drill fitted with a diamond burr. The samples of dentine were then sonicated for five minutes and fully dried.
Bone and tooth samples were cleaned with alum oxide air abrasive equipment (Bego Easyblast). A modified Longin method141 was used to extract collagen from the bone and dentine samples at the University of Otago, Dunedin, New Zealand142,143. All bone and dentine samples were soaked in 0.5 M HCl at 4 °C (changed every other day) until completely demineralized. The demineralized samples were then rinsed in MilliQ H2O until they reached a neutral pH. The samples were gelatinized at 70 °C in a pH 3 solution for 48 h, followed by filtering with 5–8-μm Ezee mesh filters (Elkay Laboratory Products) to remove any reflux-insoluble residues and then ultrafiltered with Millipore Amicon Ultra-4 centrifugal filters (30,000 nominal molecular weight limit (NMWL)) to retain molecules larger than 30 kDa (ref. 144). The purified ‘collagen’ was frozen and then lyophilized for 48 h and subsequently weighed into tin capsules before analysis by elemental analyser–isotope ratio mass spectrometry (EA–IRMS) at either IsoTrace (Dunedin, New Zealand) or Iso-Analytical (Cheshire, UK)50,51. Analytical error, calculated from replicate measurements of samples, was ± 0.1‰ for δ 13C and ± 0.2‰ for δ 15N (1 standard deviation (SD)).
The stable isotope analysis of the petrous bones analysed for AMS dating was conducted at Curt–Engelhorn–Zentrum Archäometrie gGmbH, Mannheim by the same method described above, except for the addition of a step to remove humic acids with NaOH after demineralization. Analyses were conducted in triplicates and included (1) combustion of the sample and determination of C and N per cent using a vario PYRO cube CNSOH elemental analyser (Elementar); (2) determination of isotope ratios using a precision isotope ratio mass spectrometer (Isoprime); (3) data are corrected to USGS 40 and USGS 41a using the internal software (two-point-normalization).
For the enamel carbonate analyses, the enamel surface was cleaned of surface contaminants by abrasion with a Dremel rotary tool with a diamond cutting blade, and any adhering dentine was removed. Two enamel chips weighing between 10–50 mg a piece were sampled for carbon and oxygen stable isotope analysis54 and strontium isotope analyses (below)54. For the carbon and oxygen isotope analyses (note that oxygen results not reported here), the enamel chip was ground in a clean mortar and pestle and transferred into an acid-cleaned glass vial. Following pretreatment protocols145 the samples were first soaked in 2% NaOCl (sodium hypochlorite) for 24 h, rinsed with MilliQ H2O and then soaked in 0.1 M acetic acid for another 24 h to remove organic and secondary diagenetic carbonates. Carbon stable isotope analyses from the carbonate were undertaken at the Isotrace Research Laboratory, Dunedin, New Zealand, using a Thermo Delta Plus Advantage linked to a Gasbench II via a GC PAL autosampler. Delta values were normalized and reported against the international standards Vienna Pee Dee Belemnite. Analytical precision was determined through repeated analysis of laboratory standards IA-R022 (calcium carbonate), NBS-18 (calcite) and NBS-19 (calcite) and duplicate samples. Analytical error, calculated from replicate measurements of samples, was ±0.1 for δ 13C (1 SD, n = 10).
For strontium isotope analysis, ion exchange techniques were used to isolate and purify the Sr fraction from the digested sample matrix at the Centre for Trace Element Analysis, Department of Geology, University of Otago. Enamel chips weighing 10–20 mg were sonicated in MilliQ H2O, dried and transferred into clean perfluoroalkoxy alkanes vials (Savillex) and weighed in the clean lab before digestion in 2 ml of 3 M HNO3 solution at 110 °C overnight. Once fully digested, samples were evaporated for four hours at 110 °C. Strontium was manually separated utilizing a micro- chromatographic exchange column, Eichrom Sr-SPEC resin, and the established method of column chemistry146. Only a single elution was necessary for the human samples which were then evaporated and dissolved in a 2% HNO3 solution for mass spectrometric analysis.
The 87Sr/86Sr values were measured using a Nu Plasma-HR MC-ICP-MS instrument (Nu Instruments Ltd., UK). The 87Sr/86Sr data were normalized using repeated measurement of the NIST-SRM 987 standard (n = 25, average 87Sr/86Sr = 0.710286 ± 0.000013 (1 SD) in very good agreement with the accepted value of 0.71034 ± 0.00013 (1 SD)147. An in-house sample of a giant clam (Tridacna) (ANU) carbonate control (n = 5, average 87Sr/86Sr = 0.70920 ± 0.000006 1 SD) which is consistent with expected seawater 87Sr/86Sr value of 0.7092148. Total procedural blanks for the chemical separation process were 60 ng which is negligible relative to the amount of Sr in the tooth enamel samples. Duplicate analyses were also performed on three samples which ranged between ±0.000013 and 0.000027 (2 SE). All uncertainties are reported at the 2-sigma level, unless stated otherwise.
Population genomic analysis
Principal component analyses were performed using smartpca version 13050149 with a set of populations from East Asia and the Pacific78,100,123,133,134,135 (Fig. 2 and Supplementary Table 10). Ancient individuals8,10,66,124,125,126,127 were projected onto the calculated components using the options ‘lsqproject: YES’, ‘shrinkmode: YES’ and ‘numoutlieriter: 0’. Individuals with less than 20,000 SNPs covered were not projected with the exception of WAT006.
To identify the differences on an individual basis and we used qp3Pop version 5.0 (70) and computed f3-outgroup statistics comparing all individuals to each other with Mbuti.DG serving as an outgroup. We used qpDstat version 5.0133 to run f4-statistics of the form f4(Mbuti, Ami.DG Individual 1 site X; Individual 2 site X) (Supplementary Table 10). The values close to zero for all individuals excavated from Nebira suggested a grouping of the individuals by site was reasonable for certain analyses. For individuals excavated from Eriama, the two individuals ERI004 and ERI006 produced absolute Z-scores greater than 3 with all other individuals from the site and hence were kept as a separate group. Both individuals excavated from the site Tilu and the two younger individuals excavated on Watom (WAT001 and WAT003) were also grouped on this basis, whereas WAT002, WAT005 and WAT006 were kept separate, also accounting for the long-time intervals between them (Table 1). To test which present-day populations represented best the Papuan and East-Asian-related ancestries in the individuals, we computed f4 statistics of the form f4(Mbuti.DG, Test, X, New_Guinea) and f4(Mbuti.DG, Test, X, Ami), respectively, testing in X all other populations from the region (Supplementary Table 11). To understand whether the Asian ancestry component was more similar to the Early Remote Oceanians (ERO) from Vanuatu and Tonga124 compared to ancient Austronesians from Taiwan (Suogang)125, we calculated f4(Mbuti.DG, Test; Suogang, ERO), expecting positive test scores for a higher affinity to ERO (Supplementary Table 11 and Extended Data Fig. 3b). To understand the differential affinities with respect to Near Oceanian populations, disregarding the differences in Asian ancestry, we produced a scatterplot (Extended Data Fig. 3a) based on the f4(Mbuti. DG, Ami, Test, Baining_Marabu/ New_Guinea).
We used qpWave version 410150 to test whether individuals were consistent with deriving from the same group as other individuals from the same site, relative to a set of reference groups (Mbuti, Onge, New_Guinea, Baining_Marabu, Ami, Han, English, Chukchi, Nasioi, Denisova_published.DG). To test whether some of the individuals could be modelled as consisting of a single ancestry component, we modelled the respective individual and Ami to test for exclusively Asian ancestry and New_Guinea Highlanders to test for exclusively Near Oceanian ancestry (Supplementary Table 12) using the same references detailed above, excluding the respective populations used in the test151. After identifying the individuals not consistent with deriving from one respective ancestry, we used qpAdm version 5.0133 to model all groups and the individuals in each group covered by more than 50,000 SNPs as a two-way admixture between New_Guinea Highlanders and Amis (Supplementary Table 12) and the grouped individuals as a mixture of Early Remote Oceanians and TAN002, a previously published individual with exclusively Baining-like, ancestry9. As reference groups, we used Mbuti, Onge, Han, Chukchi, English and Denisova_published.DG. To test for sex-biased admixture, we calculated the excess ancestry on the X chromosome restricting the above analysis on grouped individuals restricted to the X chromosome (Supplementary Table 12). We subtracted the value obtained from all chromosomes from that obtained from the X chromosome alone (Supplementary Table 12).
To evaluate the genetic affinity relative to other selected populations in a tree-like representation, we used qpGraph version 6450, with the parameters outpop: NULL, useallsnps: NO, blgsize: 0.05, forcezmode: YES, lsqmode: YES, diag: .0001, bigiter: 6, hires: YES, lambdascale: 1. We used this tool strictly for testing two hypotheses: (1) the genetic composition of the individuals from Guam and Saipan8,66, in this model grouped into one ‘Mariana’ group, is more similar to that of Austronesian populations in Taiwan or the Philippines, suggesting the Marianas were settled directly from Island South East Asia; (2) the ancestry is more similar to the Asian component of WAT002, suggesting that the expansion first arrived in the Bismarck Archipelago, and from there, facilitated by the prevailing winds and currents, travelled to the Mariana islands. For completeness of possible scenarios, we also tested (3) if the ERO component in Remote Oceania is ancestral to Watom and Guam/Saipan. To test these three models, we first constructed a base tree including Mbuti.DG, Ami.DG, Igorot.DG, Papuan.DG123, then adding according to their chronological age ERO9,124, followed by the ancient genomes from Guam and Saipan8,66. The best-fitting base tree models Ami.DG and Igorot.DG as sister groups shows no admixture events and is in concordance with the genetic data, reflected in a Z-Score of |Z| = 2.8 (Extended Data Fig. 6a). Trees where the Marianas were modelled to split with Amis.DG (suggesting settlement from Taiwan) or with Igorot.DG (from the Philippines) were rejected (|Z| = 10.88, |Z| = 11.303, respectively; Extended Data Fig. 6d,e). We fitted B10/WAT002 as a mixture of Papuan.DG and an Asian-related component deriving from different branches reflecting the patterns consistent with the respective hypothesis: assuming (1), B10/WAT002s Asian-related component was fixed to split as a sister group of the ERO genomes. In the alternative hypothesis (2), the Asian-related component of WAT002 was fit as branching after Ami.DG and Igorot.DG, resulting in sister groups of Guam/Saipan and ERO. The first (1) tree (Extended Data Fig. 6b) provides a better fit for the data, indicated by the worst Z-score of |Z| = 2.883. The fit of the second (2) tree (Extended Data Fig. 6c) is much worse (|Z| = 4.196) and shows a zero-drift branch from the common ancestor of ERO and the Marianas, and the ancestry present in WAT002, suggesting a poor fit. The model (3), suggesting the East-Asian-related ancestry arrived first in Remote Oceania before dispersing back to the Bismarck Islands and the Marianas, also shows a poor fit (|Z| = 4.196) (Extended Data Fig. 6f).
To assess the maternal lineage, the mtDNA-enriched sequences were processed in nf-core/eager version 2.4.0152 (https://nf-co.re/eager) using Nextflow version 21.04.3153. FastQC version 0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and fastP version 0.20.1. Fastp154 was used for sequencing quality control, adaptors were removed with AdapterRemoval version 2.3.2120. The remaining reads were aligned to the mitochondrial reference genome with circular mapper version 1.0119 and the resulting bam files filtered with SAMtools version 1.12122 for a minimum mapping quality and minimum read length of 30. Qualimap version 2.2.2-dev155 and bedtools version 2.30.0156 were used to generate mapping statistics and duplicates were removed using Picard MarkDuplicates version 2.26.0 (http://broadinstitute.github.io/picard/). Ancient DNA damage was assessed with DamageProfiler version 0.4.9157, endogenous DNA estimated using endorS.py version 0.4 (https://github.com/aidaanva/endorS.py) and the final report generated with MultiQC version 1.11158. From the resulting bam files, we built consensus sequences using the export function in Geneious version 2019.2.3159, setting ‘If no coverage call’ and to ‘Call ‘{}’ if coverage < ‘{}’’N/X with a coverage threshold of 5, aligned to ‘highest quality’ and >50% Sanger heterozygotes. Only NBR008 did not contain enough reads for sequence calling. For all other individuals, mitochondrial haplogroups were then determined with HaploGrep version 2.4.0160 Most individuals were assigned to ‘Austronesian’-associated haplogroups with 23 B4a1a*1, four B4a1a1a and one B4a1a1 haplogroup assignment. For NBR006 we identified haplogroup F1a3a, associated with present-day populations from the Philippines. For NBR017 the haplogroup P1, frequently found in Papua New Guinea and for NBR004 and NBR013 the Papuan-associated haplogroup Q1. Haplogroups could not confidently be assigned for NBR022, NBR008 and WAT006. Y-chromosomal haplogroups were identified by calling the SNPs covered on the Y chromosome of all male individuals by using the pileup from the Rsamtools161 package and assigning haplogroups by analysing the overlap with the ISOGG SNP index v.14.07 as detailed in (Extended Data Fig. 5)162.
To evaluate parental relatedness, runs of homzygosity (ROH) (Extended Data Fig. 8) were detected using hapROH v.0.3a4163 in Jupyter notebooks v.6.4.4. Using the provided down-sampled 1,000 Genomes data as a reference panel, the habs_ind command evaluated ROH on 22 chromosomes (chs=range(1,23)) with data-specific parameters (e_model = ’haploid’, p_model = ’Eigenstrat’, n_ref=2504, random_allele=True, readcounts=False, delete=False, logfile=True, combine=True) (Extended Data Fig. 8). The ROH were also used to estimate effective population size (Ne) with a maximum likelihood approach (MLE_ROH_Ne) after removing individuals with ROH segments >20 cM (indicative of very close parental relatedness).
The genetic relatedness of all individuals within a site was calculated using READ164 (Supplementary Table 14). Using this signal, we identified first- and second-degree relationships, illustrated in Extended Data Fig. 7. To corroborate these results and understand more distant genetic relatedness, genetic segments identical by descent (IBD blocks) were analysed. To prepare the data, ATLAS v.0.9 was used to call genotype likelihoods (method = MLE) for all positions from the 1,000 Genomes Phase 3 release after recalibrating the base quality scores (regions = ’88_mammals.epo_low_coverage.10M_GRCh37.masked.bed’) according to post-mortem damage (length = 50). The genotype likelihoods were then used for imputation with GLIMPSE v.1.0.0, followed by sub-setting to the 1,240 K SNPs. Imputation quality was assessed by counting the number of SNPs with genotype probabilities above 0.99 and excluding individuals with less than 600,000 high-quality SNPs (ERI003, ERI004, NBR001, NBR002, NBR005, NBR007, NBR008, NBR009, NBRO14, NBR017, NBR020, NBR022, NBR024, NBR026, WAT002, WAT005, WAT006). Finally, ancIBD v.0.2a2 (https://pypi.org/project/ancIBD/)165 was used to calculate IBD blocks per individual from the 1,240 K-extracted data using the recommended default settings. We filtered to IBD > 20 cm to detect relationships up to the sixth degree that typically share multiple such long IBD segments. The IBD segments were further integrated to a model calculating the split times between the two sites. Details on the calculations can be found in the Supplementary Information.
Ethics statement
Permissions for the samples included in this study were provided by the National Museum and Art Gallery of Papua New Guinea (NMAG) to H.B. for Nebira, Eriama and Watom and G.S. and D.G. for Tilu and Nunguri, including the ancient DNA processing. The scope of the study and the results were discussed with members of the NMAG before submission and shared with the Papua New Guinea Institute of Medical Research for consideration through their Institutional Review Board. Results of the study are communicated to the general public in Madang through personal communication of G.S. and materials made available to the NMAG, Madang University and various local high schools. The skeletal elements will be returned to the NMAG (Nunguri and Tilu) and to the Department of Anatomy, University of Otago, Aotearoa (Nebira, Eriama and Watom).
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Data availability
The raw data of the captured libraries are available at the European Nucleotide Archive (ENA) at https://www.ebi.ac.uk/ena/browser/view/PRJEB68153 under the accession number PRJEB68153. The genotypes of the newly published individuals can be sourced through the Poseidon framework community archive via Github at https://github.com/poseidon-framework/community-archive under 2024_NaegeleNatureEcologyEvolution. We ask users of this data to consider the recommendations on the use detailed in the Supplementary Information, ‘Ethical considerations for the analysis of ancient genomes’. The skeletal elements sampled will be returned to the National Museum and Art Gallery of Papua New Guinea in Port Moresby (Nunguri and Tilu) and to the Department of Anatomy, University of Otago, Aotearoa (Nebira, Eriama and Watom) by spring of 2025 to be curated with their respective skeletal assemblages.
Code availability
The code used in the site split analysis is deposited via GitHub at https://github.com/hyl317/two_island.
References
Summerhayes, G. R. et al. Human adaptation and plant use in highland New Guinea 49,000 to 44,000 years ago. Science 330, 78–81 (2010).
Groube, L., Chappell, J., Muke, J. & Price, D. A 40,000 year-old human occupation site at Huon Peninsula, Papua New Guinea. Nature 324, 453–455 (1986).
Allen, J., Gosden, C., Jones, R. & White, J. P. Pleistocene dates for the human occupation of New Ireland, northern Melanesia. Nature 331, 707–709 (1988).
Clarkson, C. et al. Human occupation of northern Australia by 65,000 years ago. Nature 547, 306–310 (2017).
David, B. et al. 45,610–52,160 years of site and landscape occupation at Nawarla Gabarnmang, Arnhem Land plateau (northern Australia). Quat. Sci. Rev. 215, 64–85 (2019).
Stoneking, M. & Delfin, F. The human genetic history of East Asia: weaving a complex tapestry. Curr. Biol. 20, R188–R193 (2010).
Oliveira, S. et al. Ancient genomes from the last three millennia support multiple human dispersals into Wallacea. Nat. Ecol. Evol. 6, 1024–1034 (2022).
Liu, Y.-C. et al. Ancient DNA reveals five streams of migration into Micronesia and matrilocality in early Pacific seafarers. Science 377, 72–79 (2022).
Posth, C. et al. Language continuity despite population replacement in Remote Oceania. Nat. Ecol. Evol. 2, 731–740 (2018).
Lipson, M. et al. Population turnover in Remote Oceania shortly after initial settlement. Curr. Biol. 28, 1157–1165. e1157 (2018).
Kayser, M. et al. Melanesian origin of Polynesian Y chromosomes. Curr. Biol. 10, 1237–1246 (2000).
Wollstein, A. et al. Demographic history of Oceania inferred from genome-wide data. Curr. Biol. 20, 1983–1992 (2010).
Bellwood, P. & Dizon, E. The Batanes archaeological project and the ‘out of Taiwan’ hypothesis for Austronesian dispersal. J. Austronesian Stud. 1, 1–33 (2005).
Bellwood, P., Fox, J. J. & Tyron, D. The Austronesians: Historical and Comparative Perspectives (ANU Press, 2006).
Gray, R. D., Drummond, A. J. & Greenhill, S. J. Language phylogenies reveal expansion pulses and pauses in Pacific settlement. Science 323, 479–483 (2009).
Summerhayes, G. et al. An early Lapita site on Emirau, New Ireland, PNG. J. Pac. Archaeol. 1, 62–75 (2010).
Kirch, P. V., Cruz, J. F. S. & Bauer, B. S. Talepakemalai: Lapita and Its Transformations in the Mussau Islands of Near Oceania (Cotsen Institute of Archeology Press at UCLA, 2021).
Kinaston, R. et al. Lapita diet in Remote Oceania: new stable isotope evidence from the 3000-year-old Teouma site, Efate Island, Vanuatu. PLoS ONE 9, e90376 (2014).
Kirch, P. V. On The Road Of The Winds: An Archaeological History of the Pacific Islands Before European Contact (Univ. of California Press, 2017).
Green, R. C. [Diversity, continuity and change in the history of the Bismarck Archipelago] The Lapita Cultural Complex: current evidence and proposed models. Bull. Indo-Pac. Prehist. Assoc. 11, 295-305 (1991).
Hammarström, H., et al. Glottolog 5.1. Zenodo https://doi.org/10.5281/zenodo.14006617 (2024).
Allen, J. Revisiting Papuan ceramic sequence changes: another look at old data. J. Archaeol. Anthropol. Soc. Victoria 4, 4–15 (2010).
Allen, J., Holdaway, S. & Fullagar, R. Identifying specialisation, production and exchange in the archaeological record: the case of shell bead manufacture on Motupore Island, Papua. Archaeol. Oceania 32, 13–38 (1997).
Dutton, T. Language and trade in central and south-east Papua. Aust. J. Anthropol. 11, 341 (1978).
Summerhayes, G. R. & Allen, J. Lapita writ small? Revisiting the Austronesian colonisation of the Papuan south coast. Terra Australis Vol. 26 Oceanic Explor. 26, 97–122 (2007).
Kirch, P. V. Lapita and Its Transformations in Near Oceania: Archaeological Investigations in the Mussau Islands, Papua New Guinea, 1985–88. (Archaeological Research Facility, University of California, 2001).
McNiven, I. J. et al. New direction in human colonisation of the Pacific: Lapita settlement of south coast New Guinea. Australian Archaeology 72, 1–6 (2011).
Shaw, B. et al. Frontier Lapita interaction with resident Papuan populations set the stage for initial peopling of the Pacific. Nature Ecology & Evolution 6, 11 (2022).
Gaffney, D., Summerhayes, G. R. & Mennis, M. A Lapita presence on Arop/Long Island, Vitiaz Strait, Papua New Guinea? Debating Lapita: Distribution, Chronology Soc. Subsistence 52, 115 (2019).
Lilley, I. Prehistoric exchange across the Vitiaz Strait, Papua New Guinea. Curr. Anthropol. 29, 513–516 (1988).
Terrell, J. E. & Schechter, E. M. Deciphering the Lapita code: the Aitape ceramic sequence and late survival of the ’Lapita face’. Cambridge Archaeol. J.l 17, 59 (2007).
Gaffney, D. et al. Earliest pottery on New Guinea mainland reveals Austronesian influences in highland environments 3000 years ago. PLoS ONE 10, e0134497 (2015).
Summerhayes, G. R. Austronesian expansions and the role of mainland New Guinea: a new perspective. Asian Perspect. 58, 250–260 (2019).
David, B. et al. A new ceramic assemblage from Caution Bay, south coast of mainland PNG: the linear shell edge-impressed tradition from Bogi 1. J. Pac. Archaeol. 3, 73–89 (2012).
Lilley, I. Flights of fancy: fractal geometry, the Lapita dispersal and punctuated colonisation in the Pacific. Isl. Inq.: Colonisation, Seafaring Archaeol. Marit. Landscapes 29, 75 (2008).
McNiven, I. J. et al. Terrestrial engagements by terminal Lapita maritime specialists on the southern Papuan coast. In Peopled landscapes: Archaeological and biogeographic approaches to landscapes (eds Haberle, S. G. & David, B.) 119–154 (2012).
Ross, M. Proto Oceanic and the Austronesian Languages of Western Melanesia (Dept. of Linguistics, Research School of Pacific Studies, 1988).
Oram, N. in Oral Tradition in Melanesia (eds Denoon, D. & Lacey, R.) 207–229 (Institute of Papua New Guinea Studies, 1981).
Swadling, P. A review of the traditional and archaeological evidence for early Motu, Koita and Koiari settlement along the central south Papuan coast. Oral Hist. 5, 37–43 (1977).
Rhoads, J. W. Prehistoric Papuan exchange Systems: The hiri and its antecedents. In The Hiri in History: Further Aspects of Long Distance Motu Trade in Central Papua (ed. Dutton, T.) 131–151 (Australian National Univeristy, 1982).
Irwin, G. Themes in the History of Coastal papua New Guinea. In Man And A Half: Essays In Pacific Anthropology And Ethnobiology In Honour Of Ralph Bulmer (ed. Pawley, A.) 503–510 (Polynesian Society, 1991).
Skelly, R. J. & David, B. Hiri: Archaeology of Long-Distance Maritime Trade Along the South Coast of Papua New Guinea (Univ. of Hawai ‘i Press, 2017).
Vilgalys, G. & Summerhayes, G. Do hiccups echo? Late Holocene interaction and ceramic production in southern Papua New Guinea. Asian Perspect., 61–88 (2016).
Bulmer, S. Prehistoric settlement patterns and pottery in the Port Moresby area. J. Papua New Guinea Soc. 5, 29–91 (1971).
Petchey, F., Allen, M. S., Addison, D. J. & Anderson, A. Stability in the South Pacific surface marine 14C reservoir over the last 750 years. Evidence from American Samoa, the southern Cook Islands and the Marquesas. J. Archaeolog. Sci. 36, 2234–2243 (2009).
Petchey, F. & Green, R. Use of three isotopes to calibrate human bone radiocarbon determinations from Kainapirina (SAC), Watom Island, Papua New Guinea. Radiocarbon 47, 181–192 (2005).
Petchey, F., Green, R., Jones, M. & Phelan, M. A local marine reservoir correction (∆R) for Watom Island, Papua New Guinea. N. Z. J. Archaeol. 26, 29–40 (2005).
Beavan-Athfield, N., Green, R., Craig, J. C., McFadgen, B. & Bickler, S. Influence of marine sources on 14C ages: isotopic data from Watom Island, Papua New Guinea inhumations and pig teeth in light of new dietary standards. J. R. Soc. N. Z. 38, 1–23 (2008).
Petchey, F. New evidence for a mid-to late-Holocene change in the marine reservoir effect across the South Pacific Gyre. Radiocarbon 62, 127–139 (2020).
Kinaston, R. L. et al. Late-Lapita diet and subsistence strategies on Watom Island, Papua New Guinea: new stable isotope evidence from humans and animals. Am. J. Phys. Anthropol. 157, 30–41 (2015).
Kinaston, R. L., Buckley, H. R., Gray, A., Shaw, B. & Mandui, H. Exploring subsistence and cultural complexes on the south coast of Papua New Guinea using palaeodietary analyses. J. Archaeol. Sci. 40, 904–913 (2013).
Kinaston, R. L., Bedford, S., Spriggs, M. & Buckley, H. in The Routledge Handbook of Bioarchaeology in Southeast Asia and the Pacific (eds Oxenham, M. & Buckley, H.) 427–461 (Routledge, 2016).
Shaw, B. J. et al. Migration and mobility at the Late Lapita site of Reber-Rakival (SAC), Watom Island using isotope and trace element analysis: a new insight into Lapita interaction in the Bismarck Archipelago. J. Archaeol. Sci. 37, 605–613 (2010).
Shaw, B. J., Buckley, H., Summerhayes, G., Stirling, C. & Reid, M. Prehistoric migration at Nebira, South Coast of Papua New Guinea: new insights into interaction using isotope and trace element concentration analyses. J. Anthropol. Archaeol. 30, 344–358 (2011).
Gaffney, D. et al. Tropical foodways and exchange along the coastal margin of northeastern New Guinea. J. Field Archaeol. 45, 498–511 (2020).
Smith, I. Terrestrial fauna from excavations at the Kainapirina {SAC} locality, Watom Island, Papua New Guinea. N. Z. J. Archaeol. 20, 137–147 (2000).
Tromp, M. et al. Exploitation and utilization of tropical rainforests indicated in dental calculus of ancient oceanic Lapita culture colonists. Nat. Hum. Behav. 4, 489–495 (2020).
Duong, N. T. et al. Complete human mtDNA genome sequences from Vietnam and the phylogeography of Mainland Southeast Asia. Sci. Rep. 8, 11651 (2018).
Delfin, F. et al. Complete mtDNA genomes of Filipino ethnolinguistic groups: a melting pot of recent and ancient lineages in the Asia-Pacific region. Eur. J. Hum. Genet. 22, 228–237 (2014).
Ko, A. M.-S. et al. Early Austronesians: into and out of Taiwan. Am. J. Hum. Genet. 94, 426–436 (2014).
Petchey, F. & Clark, G. Clarifying the age of initial settlement horizon in the Mariana Islands and the impact of hard water: a response to Carson (2020). Radiocarbon 63, 905–913 (2021).
Carson, M. T. Peopling of Oceania: vlarifying an initial settlement horizon in the Mariana Islands at 1500 BC. Radiocarbon 62, 1733–1754 (2020).
Hung, H.-C. et al. The first settlement of Remote Oceania: the Philippines to the Marianas. Antiquity 85, 909–926 (2011).
Winter, O., Clark, G., Anderson, A. & Lindahl, A. Austronesian sailing to the northern Marianas, a comment on Hung et al. (2011). Antiquity 86, 898–910 (2012).
Fitzpatrick, S. M. & Callaghan, R. T. Estimating trajectories of colonisation to the Mariana Islands, western Pacific. Antiquity 87, 840–853 (2013).
Pugach, I. et al. Ancient DNA from Guam and the peopling of the Pacific. Proc. Natl Acad. Sci. USA 118, e2022112118 (2021).
Bulmer, S. Prehistoric Culture Change In The Port Moresby region (Univ. of Papua New Guinea, 1978).
Groves, M. Western Motu descent groups. Ethnology 2, 15–30 (1963).
Petchey, F. et al. Testing the human factor: radiocarbon dating the first peoples of the South Pacific. J. Archaeolog. Sci. 38, 29–44 (2011).
Petchey, P., Buckley, H., Walter, R., Anson, D. & Kinaston, R. The 2008–2009 excavations at the SAC locality, Reber-Rakival Lapita site, Watom Island, Papua New Guinea. J. Indo-Pac. Archaeol. 40, 12–31 (2016).
Blackwood, B. & Danby, P. A study of artificial cranial deformation in New Britain. J. R. Anthropol. Inst. G.B. Irel. 85, 173–191 (1955).
Parkinson, R. Thirty Years in the South Seas: Land and People, Customs and Traditions in the Bismarck Archipelago and on the German Solomon Islands (Sydney Univ. Press, 2010).
Garanger, J. Incised and applied-relief pottery, its chronology and development in southeastern Melanesia, and extra areal comparisons. Stud. Oceanic Cult. Hist. 2, 53–66 (1971).
Specht, J. in Australian Natural History Vol. 17 (ed. Talbot, F. H.) (Australian Museum, 1973).
Lea, D. in Australian Natural History Vol. 17 (ed. Talbot, F. H.) (Australian Museum, 1973).
Spriggs, M. The Island Melanesians (Blackwell, 1997).
Neumann, K. Not the way it really was: writing a history of the Tolai (Papua New Guinea). J. Pac. Hist. 24, 209–220 (1989).
Bergström, A. et al. A Neolithic expansion, but strong genetic structure, in the independent history of New Guinea. Science 357, 1160–1163 (2017).
Gaffney, D. et al. Archaeological investigations into the origins of Bel trading groups around the Madang coast, northeast New Guinea. J. Isl. Coastal Archaeol. 13, 501–530 (2018).
Ross, M. A history of metatypy in the Bel languages. J. Lang. Contact 2, 149–164 (2008).
Gaffney, D. & Summerhayes, G. R. Coastal mobility and lithic supply lines in northeast New Guinea. Archaeol. Anthropol. Sci. 11, 2849–2878 (2019).
Gaffney, D. Materialising Ancestral Madang: Pottery Production and Subsistence Trading on the Northeast Coast of New Guinea (Univ. of Otago Studies in Archaeology, 2020).
Gaffney, D. & Summerhayes, G. R. An archaeological survey of inland Madang, northeast Papua New Guinea. Archaeology 12, 12–26 (2018).
Posth, C. et al. Response to ‘Ancient DNA and its contribution to understanding the human history of the Pacific Islands’ (Bedford et al. 2018). Archaeol. Oceania 54, 57–61 (2019).
Petchey, F., Spriggs, M., Bedford, S., Valentin, F. & Buckley, H. Radiocarbon dating of burials from the Teouma Lapita cemetery, Efate, Vanuatu. J. Archaeolog. Sci. 50, 227–242 (2014).
David, B. Rethinking cultural chronologies and past landscape engagement in the Kopi region, Gulf Province, Papua New Guinea. Holocene 18, 463–479 (2008).
Frankel, D. & Rhoads, J. W. (eds). Archaeology of a Coastal Exchange System: Sites and Ceramics of the Papuan Gulf (Australian National University, 1994).
Bickler, S. H. Early pottery exchange along the south coast of Papua New Guinea. Archaeol. Oceania 32, 151–162 (1997).
Skelly, R. J. From Lapita to the Hiri: Archaeology of the Kouri Lowlands, Gulf of Papua, Papua New Guinea (Monash Univ., 2014).
Lüning, S., Gałka, M., García-Rodríguez, F. & Vahrenholt, F. The Medieval climate anomaly in Oceania. Environ. Rev. 28, 45–54 (2020).
Lawman, A. E. et al. A century of reduced ENSO variability during the medieval climate anomaly. Paleoceanogr. Paleoclimatol. 35, e2019PA003742 (2020).
Shaw, B. et al. 2500-year cultural sequence in the Massim region of eastern Papua New Guinea reflects adaptive strategies to small islands and changing climate regimes since Lapita settlement. Holocene 30, 1075–1090 (2020).
Bulmer, S. Settlement and economy in prehistoric Papua New Guinea: a review of the archeological evidence. Journal de la Société des Océanistes 31, 7–75 (1975).
David, B. et al. Kumukumu 1, a hilltop site in the Aird Hills: implications for occupational trends and dynamics in the Kikori River delta, south coast of Papua New Guinea. Quat. Int. 385, 7–26 (2015).
Allen, J. Nebira 4: an early Austronesian site in central Papua. Archaeol. Phys. Anthropol. Oceania 7, 92–124 (1972).
Allen, J. Management of resources in prehistoric coastal Papua. Melanesian Environ., 35–44 (1977).
Allen, J., Golson, J. and Jones, R.M. (eds) Sunda and Sahul: prehistoric studies in southeast Asia, Melanesia and Australia (Academic Press, 1977).
Vanderwal, R. L. Prehistoric Studies in Central Coastal Papua (Australian National Univ., 1973).
Bulmer, S. West of Bootless Inlet: archaeological evidence for prehistoric trade in the Port Moresby area and the origins of the Hiri. Hiri Hist.: Further Aspects Long Distance Motu Trade Central Papua 11, 7–130 (1982).
Liu, D., Peter, B. M., Schiefenhövel, W., Kayser, M. & Stoneking, M. Assessing human genome-wide variation in the Massim region of Papua New Guinea and implications for the Kula trading tradition. Mol. Biol. Evol. 39, msac165 (2022).
Irwin, G., Shaw, B. & Mcalister, A. The origins of the Kula Ring: archaeological and maritime perspectives from the southern Massim and Mailu areas of Papua New Guinea. Archaeol. Oceania 54, 1–16 (2019).
Blust, R. Three recurrent changes in Oceanic languages. In Pacific Island Languages: Essays in Honour of G.B. Milner (ed. Davidson, J. H. C. S.) 7–28 (University of Hawaii Press, 1990).
Ross, M. D. Central Papuan culture history: some lexical evidence. In Austronesian Terminologies: Continuity and Change (ed. Ross, M. D.) 389–479 (Pacific Linguistics, 1994).
Dutton, T. E. Reconstructing Proto Koiarian: The history of a Papuan Language Family (Australian National University, 2010).
Dutton, T. E. The Peopling of Central Papua: Some Preliminary Observations (Ausralian National University, 1969).
Swadling, P. The settlement history of the Motu and Koita speaking people of the Central Province, Papua New Guinea. In Oral Tradition in Melanesia (eds Denoon, D. & Lacey, R.) 240–251 (University of Papua New Guinea, 1981).
Seligman, C. G., Barton, F. R. & Giblin, E. The Melanesians of British New Guinea (Cambridge Univ. Press, 1910).
Pinhasi, R. et al. Optimal ancient DNA yields from the inner ear part of the human petrous bone. PLoS ONE 10, e0129102 (2015).
Dabney, J. et al. Complete mitochondrial genome sequence of a Middle Pleistocene cave bear reconstructed from ultrashort DNA fragments. Proc. Natl. Acad. Sci. USA 110, 15758–15763 (2013).
Rohland, N. & Hofreiter, M. Ancient DNA extraction from bones and teeth. Nat. Protoc. 2, 1756–1762 (2007).
Rohland, N., Harney, E., Mallick, S., Nordenfelt, S. & Reich, D. Partial uracil-DNA-glycosylase treatment for screening of ancient DNA. Philos. Trans. R. Soc. Lond. B Biol. Sci. 370, 20130624 (2015).
Meyer, M. & Kircher, M. Illumina sequencing library preparation for highly multiplexed target capture and sequencing. Cold Spring Harbor Protoc. 2010, pdb prot5448 (2010).
Kircher, M., Sawyer, S. & Meyer, M. Double indexing overcomes inaccuracies in multiplex sequencing on the Illumina platform. Nucleic Acids Res. 40, e3 (2012).
Gansauge, M.-T. et al. Single-stranded DNA library preparation from highly degraded DNA usingT4DNA ligase. Nucleic Acids Res. https://doi.org/10.1093/nar/gkx033 (2017).
Gansauge, M.-T., Aximu-Petri, A., Nagel, S. & Meyer, M. Manual and automated preparation of single-stranded DNA libraries for the sequencing of DNA from ancient biological remains and other sources of highly degraded DNA. Nat. Protoc. 15, 2279–2300 (2020).
Fu, Q. et al. A revised timescale for human evolution based on ancient mitochondrial genomes. Curr. Biol. 23, 553–559 (2013).
Fu, Q. et al. An early modern human from Romania with a recent Neanderthal ancestor. Nature 524, 216–219 (2015).
Ghodsi, M., Liu, B. & Pop, M. DNACLUST: accurate and efficient clustering of phylogenetic marker genes. BMC Bioinf. 12, 271 (2011).
Peltzer, A. et al. EAGER: efficient ancient genome reconstruction. Genome Biol. 17, 60 (2016).
Schubert, M., Lindgreen, S. & Orlando, L. AdapterRemoval v2: rapid adapter trimming, identification, and read merging. BMC Res. Notes 9, 88 (2016).
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25, 1754–1760 (2009).
Li, H. et al. The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079 (2009).
Mallick, S. et al. The Simons genome diversity project: 300 genomes from 142 diverse populations. Nature 538, 201–206 (2016).
Skoglund, P. et al. Genomic insights into the peopling of the Southwest Pacific. Nature 538, 510–513 (2016).
Yang, M. A. et al. Ancient DNA indicates human population shifts and admixture in northern and southern China. Science, 282–288 (2020).
Posth, C. et al. Reconstructing the deep population history of Central and South America. Cell 175, 1185–1197 e1122 (2018).
Lipson, M. et al. Three phases of ancient migration shaped the ancestry of human populations in Vanuatu. Curr. Biol. 30, 4846–4856. e4846 (2020).
Fu, Q. et al. The genetic history of Ice Age Europe. Nature 534, 200–205 (2016).
Korneliussen, T. S. A. & Nielsen, A. Rasmus. ANGSD: analysis of next generation sequencing data. BMC Bioinf. 15, 13 (2014).
Renaud, G., Slon, V., Duggan, A. T. & Kelso, J. Schmutzi: estimation of contamination and endogenous mitochondrial consensus calling for ancient DNA. Genome Biol. 16, 224 (2015).
Key, F. M., Posth, C., Krause, J., Herbig, A. & Bos, K. I. Mining metagenomic data sets for ancient DNA: recommended protocols for authentication. Trends Genet. 33, 508–520 (2017).
Furtwangler, A. et al. Ratio of mitochondrial to nuclear DNA affects contamination estimates in ancient DNA analysis. Sci. Rep. 8, 14075 (2018).
Patterson, N. et al. Ancient admixture in human history. Genetics 192, 1065–1093 (2012).
Skoglund, P. et al. Genetic evidence for two founding populations of the Americas. Nature 525, 104–108 (2015).
Bergström, A. et al. Insights into human genetic variation and population history from 929 diverse genomes. Science 367, eaay5012 (2020).
Reimer, P. J. et al. IntCal13 and marine13 radiocarbon age calibration curves 0–50,000 years cal BP. Radiocarbon 55, 1869–1887 (2013).
Petchey, F. & Ulm, S. Marine reservoir variation in the Bismarck region: an evaluation of spatial and temporal change in ΔR and R over the last 3000 years. Radiocarbon 54, 45–58 (2012).
Smith, B. H. Standards of human tooth formation and dental age assessment. In Advances in Dental Anthropology (ed. Larsen, K.) 43–168 (Wiley-Liss, 1991).
Moorrees, C. F., Fanning, E. A. & Hunt, E. E. Jr Age variation of formation stages for ten permanent teeth. J. Dent. Res. 42, 1490–1502 (1963).
Fuller, B. T., Richards, M. P. & Mays, S. A. Stable carbon and nitrogen isotope variations in tooth dentine serial sections from Wharram Percy. J. Archaeol. Sci. 30, 1673–1684 (2003).
Longin, R. New method of collagen extraction for radiocarbon dating. Nature 230, 241–242 (1971).
Brown, T. A., Nelson, D. E., Vogel, J. S. & Southon, J. R. Improved collagen extraction by modified Longin method. Radiocarbon 30, 171–177 (1988).
Collins, M. & Galley, P. Towards an optimal method of archaeological collagen extraction: the influence of pH and grinding. Ancient Biomol. 2, 209–223 (1998).
Jacobi, R. M., Higham, T. F. G. & Bronk Ramsey, C. AMS radiocarbon dating of Middle and Upper Palaeolithic bone in the British Isles: improved reliability using ultrafiltration. J. Quat. Sci. 21, 557–573 (2006).
Koch, P. L., Tuross, N. & Fogel, M. L. The effects of sample treatment and diagenesis on the isotopic integrity of carbonate in biogenic hydroxylapatite. J. Archaeol. Sci. 24, 417–429 (1997).
Pin, C. & Bassin, C. Evaluation of a strontium-specific extraction chromatographic method for isotopic analysis in geological materials. Anal. Chim. Acta 269, 249–255 (1992).
Moore, L., Murphy, T., Barnes, I. & Paulsen, P. Absolute isotopic abundance ratios and atomic weight of a reference sample of strontium. J. Res. Natl Bur. Stand. 87, 1 (1982).
Burke, W. et al. Variation of seawater 87Sr/86Sr throughout Phanerozoic time. Geology 10, 516–519 (1982).
Patterson, N., Price, A. L. & Reich, D. Population structure and eigenanalysis. PLoS Genet. 2, e190 (2006).
Reich, D. et al. Reconstructing Native American population history. Nature 488, 370–374 (2012).
Harney, É., Patterson, N., Reich, D. & Wakeley, J. Assessing the Performance of qpAdm: a statistical tool for studying population admixture. Genetics 217, iyaa045 (2020).
Fellows Yates, J. A. et al. Reproducible, portable, and efficient ancient genome reconstruction with nf-core/eager. PeerJ 9, e10947 (2021).
Di Tommaso, P. et al. Nextflow enables reproducible computational workflows. Nat. Biotechnol. 35, 316–319 (2017).
Chen, S., Zhou, Y., Chen, Y. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, i884–i890 (2018).
Okonechnikov, K., Conesa, A. & García-Alcalde, F. Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics 32, 292–294 (2016).
Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics 26, 841–842 (2010).
Neukamm, J., Peltzer, A. & Nieselt, K. DamageProfiler: fast damage pattern calculation for ancient DNA. Bioinformatics 37, 3652–3653 (2021).
Ewels, P., Magnusson, M., Lundin, S. & Käller, M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32, 3047–3048 (2016).
Kearse, M. et al. Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649 (2012).
Weissensteiner, H. et al. HaploGrep 2: mitochondrial haplogroup classification in the era of high-throughput sequencing. Nucleic Acids Res. 44, W58–W63 (2016).
Morgan, M., Pagès, H., Obenchain, V., & Hayden, N. R: samtools: Binary alignment (BAM), FASTA, variant call (BCF), and tabix file import. R version 1.34.1 (2019); https://bioconductor.org/packages/release/bioc/html/Rsamtools.html
Rohrlach, A. B. et al. Using Y-chromosome capture enrichment to resolve haplogroup H2 shows new evidence for a two-path Neolithic expansion to Western Europe. Sci. Rep. 11, 15005 (2021).
Ringbauer, H., Novembre, J. & Steinrücken, M. Parental relatedness through time revealed by runs of homozygosity in ancient DNA. Nat. Commun. 12, 5425 (2021).
Monroy Kuhn, J. M., Jakobsson, M. & Gunther, T. Estimating genetic kin relationships in prehistoric populations. PLoS ONE 13, e0195491 (2018).
Ringbauer, H. et al. ancIBD-Screening for identity by descent segments in human ancient DNA. Preprint at bioRxiv (2023).
Acknowledgements
We thank the National Museum and Art Gallery of Papua New Guinea for the permission to analyse the individuals excavated from Watom, Nebira, Eriama Nunguri and Tilu. We specifically thank A. Kuaso for his consultations and input on the final results. We acknowledge the pioneering work of the late S. Bulmer. She led extensive surveys and excavations on the south coast and highlands of Papua New Guinea. We thank the Tilu community at Malmal village and the Watom communities for their assistance in the research. Many thanks to F. Aron and L. Papac for support in the wet lab. We thank S. Halcrow for assistance with the repatriation of the studied elements. This study was funded through the ‘WAVES’ starting grant, granted by the European Research Council (ERC) (#ERC758967) awarded to A.P., a University of Otago Research Grant, a University of Otago Doctoral Scholarship and the Max Planck Society.
Funding
Open access funding provided by Max Planck Society.
Author information
Authors and Affiliations
Contributions
R.K., H.B., D.G., G.S., M.T., P.P. and D.A. contributed archaeological material and contextualized the genetic findings. R.K. analysed stable isotope data with critical input from B.S., C.S., M.R. and D.B. F.P. performed the corrections for the marine reservoir effect. M.T. analysed microparticle data. K.N., E.B. and R.R. performed ancient DNA laboratory work. K.N., Y.H., A.B.R. and S.C. performed population genetic analysis. K.N. and R.K. interpreted the data with critical input from C.P., J.K., H.R., A.B.R., archaeological contextualization by D.G., G.S., M.T., R.K., B.S. and H.B. and linguistic contextualization by M.W. K.N. wrote the paper with critical input from C.P., R.K., J.K., A.P., G.S., D.G. and the other authors. K.N., Y.H., S.C. and R.K. produced the figures. A.P., G.S., H.B. and R.K. organized sample collection. A.P., J.K. and C.P. conceived of and coordinated the study.
Corresponding authors
Ethics declarations
Competing interests
The authors declare no competing interests.
Peer review
Peer review information
Nature Ecology & Evolution thanks João Teixeira and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Extended data
Extended Data Fig. 1 Isotope analysis.
(a) Comparison of mean ( ± 1 SD) human δ¹³Ccollagen and δ¹⁵Ncollagen values for each site; (b) Ward’s hierarchical cluster analysis of the δ¹³Ccollagen δ¹⁵Ncollagen values of the humans from all sites; (c) Strontium (87Sr/86Sr) results for the humans from Watom, compared to the human mean 87Sr/86Sr value (±2 SD) from the site (the non-Lapita burials [SAD1 and B12] are identified on the graph); (d) Strontium (87Sr/86Sr) results for the humans from Nebira compared to the mean wallaby 87Sr/86Sr value (±2 SD) (the non-local burials are identified on the graph). (e) Ward’s hierarchical cluster analysis of the δ¹³Ccollagen, δ¹⁵Ncollagen, δ¹³Cdentine, δ¹⁵Ndentine, δ¹³Cenamel, and 87Sr/86Sr values of the humans from Nebira (the red rectangle delineates the non-local individuals as identified from the strontium isotope analysis, excluding Burial 42, who did not provide δ¹³Ccollagen, δ¹⁵Ncollagen data). (f) δ¹³Cdentine and δ¹³Cenamel results for the Nebira and Watom humans and Nebira and Eriama wallabies are presented in reference to experimental dietary regression lines presented in 166.
Extended Data Fig. 3 f4-statistics.
(a) f4-scatterplot investigating the affinities to New Guinean highlanders (x-axis) and Baining from New Britain in the Bismarck Archipelago (y-axis). (b) f4-statistic showing affinity of the Asian component to ancient Taiwanese (Suogang, negative test scores) or Early Remote Oceanians (ERO) from Vanuatu and Tonga (positive test scores). Tests with overlapping SNPs < 1000 are not shown. Points indicate f4-value; black lines indicate one standard error determined by block jackknive.
Extended Data Fig. 4
LD decay curves for the Admixture dating models.
Extended Data Fig. 5 Uniparentally inherited markers.
Showing the LD decay between segments of Asian (Han.DG, Ami.DG, Atayal.DG, Igorot.DG, Kinh.DG, She.DG, Dai.DG) and Near Oceanic (Papuan.DG) ancestry in the individuals included in the analysis (Fig. 3d). (a) Mitochondrial haplogroups. (b) Y-chromosomal haplogroups.
Extended Data Fig. 6 qpGraph model.
Testing relative affinities between the included populations. The constructed base tree (a) with the best tree modelling Amis and Kankanaey as sister groups, rather than a split of the Mariana Island ancestry from Amis (d) or Kankanaey (e); testing hypothesis i: settlement of the Mariana islands directly from Island South East Asia (b); and the alternative hypothesis ii: settlement of the Marianas after the settlement of the Bismarck Archipelago by Austronesian-related groups (c); for completness, model iii was tested, to exclude a back dispersal from Remote Oceania (f).
Extended Data Fig. 7 Burial patterns and genetic relationships at the site of Nebira.
Position and orientation of graves from the early and late burial phase at Nebira (adapted from66) (a). Coloured graves indicate individuals with close genetic relationships. Individuals with black numbers indicate a close genetic relationship; individuals with grey numbers are included in the genetic analysis but do not indicate close genetic relations. Numbers correspond to burial numbers in66. Genetic relationship network (b), showing first (blue lines) and second degree (green lines) relationships, mitochondrial and Y-chromosomal haplogroups where available, morphological age and grave goods found with the individuals. (c) IDB block analysis for pairs of all individuals with successfully imputed genomes evaluating the total length of IBD blocks longer than 20 cM against the number of IBD blocks longer than 20 cm (d) and zoom in to the cluster of pairs with smaller sharing of blocks. Pairs are coloured by the pairing of sites as indicated in the legend. Pairs sharing IBD between the sites are labelled with the respective pair.
Extended Data Fig. 8 Sum of Runs of Homozygosity (ROH).
ROH of different sizes: 4–8 cm in dark blue, 8–12 cm in light blue, 12–20 cm in yellow and 20–300 cm in red.
Supplementary information
Supplementary Information
Modelling of the split time between Eriama and Nebira and Supplementary Tables 1–14 index.
Supplementary Tables 1–14
Supplementary Tables 1–14 containing the full test results and raw data described in the paper.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Nägele, K., Kinaston, R., Gaffney, D. et al. The impact of human dispersals and local interactions on the genetic diversity of coastal Papua New Guinea over the past 2,500 years. Nat Ecol Evol 9, 908–923 (2025). https://doi.org/10.1038/s41559-025-02710-x
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1038/s41559-025-02710-x
