Ancient mitochondrial DNA provides high-resolution time scale of the peopling of the Americas.

The exact timing, route, and process of the initial peopling of the Americas remains uncertain despite much research. Archaeological evidence indicates the presence of humans as far as southern Chile by 14.6 thousand years ago (ka), shortly after the Pleistocene ice sheets blocking access from eastern Beringia began to retreat. Genetic estimates of the timing and route of entry have been constrained by the lack of suitable calibration points and low genetic diversity of Native Americans. We sequenced 92 whole mitochondrial genomes from pre-Columbian South American skeletons dating from 8.6 to 0.5 ka, allowing a detailed, temporally calibrated reconstruction of the peopling of the Americas in a Bayesian coalescent analysis. The data suggest that a small population entered the Americas via a coastal route around 16.0 ka, following previous isolation in eastern Beringia for ~2.4 to 9 thousand years after separation from eastern Siberian populations. Following a rapid movement throughout the Americas, limited gene flow in South America resulted in a marked phylogeographic structure of populations, which persisted through time. All of the ancient mitochondrial lineages detected in this study were absent from modern data sets, suggesting a high extinction rate. To investigate this further, we applied a novel principal components multiple logistic regression test to Bayesian serial coalescent simulations. The analysis supported a scenario in which European colonization caused a substantial loss of pre-Columbian lineages.

INTRODUCTION

The geographic isolation of the Americas delayed human settlement until the end of the Pleistocene [20 to 10 thousand years ago (ka)]; however, despite this relatively recent date, the specific time, place, and route of entry remain uncertain. It is likely that the first peoples moved from Asia across the Bering Land Bridge (, ), the landmass between Eurasia and America exposed by lowered sea levels during the Last Glacial Maximum (LGM). However, at this time, much of northern North America was covered by the Cordilleran and Laurentide ice sheets, which blocked access from eastern Beringia (Bering Land Bridge and Alaska/Yukon) southward to the rest of the Americas (Fig. 1A). Shortly after the Cordilleran ice sheet began to retreat ~17 ka (), a potential Pacific coastal route became available ~15 ka (Fig. 1B) (, ), whereas an alternative route through an inland ice-free corridor along the eastern side of the Rocky Mountains opened around ~11.5 to 11 ka (–). The timing and route used in the migration event are important in understanding the size, number, and speed of the first migratory wave(s). Timing and route are also pivotal in resolving contentious issues such as the nature of peoples before Clovis—the first widespread archaeologically recognized culture in North America (13.2 to 12.8 ka) ().

Fig. 1

Eastern Beringia during the LGM and retreat of the ice sheets.

(A) Exposed land when sea levels were lowest (light green), modern-day landmass (dark green), and ice sheets (white). At the height of the LGM, the Laurentide and Cordilleran ice sheets blocked access to the Americas from eastern Beringia (that is, the Bering Land Bridge and Alaska/Yukon) (). Populations west of the Bering Land Bridge were able to migrate southward during the LGM, but those on the Bering Land Bridge were unable to retreat farther than the Aleutian ice belt (arrows). The last point of detectable gene flow between Siberian and Native American ancestral populations (24.9 ka) and the geographic isolation marked by the formation of Native American founder lineages (18.4 ka) are shown (see Fig. 2B for details). The Yana Rhinoceros Horn site (32 ka) and the Swan Point site (14 ka) illustrate the temporal and geographic gaps in the Beringian archaeological record. (B) The ice sheets that began to retreat ~17 ka, opening a potential Pacific coastal route by ~15 ka (arrow). The rapid population expansion (16.0 ka) likely marks the movement south of the ice (see Fig. 3C for details).

Fig. 2

Comparison of Bayesian estimates of the TMRCA of the Native American founder haplogroups and of the divergence from Siberian lineages.

(A) Mean age (symbols) and 95% highest posterior density (HPD) (error bars) for the TMRCA of each of the Native American haplogroups. Shading indicates the period between the oldest lower bound of any 95% HPD and the youngest upper bound of any 95% HPD for each data set. The purple dotted lines show the TMRCA bounds based on tip calibration; the blue dotted lines show the extreme TMRCA bounds from previous publications (26.3 to 9.7 ka) (, ). (B) The isolation of Native American populations estimated to have occurred after the last observable divergence between Siberian and Native American lineages (24.9 ka based on the lowest 95% HPD upper bound) and before the oldest date at which all Native American founder haplogroups formed (18.4 ka based on the lowest 95% HPD upper bound). See section S5 for detailed methods.

Fig. 3

Dated Bayesian mitogenomic tree and reconstruction of past effective female population size.

The mitogenomic tree and the demographic plot are based on replicate data set 1, which is representative of the three replicate data sets (fig. S7). (A) Complete tree showing the relationships between the main Native American haplogroups A, B, C, and D, as well as their TMRCA (colored circles). Black circles show the divergences between Siberian and Native American lineages. Siberian clades are shown in black and Native American clades are shown in gray. (B) Detailed tree with Siberian clades (black), modern Native Americans (blue), and ancient Native Americans (red). Colored and black circles as in (A). Gray shadings and empty black circles highlight shared ancestry for individuals from the same geographic location or from the same cultural background. The filled black triangle (haplogroup A2) is the most recent common ancestor between an ancient haplotype and a modern haplotype at ~9 ka. (C) Extended Bayesian skyline plot of female effective population size, based on a generation time of 25 years.

Genetic studies of Native American populations are complicated by the demographic collapse and presumed major loss of genetic diversity following European colonization at the end of the 15th century (). However, geographically widespread signals of low diversity and shared ancestry (–)—particularly striking in maternally inherited mitochondrial and paternally inherited Y-chromosome sequence data—suggest that small founding groups possibly initially entered the Americas in a single migration event that gave rise to most of the ancestry of Native Americans today (, , ). In contrast, the distribution of some of the rare founding mitochondrial haplogroups (D4h3a along the Pacific coast of North and South America, and X2a in northwestern North America) suggests that distinct migrations along the coastal route and the ice-free corridor occurred within less than 2000 years (). Recent studies have identified a weak Australasian genomic signature in several Native American groups from the Amazon, compatible with two founding migrations (), although the Australasian gene flow may have occurred after the initial peopling (). Irrespective of the number of migration waves, the founding population appears to have rapidly grown and expanded southward (, , ), with low levels of gene flow between areas following initial dispersion (, ).

Unfortunately, the precision of molecular clock studies in the Americas to date has been limited by the low genetic diversity and lack of appropriate calibration points to accurately estimate rates of molecular evolution. As a result, current mitochondrial molecular clock estimates of the initial entry into the Americas, which assume that the event corresponds to the initial diversification of Native American genetic lineages, range from 26.3 to 9.7 ka (Fig. 2A). This broad range spans most of the time frame over which the Bering Land Bridge route was feasible. Given the narrow temporal span of the actual diversification and migration into the Americas, much greater precision is needed to distinguish between different migration routes and hypotheses.

Molecular clocks highly depend on the quality of calibration points to accurately estimate rates of molecular evolution (). In the Americas, the scarce evidence of early human occupation and the absence of sites in eastern Beringia for most of the late Pleistocene hinder reliable calibration. An additional major challenge is the temporal dependence of molecular rate estimates, whereby molecular evolution appears more rapid when measured over short time intervals (). In humans, this problem is most apparent when recent time scales (for example, the human settlement of the Americas) are analyzed using deep fossil calibrations such as the human-chimpanzee split ~6 to 7 million years ago (). Accurate molecular rate estimates require a distribution of calibration points close to the age of events under study (, ); in this regard, ancient DNA sequences from dated skeletons provide suitable tip calibrations for studying recent evolutionary events ().

To generate a detailed view of human matrilineal genetic diversity in the Americas through time, with multiple calibration points, we sequenced the complete mitochondrial genome (mitogenome) of 92 pre-Columbian individuals, ranging in age from ~8.6 to 0.5 ka (table S2). We then used Bayesian approaches to estimate coalescence times, reconstruct a demographic history, and simulate and test population scenarios.

RESULTS AND DISCUSSION

Novel genetic diversity in pre-Columbian times

The 92 pre-Columbian mitogenomes were sequenced to an average coverage depth of 112× (5.6× to 854.2×; table S2). Sequences were assigned to 84 distinct haplotypes, which fell within the expected overall mitochondrial diversity of Native South Americans (), that is, haplogroups A2, B2, C1b, C1c, C1d, and D1 (figs. S2 to S5). The Native South American haplogroup D4h3a was not observed in our ancient data set, although we sampled the South American southern cone (Arroyo Seco 2, Argentina) where this lineage is common today (). None of the 84 haplotypes identified from ancient samples are represented in the existing genealogy of global human mitochondrial diversity [that is, PhyloTree mt; ()] (figs. S2 to S5) or in the literature (fig. S6). Although modern Native American genetic diversity is not well characterized, this result clearly illustrates the importance of sampling pre-Columbian specimens to fully measure the past genetic diversity and to reconstruct the process of the peopling of the Americas.

Marked synchronicity of the Native South American haplogroup times to most recent common ancestor

The estimated times to most recent common ancestor (TMRCA) for haplogroups A2, B2, C1, D1, and D4h3a were highly synchronous (Fig. 3 and fig. S8), confirming previous interpretations that all five haplogroups were part of one initial population (). The TMRCA fell within the range of previous molecular date estimates, although the narrower 95% credible intervals considerably increased the precision (Fig. 2A). Older dates for the initial diversification within each haplogroup have been previously calculated using the human-chimpanzee calibration (, ), whereas much younger dates resulted from calibrations using non–Native American mitochondrial lineages associated with biogeographic events (). This clearly illustrates the impact of time-dependent rate estimates and the critical influence of the calibration framework (). Our use of a large number of temporally and phylogenetically distributed tip dates provides an accurate calibration of the molecular rates relevant to Native American early history (), allowing a uniquely precise timeline for the peopling of the Americas.

Separation from Siberian populations during the LGM

The most recent genetic divergence observed between the ancestors of Siberians and Native Americans (24.9 ka; section S5 and Fig. 2B) is the last point at which we can detect apparent gene flow (that is, a shared lineage) between the Siberian population and the ancestral Native American population. We can assume that the real population divergence occurred after this point. In addition, if we accept that the estimated TMRCA of each of the five Native American haplogroups provides an independent estimate of the timing of the same small population’s isolation, we can use the 95% credible intervals to constrain the lower bound (section S5). The resulting estimate that the two populations became fully isolated between 24.9 and 18.4 ka is in accordance with calculations from modern complete genomes which indicate that Siberians and Native Americans split no later than ~23 ka (). Gene flow to and from east Siberia certainly appears to have ceased by the height of the LGM (18.4 ka; Fig. 2B).

Eastern Beringia as a sustainable refugium for ancestral Native Americans

Our data cannot determine whether the separation between Siberian and Native American ancestral populations occurred in Siberia or Beringia. However, the start of isolation (24.9 to 18.4 ka) closely coincides with the LGM. We hypothesize that cold arid conditions drove populations on the western (that is, Siberian) margins of the Bering Land Bridge to migrate to southern refugia (Fig. 1A), as suggested by the absence of megafauna kill sites younger than the far north Yana Rhinoceros Horn site 32 ka (). In contrast, any populations east of the Kamchatka and Chukotka Peninsulas would not have been able to retreat farther south than the Aleutian ice belt and would thus remain isolated in eastern Beringia (Fig. 1A). We cannot accurately estimate the size of this founding population, but the effective female population that subsequently entered the Americas appears to be ~2000, which accords well with previous studies (, , ). Although this number cannot be directly translated into census population size, it suggests that the human population isolated in eastern Beringia was relatively small, probably not exceeding a few tens of thousands of people (section S6). The presence of large numbers of megafauna in eastern Beringia during the late Pleistocene, including the LGM, indicates an ice-free region dominated by shrub tundra (), which would have been more than capable of sustaining such a population size (section S6). Thus, our observations are consistent with the idea that the founding Native American population used the exposed Bering Land Bridge and adjacent regions in Alaska/Yukon as a refugium during the height of the LGM, before climatic change and the retreat of the ice sheets allowed access to the remainder of the Americas. Unfortunately, the large temporal and geographic gaps in the archaeological record between the Yana Rhinoceros Horn site (~32 ka, western Beringia) and the Swan Point site (~14 ka, eastern Beringia) provide little additional information about this process (Fig. 1A) () or how the ancestral Native Americans were isolated from their Asian counterparts.

The Beringian Standstill (~2.4 to 9 thousand years)

The scenario of an Eastern Beringia refugium is consistent with the Beringian standstill hypothesis, which suggests that the ancestral Native Americans were isolated in the area for up to 15 thousand years (ky) (, , ). Our large data set of dated mitogenomes provides tight estimates for the duration of the standstill and the subsequent movement out of the area. The mitogenomic tree shows a sudden burst of lineage diversification starting ~16.0 to 13.0 ka (Fig. 3B). This is followed by a steep increase in the mean female effective population size (>10%) between adjacent time bins starting 16.0 ka (Fig. 3C). Overall, the population underwent a 60-fold increase between 16.0 and 13.0 ka, suggesting that 16.0 ka represents the initial entry into the Americas, where population size significantly increased in a more favorable environment. Considering the time between isolation (24.9 to 18.4 ka) and entry (16.0 ka), the improved temporal resolution provided by our data suggests that the Beringian Standstill could have been as short as ~2.4 ky, and no longer than ~9 ky, consistent with recent estimates based on autosomal data from complete modern-day genomes ().

A coastal entry route

The population burst at 16.0 ka is contemporaneous with the rapid retreat of coastal glaciers along the northwest Pacific coast associated with a phase of stepwise ocean warming (2° to 3°C) in the region (). This date considerably predates the opening of the inland ice-free corridor ~11.5 to 11.0 ka (–) and indicates that the initial entry into the Americas took place via a southward expansion along recently emerged northwest Pacific coastal land (Fig. 1B) (, , , , ). Given the early archaeological sites in Monte Verde in southern Chile at 14.6 ka (), the mitogenome data indicate that the transit of the full length of the Americas took around 1.4 ky.

Early geographic structure after entry into the Americas

The phylogenetic trees feature multiple long branches stemming from the initial expansion period (Fig. 3B), irrespective of whether the lineages lead to pre-Columbian (ancient) or modern-day individuals. This topology appears to reflect the swift migration and expansion of a population, which contained each major haplogroup, into and throughout the Americas (, , ). Subsequent lineage diversification within each haplogroup appears to be constrained to within specific geographic regions or shared cultural backgrounds (Fig. 3B, gray shadings), which is consistent with suggestions that geographic structure was rapidly established after colonization and was thereafter followed by limited gene flow between populations from diverse regions (, ).

Extinction of ancient lineages after European colonization

It has been suggested that European colonization resulted in a bottleneck in Native American genetic diversity (, ). In the present data set, it is notable that no ancient haplotype (Fig. 3B, red lineages) shared a common ancestor with a modern haplotype (Fig. 3B, blue lineages) more recently than ~9.0 ka (Fig. 3B, black triangle in haplogroup A2), despite the number of samples examined from both populations. Our ancient samples were principally derived from large population centers along the western coast of South America, which experienced high extinction rates following European colonization. Historic accounts have reported that the population decline was more rapid and intense in the Gulf of Mexico and the Pacific coast of Peru than in other areas such as the Mesoamerican plateau or the Andean highlands (). The pronounced phylogeographic structure observed in Native American mitogenomes suggests that such demographic events could have removed major portions of genetic diversity in certain areas. As a result, it is possible that the lack of overlap between ancient and modern haplotypes is influenced by the undersampling of modern-day diversity in the source regions of the archaeological material.

To further examine the issue of local extinction, we used BayeSSC () to simulate 15,000 coalescence trees for each of seven potential population scenarios (section S7; Fig. 4 and fig. S10). The Akaike information criterion for model selection, which is included in the BayeSSC package, informs about the relative fit of the models but does not inform about their absolute quality. As an alternative, we applied a principal components multiple logistic regression (PCMLR) model () to the BayeSSC summary statistics to predict which of the population scenarios best fits the observed data. We introduce this novel approach in approximate Bayesian computation as a powerful way to estimate both the relative [graphical representation of principal components analysis (PCA)] (Fig. 4 and fig. S10) and the absolute (Bayesian model verification; table S6) qualities of the tested models (section S7).

Fig. 4

Effects of population structure and European colonization on South American mitochondrial diversity.

Seven population scenarios simulated with BayeSSC are represented (models A to G). The arrowheads in demographic models C to G represent the time of separation (9 ka) between the population carrying modern haplotypes (Population 0) and the population carrying ancient haplotypes (Population 1) (that is, meta-demes with a pattern of localized small-scale separation across the continent). The bottom-right panel is a PCA plot of summary statistics for the 15,000 simulations for each of the seven models (25% of simulated data are reported; see fig. S10 for the full data set). Observed data from the three replicates are shown as black circles and fall closest to results from model C (green points).

We found that demographic models with one panmictic population (model A in Fig. 4 and fig. S10) or with a simple geographic separation of populations harboring the ancient and modern lineages (that is, a localized small-scale separation across the continent; model D) were poor fits for the data (table S6). The addition of a post-European bottleneck with survival of both ancient and modern lineages (models B and E) did not improve the goodness of fit. More complex models combined geographic separation and early and ongoing loss of ancient lineages—accelerated (model F) or not (model G) by the European landfall—to account for demographic processes such as serial founder effects and drift, or population turnover and elimination during the expansion of Central Andean empires. Such models were also poor fits for the data (table S6). Only a model that combined both the geographic separation of populations harboring modern and ancient haplotypes and the subsequent rapid extinction of the ancient lineages following European colonization fit our empirical observations (model C, minimum probability of 0.995; table S6). As a further test, when the demographic model C was removed from the analysis, the principal components remained mostly unchanged (fig. S11). This observation supports the robustness of the PCMLR model in correctly discriminating between the different demographic alternatives (). However, none of the remaining demographic models (model A, B, D, E, F, or G) strongly fit all empirical observations (table S8), showing that model C is the only strong fit for the empirical data.

As a result, our ancient mitochondrial data suggest that European colonization was followed by local mass mortality and extinction of lineages associated with major population centers of the pre-Columbian past. Our results contrast with previous observations that Native American genetic diversity has been temporally and geographically stable for at least the past 2000 years (). However, the apparent contradiction between our study and earlier work is likely attributable to a significant improvement of sequence resolution. We also caution that the demographic models tested in this study are likely too simple to encompass the particulars of the colonization of the Americas—despite exploring scenarios with early emergence of geographic structure, as well as ongoing haplotype loss resulting from either serial founder effects and drift, or population displacement and elimination during the expansion of major Central Andean imperial states. Finally, more complete mitochondrial sequence data from early and modern-day populations will be needed to further refine the demographic models. In particular, it will be critical to ensure a strong spatial overlap between archaeological sites and present-day population locations (for example, Peru).

CONCLUSION

The additional resolving power gained from widespread temporal calibrations makes ancient mitochondrial DNA data from South America a very cost-effective approach to examining the precise timing, route, and sequence of events that led to the peopling of the Americas. In addition, this data set provides a genetic view of the possible role played by European colonization in reducing the overall Native American genetic diversity to the low levels observed today. In the absence of archaeological evidence in southern and eastern Beringia that is relevant to the peopling of the Americas, further resolution of the entry and routes of dispersal of the early settlers will be possible with ancient genomic data from South America.

MATERIALS AND METHODS

Archaeological samples

Most of the archaeological samples were collected in the Central Andes from sites in Peru (n = 70), western Bolivia (n = 9), and northern Chile (n = 6), whereas the remaining samples came from northern Mexico (n = 5) and the Argentinian Pampas (n = 2) (fig. S1 and table S2). Permissions from respective National Heritage organizations were obtained for all archaeological material analyzed in the study (section S1) and are available upon request. Radiocarbon dates were either published previously (, ) or obtained from the Oxford Radiocarbon Unit to confirm the archaeological assignation based on stratification and examination of artefacts (table S1).

Ancient DNA analysis

All ancient DNA work was conducted in clean-room facilities of the Australian Centre for Ancient DNA (ACAD) at the University of Adelaide, the Human Paleogenomics laboratory at the University of California at Santa Cruz (UCSC), and the Reich laboratory at Harvard Medical School (HMS). DNA contamination control measures are detailed in section S1. DNA was extracted using silica in solution (, ) and converted into double-stranded libraries with truncated Illumina adapters (). RNA or DNA baits were used for hybridization capture, followed by Illumina high-throughput sequencing (sections S1 to S4). Automated variant calls were confirmed by eye directly from the read pileups to avoid errors attributable to DNA damage (sections S2 and S4).

BEAST analysis

We used a Bayesian approach to estimate mitogenome coalescence times from the complete data using BEAST 1.8.0 () and reconstructed a demographic history via the extended Bayesian skyline method (), which infers effective population size. Three random subsets of 87 modern South American and Mexican mitogenomes were used with the ancient data in replicate analyses. To estimate the date of separation from Old World populations, we also added 20 of the most closely related east Siberian mitogenomes. The previously published 12.6-ky Anzick-1 genome from Montana, United States (), was included to provide genetic diversity for haplogroup D4h3a and to increase the depth of calibration points. We combined modern and dated ancient mitogenomes and analyzed three replicate data sets of 200 sequences. The 93 ancient samples (including 12 sequences older than 3.5 ky) ensured temporal (12.6 to 0.5 ka) and phylogenetic coverage (Fig. 3B), allowing internal calibration of the substitution rate for each South American haplogroup (fig. S9). Details are provided in section S5.

BayeSSC and model testing

We used BayeSSC () to simulate data for seven population scenarios (Fig. 4): one panmictic population without (model A) or with a post-European bottleneck (model B), or geographic separation (9 ka) of parallel populations featuring modern or ancient haplotypes (model D); geographic separation can be followed by either a post-European bottleneck for both populations (model E), a complete extinction of ancient lineages after contact with Europeans (model C), an early and progressive loss of lineages leading to extinction of ancient lineages (model G), or an early and progressive loss of ancient lineages followed by a rapid post-European extinction (model F). We calculated within-deme (Tajima’s D, haplotype diversity, number of haplotypes, etc.) and between-deme [Fixation index (FST), average pairwise distances] statistics for the observed and simulated data. We then performed a PCA to remove colinearity between the summary statistics, because they are likely to be correlated (, ), using the princomp function in R (). Finally, we fit a multiple logistic regression model () in which the response variable was the model under which the data were simulated and the explanatory variables were the first five principal components. This approach predicted which scenario was the most likely to have produced data like those that we observed. Details are provided in section S7.