A Genomic History of Aboriginal Australia

Dataset

We collected saliva samples for DNA sequencing in collaboration with Aboriginal Australian communities and individuals in Australia (S01). We sequenced genomes at high-depth (average of 60X, range 20X-100X) from 83 Aboriginal Australian individuals widely distributed geographically and linguistically (Figure 1, Extended Data Table 1, S02-S04). Additionally, we sequenced 25 Highland Papuan genomes (38X-53X; S03, S04) from individuals representative of five linguistic groups, and generated genotype data for 45 additional Papuans living or originating in the highlands (Figure 1). These datasets were combined with previously published genomes and SNP array genotype data, including Aboriginal Australian data from Arnhem Land, and from a human diversity cell line panel from the European Collection of Cell Cultures²⁰ (ECCAC, Figure 1, S04).

We explored the extent of admixture in the Aboriginal Australian autosomal gene pool by estimating ancestry proportions with an approach based on sparse nonnegative matrix factorization (sNMF)²¹. We found that the genomic diversity of Aboriginal Australian populations is best modelled as a mixture of four main genetic ancestries that can be assigned to four geographic regions based on their relative frequencies: Europe, East Asia, New Guinea and Australia (Figure 2a, Extended Data Figure 1, S05). The degree of admixture varies among groups (S05) with the Ngaanyatjarra speakers from central Australia (WCD) having a significantly higher “Aboriginal Australian component” (median value = 0.95) in their genomes than the other groups sampled (median value = 0.64; Mann-Whitney rank sum test, one tail p-value = 3.55e-07). The “East Asian” and “New Guinean” components are mostly present in northeastern Aboriginal Australian populations, while the “European component” is widely distributed across groups (Figure 2a, Extended Data Figure 1, S05). In most of the subsequent analyses, we either selected specific samples or groups according to their level of Aboriginal Australian ancestry, or masked the data for the non-Aboriginal Australian ancestry genomic components (S06).

Colonisation of Sahul

The origin of Aboriginal Australians is a source of much debate, as is the nature of the relationships among Aboriginal Australians, and between Aboriginal Australians and Papuans. Using f₃ statistics²², estimates of genomic ancestry proportions and classical multidimensional scaling (MDS) analyses, we find that Aboriginal Australians and Papuans are closer to each other than to any other present-day worldwide population considered in our study (Figure 2b, c, S05). This is consistent with Aboriginal Australians and Papuans originating from a common ancestral population which initially colonised Sahul. Moreover, outgroup f₃ statistics do not reveal any significant differences between Papuan populations (Highland Papuan groups and HGDP-Papuans) in their genetic affinities to Aboriginal Australians (Extended Data Figure 2a), suggesting that Papuans share a common ancestor after or at the same time as the divergence between Aboriginal Australians and Papuans.

To investigate the number of founding waves into Australia, we contrasted alternative models of settlement history through a composite likelihood method that compares the observed joint Site Frequency Spectrum (SFS) to that predicted under specific demographic models²³ (Figure 3, S07). We compared HGDP-Papuans to four Aboriginal Australian population samples with low levels of European admixture (Extended Data Figure 1) from both northeastern (CAI and WPA) and southwestern desert (WON and WCD) Australia. We compared one and two-wave models, where each Australian region was either colonised independently, or by descendants of a single Australian founding population after its divergence from Papuans. The one-wave model is a better fit to the observed SFS, suggesting that the ancestors of the sampled Aboriginal Australians diverged from a single ancestral population. This conclusion is also supported by MDS analyses (Figure 2b), as well as by estimation of ancestry proportion where all Aboriginal Australians form a cluster distinct from the Papuan populations (Extended Data Figure 1, S05). Additionally, it is supported by outgroup f₃ analyses, where all Aboriginal Australians are largely equidistant from Papuans when adjusting for recent admixture (Extended Data Figure 2b). Thus, our results, based on 83 Pama-Nyungan speakers, do not support earlier claims of multiple ancestral migrations into Australia giving rise to contemporary Aboriginal Australian diversity²⁴.

The SFS analysis indicates that there was a bottleneck in the ancestral Australo-Papuan population ~50 kya (95% CI 35-54 kya, S07), which overlaps with archaeological evidence for the earliest occupation of both Sunda and Sahul 47.5-55 kya^2,3. We further infer that the ancestors of Pama-Nyungan speakers and Highland Papuans diverged ~37 kya (95% CI 25-40 kya, Figure 3, S07), which is in close agreement with results from MSMC analyses (Extended Data Figure 2c, S08), a method estimating cross coalescence rates between pairs of populations based on individuals’ haplotypes²⁵. It is also in agreement with previous estimates, e.g., based on SNP array data¹⁸.

Archaic admixture

We characterised the number, timing and intensity of archaic gene flow events using three complementary approaches: SFS-based (S07), a goodness-of-fit analysis combining D-statistics (S09), and a method that infers putatively derived archaic ‘haplotypes’ (S10). Aboriginal Australians and Papuan genomes show an excess of putative Denisovan introgressed sites (Extended Data Figure 3a, S11), as well as substantially more putative Denisovan-derived haplotypes (PDHs) than other non-Africans (Extended Data Figure 3b, S10). The number and total length of those putative haplotypes vary considerably across samples. However, the estimated number of PDHs correlates almost perfectly (r² = 0.96) with the estimated proportion of Australo-Papuan ancestry in each individual (Extended Data Figure 3c). We found no significant difference in the distribution of the number of PDHs or the average length of PDHs between putatively unadmixed Aboriginal Australians and Papuans (Mann-Whitney U test, p>0.05). Moreover, the genetic differentiation between WCD and Papuans was also similar for both autosomal SNPs and PDHs with F_ST values around 0.12. Taken together, these observations provide evidence for Denisovan admixture predating the population split between Aboriginal Australians and Papuans (see also²⁶) and widespread recent Eurasian admixture in Aboriginal Australians (Figure 2a-b, S05). By constraining Denisovan admixture to have occurred before the Aboriginal Australian-Papuan divergence, the SFS-based approach results in an admixture estimate of ~4.0% (95% CI 3.3-5.0%, Figure 4, S07), similar to that obtained by D-statistics (~5%, S09). The SFS analyses further suggest that Denisovan/Australo-Papuan admixture took place ~44 kya (95% CI 31-50 kya, S07).

The SFS analysis also provides evidence for a primary Neanderthal admixture event (~2.3%, 95% CI 1.1-3.5%, Figure 4, S07) taking place in the ancestral population of all non-Africans ~60 kya (95% CI 55-84 kya, Figure 4, S07). Although we cannot estimate absolute dates of archaic admixture from the lengths of PDHs and putative Neanderthal-derived haplotypes (PNHs) in our samples, we can obtain a relative date. We found that for putatively unadmixed Aboriginal Australians and HGDP-Papuans, the average PNH and PDH lengths are 33.8 kb and 37.4 kb, respectively (Extended Data Figure 3b). These are significantly different from each other (p = 9.65 x10^-6 using a conservative sign test), and suggest that the time since Neanderthal admixture was about 11% greater than the time since Denisovan admixture, roughly in line with our SFS based estimates for the Denisovan pulse (31-50 kya) versus the primary pulse of Neanderthal admixture (55-84 kya). The SFS analysis also indicates that the main Neanderthal pulse was followed by a further 1.1% (95% CI 0.2-2.7%, Figure 4, S07) pulse of Neanderthal gene flow into the ancestors of Eurasians. Finally, using our SFS and haplotype based approaches, we explored additional models involving complex structure among the archaic populations. We found suggestive evidence that the archaic contribution could be more complex than a model involving the discrete Denisovan and Neanderthal admixture pulses shown in Figure 4^8,9 (S07, S10).

Out of Africa

To investigate the relationship of Australo-Papuan ancestors with other world populations, we computed D-statistics²² of the form ((H1=Aboriginal Australian,H2=Eurasian), H3=African) and ((H1=Aboriginal Australian,H2=Eurasian), H3=Ust’-Ishim). Several of these were significantly positive (S09), suggesting that Africans and Ust’-Ishim - a ~45 kya modern human from Asia²⁷ - are both closer to Eurasians than to Aboriginal Australians. These findings are in agreement with a model of Eurasians and Australo-Papuan ancestors dispersing from Africa in two independent waves. However, when correcting for a moderate amount of Denisovan admixture, Aboriginal Australians and Eurasians become equally close to Ust’-Ishim, as expected in a single OoA scenario (S09). Similarly, the D-statistics for ((H1=Aboriginal Australian,H2=Eurasian), H3=African) became much smaller after correcting for Denisovan admixture. Additionally, a goodness-of-fit approach combining D-statistics across worldwide populations indicates stronger support for two waves OoA, but when taking Denisovan admixture into account, a one-wave scenario fits the observed D-statistics equally well (Extended Data Figure 4a-b, S09).

To further investigate the timing and number of OoA events giving rise to present-day Australo-Papuans and Eurasians we used the observed SFS in a model-based composite likelihood framework. When considering only modern human genomes, we find evidence for two waves OoA, with a dispersal of Australo-Papuans ~14 ky before Eurasians (S07). However, when explicitly taking into account archaic Neanderthal and Denisovan introgression into modern humans^9,20, the SFS analysis supports a single origin for the OoA populations marked by a bottleneck ~72 kya (95% CI 60-104 kya, Figure 4, S07). This scenario is reinforced by the observation that the ancestors of Australo-Papuans and Eurasians share a 2.3% (95% CI 1.1-3.5%) Neanderthal admixture pulse. Furthermore, modern humans have both an LD decay rate and a number of predicted deleterious homozygous mutations (recessive genetic load) that correlates with distance from Africa (S05, S11, Extended Data Figure 5), again consistent with a single African origin. The model estimated from the SFS analysis also suggests an early divergence of Australo-Papuans from the ancestors of all non-Africans, in agreement with two colonisation waves across Asia^8,9,18. Under our best model, Australo-Papuans began to diverge from Eurasians ~58 kya (95% CI 51-72 kya, Figure 4, S07), whereas Europeans and East Asians diverged from each other ~42 kya (95% CI 29-55 kya, Figure 4, S07), in agreement with previous estimates^7,18,28. We find evidence for high levels of gene flow between the ancestors of Eurasians and Australo-Papuans, suggesting that, after the fragmentation of the OoA population (“Ghost” in Figure 4) 57-58 kya, the groups remained in close geographical proximity for some time before Australo-Papuan ancestors dispersed eastwards. Furthermore, we infer multiple gene flow events between sub-Saharan Africans and Western Eurasians after ~42 kya, in agreement with previous findings of gene flow between African and non-African populations²⁸.

MSMC analyses suggest that the Yoruba/Australo-Papuans and the Yoruba/Eurasians cross-coalescence rates are distinct, implying that the Yoruba and Eurasian gene trees across the genome have, on average, more recent common ancestors (Extended Data Figure 4c, S08). We show through simulations that these differences cannot be explained by typical amounts of archaic admixture (<20%, Extended Data Figure 4d). Moreover, the expected difference in phasing quality among genomes is not sufficient to fully explain this pattern (S08). While a similar separation in cross coalescence rate curves is obtained when comparing Eurasians and Australo-Papuans with Dinka, we find that, when comparing Australo-Papuans and Eurasians with San, the cross coalescence curves overlap (Extended Data Figure 4c). We also find that the inferred changes in effective population size through time of Aboriginal Australians, Papuans, and East Asians are very similar until around 50 kya, including a deep bottleneck around 60 kya (Extended Data Figure 6). Taken together, these MSMC results are consistent with a split of both Australo-Papuans and Eurasians from a single African ancestral population, combined with gene flow between the ancestors of Yoruba or Dinka (but not San) and the ancestors of Eurasians that is not shared with Australo-Papuans. These results are qualitatively in line with the SFS-based analyses (see e.g., Figure 4). While our results do not exclude the possibility of an earlier OoA expansion, they do indicate that any such event left little or no trace in the genomes of modern Australo-Papuans.

Genetic structure of Aboriginal Australians

Uniparental haplogroup diversity in this dataset (Extended Data Table 1, S12) is consistent with previous studies of mitochondrial DNA (mtDNA) and Y chromosome variation in Australia and Oceania²⁹, including the presence of typically European, Southeast and East Asian lineages³⁰. The combined results provide important insights into the social structure of Aboriginal Australian societies. Aboriginal Australians exhibit greater between-group variation for mtDNA (16.8%) than for the Y chromosome (11.3%), in contrast to the pattern for most human populations³¹. This result suggests higher levels of male- than female-mediated migration, and may reflect the complex marriage and post-marital residence patterns among Pama-Nyungan Australian groups³². As expected (S02), the inferred European ancestry for the Y chromosome is much greater than that for mtDNA (31.8% vs. 2.4%), reflecting male-biased European gene flow into Aboriginal Australian groups during the colonial era.

On an autosomal level, we find that genetic relationships within Australia reflect geography, with a significant correlation (r_GEN,GEO = 0.77, p-value < 0.0005, Extended Data Figure 7b) between the first two dimensions of an MDS analysis on masked genomes with geographical location (S13). Populations from the centre of the continent occupy genetically intermediate positions (Extended Data Figure 7a-b). A similar result is observed with an F_ST-based tree for the masked data (Extended Data Figure 7c, S05) as well as in analyses of genetic affinity based on f₃ statistics (Extended Data Figure 2a), suggesting a population division between northeastern and southwestern groups. This structure is further supported by SFS analyses showing that populations from southwestern desert and northeastern regions diverged as early as ~31 kya (95% CI 10-32 kya), followed by limited gene flow (estimated 2Nm<0.01, 95% CI 0.00<Nm< 11.25). The analysis of the major routes of gene flow within the continent supports a model in which the Australian interior acted as a barrier to migration. Using a model inspired by principles of electrical engineering where gene flow is represented as a current flowing through the Australian continent and using observed F_ST values as a proxy for resistance, we infer that gene flow occurred preferentially along the coasts of Australia (Extended Data Figure 7e-g, S13). These findings are consistent with a model of expansion followed by population fragmentation when the extreme aridity in the interior of Australia formed barriers to population movements during the LGM³³.

We used MSMC on autosomal data and mtDNA Bayesian Skyline Plots³⁴ (BSP) to estimate changes in effective population size within Australia. The MSMC analyses provide evidence of a population expansion starting ~10 kya in the northeast, while both MSMC and BSP indicate a bottleneck in the southwestern desert populations taking place during the past ~10 kya (Extended Data Figure 6, S08, S12). This is consistent with archaeological evidence for a population expansion associated with significant changes in socio-economic and subsistence strategies in the Holocene³⁵ that impacted Australia.

European admixture almost certainly had not occurred before the late 18^th century, but earlier East Asian and/or New Guinean gene flow into Australia could have taken place. We characterised the mode and tempo of gene flow into Aboriginal Australians using three different approaches (S06, S07, S14). We used approximate Bayesian computation (ABC) to compare the observed mean and variance in the proportion of European, East Asian and Papuan admixture among Aboriginal Australian individuals, to that computed from simulated datasets under various models of gene flow. We estimated European and East Asian admixture to have occurred on the order of ten generations ago (S14), consistent with historical and ethnographic records. Consistent with this, the local ancestry approach suggests that European and East Asian admixture is more recent than Papuan admixture (Extended Data Figure 4a, S06). In addition, both ABC and SFS analyses indicate that the best-fitting model for the Aboriginal Australian-Papuan data is one of continuous but modest gene flow, mostly unidirectional from Papuans to Aboriginal Australians, and geographically restricted to northeast Aboriginal Australians (2Nm=0.41, 95% CI 0.00-20.35, Figure 3, S07).

To further investigate gene flow from New Guinea, we conducted analyses on the Papuan ancestry tracts obtained from the local ancestry analysis. We inferred local ancestry as the result of admixture between four components: European, East Asian, Papuan and Aboriginal Australian (S06). Papuan tract length distribution shows a clear geographic pattern (Extended Data Figure 8); we find a strong correlation of Papuan tract length variance with distance from WCD to other Aboriginal Australian groups (r=0.64, p-value<0.0001). The prevalence of short ancestry tracts of Papuan origin, compared to longer tracts of East Asian and European origin, suggests that a large fraction of the Papuan gene flow is much older than that from Europe and Asia, which is consistent with the ABC analysis (S14). We also investigated possible South Asian (Indian related) gene flow into Aboriginal Australians, as reported recently¹⁸. However, we found no evidence of a component that can be uniquely assigned to Indian populations in the Aboriginal Australian gene pool using either admixture analyses or f₃ and D-statistics (S05), even when including the original Aboriginal Australian genotype data from Arnhem Land. The different size and nature of the comparative datasets may account for the discrepancy in the results.

Pama-Nyungan languages and genetic structure

To investigate whether or not linguistic relationships reflect genetic relationships among Aboriginal Australian populations, we built a Bayesian phylogenetic tree for the 28 different Pama-Nyungan languages represented in this sample¹³ (Extended Data Table 1, S15). The resulting linguistic and F_ST-based genetic trees (Extended Data Figure 7c,d) share several well-supported partitions. For example, both trees indicate that the northeastern (CAI and WPA) and southwestern groups (ENY, NGA, WCD and WON) form two distinct clusters, while PIL, BDV and RIV are intermediate. A distance matrix between pairs of languages, computed from the language-based tree, is significantly correlated with geographic distances (r_GEO,LAN = 0.83, Mantel test two-tail p-value on 9,999 permutations = 0.0001, S13). This suggests that differentiation among Pama-Nyungan languages in Australia follows geographic patterns, as observed in other language families elsewhere in the world⁶. Furthermore, we find a correlation between linguistics and genetics (r_GEN,LAN= 0.43, Mantel test p-value < 0.0005, S13) that remains significant when controlling for geography (r_GEN,LAN.GEO= 0.26, partial Mantel test p-value < 0.0005, S13). This is consistent with language differentiation after populations loose (genetic) contact with one another. The correlation between the linguistic and genetic trees is all the more striking given the difference in time scales: the Pama-Nyungan family is generally accepted to have diversified within the last 6 ky³⁶, while the genetic estimates are two to five times that age. The linguistic tree thus cannot simply reflect initial population dispersals, but rather reflects a genetic structure that has a complex history, with initial differentiation 10-32 kya, localised population expansions (northeast) and bottlenecks (southwest) ~10 kya, and subsequent limited gene flow from the northeast to the southwest. The latter may be the genetic signature that tracks the divergence of the Pama-Nyungan language family.

Selection in Aboriginal Australians

To identify selection signatures specific to Aboriginal Australians, we used two different methods based on the identification of SNPs with high allele frequency differences between Aboriginal Australians and other groups, similar to the Population-Branch Statistics³⁷ (PBS, S16). First, we scanned the Aboriginal Australian genomes for loci with unusually large changes in allele frequency since divergence from Papuans, taking recent admixture with Europeans and Asians into account (“global scan”). Second we identified genomic regions showing high differentiation associated with different ecological regions within Australia (“local scan”, S16). Among the top ranked peaks (Extended Data Table 2) we found genes associated with the thyroid system (NETO1, 7^th peak in the global scan, and KCNJ2, 1^st peak in the local scan) and serum urate levels (8^th peak in the global scan). Thyroid hormone levels are associated with Aboriginal-Australian-specific adaptations to desert cold³⁸ and elevated serum urate levels with dehydration³⁹. These genes are therefore candidates for potential adaptation to life in the desert. However, further studies are needed to associate putative selected genetic variants with specific phenotypic adaptations in Aboriginal Australians.

Discussion

Australia has one of the longest histories of continuous human occupation outside Africa, raising questions of origins, relatedness to other populations, differentiation and adaptation. Our large scale genomic data and analyses provide some answers but also raise new questions. We find that Aboriginal Australians and Eurasians share genomic signatures of an OoA dispersal - a common African ancestor, a bottleneck and a primary pulse of Neanderthal admixture. However, Aboriginal Australian population history diverged from that of other Eurasians shortly after the OoA event, and included private admixture with another archaic hominin.

Our genetic-based time estimates are relative, and to obtain absolute dates we relied on two rescaling-parameters: the human mutation rate and generation time (assumed to be 1.25×10^-8/generation/site and 29 years, respectively, based on recent estimates^40,41). Although the absolute estimates we report would need to be revised if these parameters were to change, the current values can be the starting point of future research and should be contextualized.

We find a relatively old divergence between the ancestors of Pama-Nyungan speakers and Highland Papuans, only ~10% younger than the European-East Asian split time. With the assumed rescaling-parameters this corresponds to ~37 kya (95% CI 25-40 kya) implying that the divergence between sampled Papuans and Aboriginal Australians is older than the disappearance of the land bridge between New Guinea and Australia about 7-14.5 kya, and thus suggests ancient genetic structure in Sahul. Such structure may be related to palaeo-environmental changes leading up to the LGM. Sedimentary studies show that the large Lake Carpentaria (500 x 250 km, Figure 1) formed ~40 kya, when sea-levels fell below the 53 m-deep Arafura Sill⁴². Although Australia and New Guinea remained connected until the early Holocene, the flooding of the Carpentaria basin and its increasing salinity^{46 24 45}may have thus promoted population isolation.

Our results also indicate that the population that diverged from Papuans ~37 kya was ancestral to all Aboriginal Australian groups sampled in this study; yet, archaeological evidence shows that by 40-45 kya, humans were widespread within Australia (Figure 1). Three non-exclusive scenarios could account for this observation: 1) the Aboriginal Australian ancestral population was widespread prior to the divergence from Papuans, maintaining gene flow across the continent; 2) it was deeply structured, and only one group survived to give rise to modern Aboriginal Australians; and 3) other groups survived, but the descendants are not represented in our sample. Additional modern genomes, especially from Tasmania and the non-Pama-Nyungan regions of the Northern Territory and Kimberley, as well as ancient genomes pre-dating European contact in Australia and other expansions across Southeast Asia¹⁷, may help resolve these questions in the future.

Data access

The Aboriginal Australian whole genome sequence data and the SNP array data generated in this study are available upon request by contacting the research ethics and integrity service at Griffith University (rick.williams@griffith.edu.au). The Papuan whole genome sequence data generated in this study are available under managed access through the EGA database (https://www.ebi.ac.uk/ega) under study accession number EGAS00001001247.

Extended Data

Supplementary Material