Evolutionary population history of early Paleoamerican cranial morphology.

The nature and timing of the peopling of the Americas is a subject of intense debate. In particular, it is unclear whether high levels of between-group craniometric diversity in South America result from multiple migrations or from local diversification processes. Previous attempts to explain this diversity have largely focused on testing alternative dispersal or gene flow models, reaching conflicting or inconclusive results. Here, a novel analytical framework is applied to three-dimensional geometric morphometric data to partition the effects of population divergence from geographically mediated gene flow to understand the ancestry of the early South Americans in the context of global human history. The results show that Paleoamericans share a last common ancestor with contemporary Native American groups outside, rather than inside, the Americas. Therefore, and in accordance with some recent genomic studies, craniometric data suggest that the New World was populated by multiple waves of dispersion from northeast Asia throughout the late Pleistocene and early Holocene.

INTRODUCTION

South America was the last major continent to be colonized by modern humans (, ), yet it has unusually high among-population cranial differentiation relative to other global continents (–). This seems counterintuitive, given that within-group neutral genetic and craniometric diversity decreases with distance from sub-Saharan Africa, due to serial founder effects as humans dispersed out of Africa (–). However, populations can exhibit low within-group variation yet still show high between-group differentiation due to population isolation (reduced gene flow) and pervasive genetic drift, which is likely to be the case in South America (). High levels of among-population differentiation in the Americas have also been noted for linguistic () and neutral genetic data (). Whereas a concordant larger-than-expected cranial diversity is observable among late Holocene “Amerindian” populations (, , ), among-group differentiation is further exaggerated by the distinct cranial morphology of early “Paleoamerican” crania compared to the morphology of contemporary Native Americans (–), which has generated a long-standing debate about the origin of morphological diversity in the continent (, , , ).

Debates regarding the cause of this high between-group differentiation have centered on two main competing hypotheses. One possibility is that the observed diversity in South America is the result of in situ processes during the Holocene, whereby high within-group variation among early Americans became subdivided among descendent populations due to the rapid colonization of the Americas and/or as a result of genetic drift or natural selection acting in small isolated populations (, , –). Variants of this model emphasize the importance of recurrent gene flow between Asia and the Americas following the initial colonization of the continent (). However, whereas among-group cranial differentiation in South America is extraordinarily high, within-group variation for early Paleoamerican samples is not excessive and is within the range expressed by contemporary global populations (). This finding argues against the notion that the earliest migrants into the Americas were the source of all subsequent among-group biological diversity. The other main hypothesis proposed is that the observed cranial diversity is the result of multiple waves of dispersion into the Americas from northeast Asia over the course of several thousand years, with each wave of migrants introducing new sources of biological diversity. This argument is largely based on the empirical observation that the average cranial shape of the earliest South Americans bears stronger affinities with Australasian and Polynesian populations than it does with East Asian or later Native American groups (, , ). Recently, Hubbe et al. () suggested that early Paleoamerican groups retain the generalized ancestral morphology that characterized late Pleistocene Eurasian populations, as represented by fossils such as the Upper Cave specimen from Zhoukoudian (China) and Upper Paleolithic European specimens. If this ancestral morphology is also shared with contemporary Oceanic populations, then this would explain the apparent connection between Australasia and South America, despite their large geographic separation. Under such a scenario, subsequent population differentiation occurred in Asia following the initial settlement of the Americas, with later migrants into the New World resembling the derived “East Asian” morphology more closely.

Genomic data from early and contemporary Native Americans paint a complex picture of dispersal into and gene flow within the Americas (). Many studies support a single migration into subarctic America with subsequent population divergence and gene flow (, –), whereas others suggest up to four separate waves of migration from northeast Asia (–). However, there is widespread agreement that at least one separate wave of dispersal from outside the Americas is required to explain the genetic variation of paleoarctic and modern arctic groups (, –), although the current consensus suggests that all indigenous South Americans descend from only one wave of dispersion (, ). Although (paleo)genomic data may provide important context for the debates surrounding the origins of cranial diversity, minimal direct genomic data are currently available for Paleoamerican specimens in South America [see the study of Fehren-Schmitz et al. ()]. It should also be noted that, although all cranial populations have ancestors (that is, population history), not all fossil populations leave descendants, so Paleoamericans did not necessarily contribute to contemporary Native American genetic history. For these reasons, we take an explicit evolutionary population history approach based exclusively on the analysis of three-dimensional (3D) craniometric shape data, which have been repeatedly shown to provide an accurate proxy for neutral genetic expectation [see, for example, the study of von Cramon-Taubadel () for a recent review].

Previous attempts to explain the origins of South American cranial diversity have largely focused on testing alternative dispersal models (–, , ), with conflicting or inconclusive results. This approach is limited by the lack of a testable null hypothesis and the number of a priori assumptions required when constructing alternative dispersal models. Here, we overcome these limitations and use a novel methodological approach to investigate the population history of a well-known Paleoamerican series (Lagoa Santa, Brazil) in the context of contemporary global craniometric diversity. We use a multiple-effects apportionment model to disentangle two distinct geographically mediated evolutionary processes: (i) population divergence (“history”) as modeled using a bifurcating tree of hierarchical common ancestry and (ii) gene flow between contiguous populations as measured by between-group geographic distance (Fig. 1). When fossil specimens are discovered, there are only two guaranteed pieces of information available: their geographic location and their morphological properties. Hence, population (phylogenetic) history must be inferred via reference to other populations. Rather than construct arbitrary alternative models of population history and dispersal, we start with the conservative null hypothesis that the Paleoamericans share a most recent common ancestor (MRCA) with their nearest geographic neighbor (Chubut, Patagonia). If the subarctic American continent was populated by a single lineage arising in northeast Asia, then we would expect the null hypothesis of most recent shared common ancestry between Paleoamericans and contemporary South Americans to best explain the overall craniometric affinity patterns. Thereafter, using a comparative global craniometric database, we test every possible history model of most recent shared common ancestry for the Paleoamerican group. In each iteration, the effects of gene flow, as mediated by geographic distance, are kept constant. With a bootstrapping procedure, the null hypothesis is rejected if any alternative histories explain the patterns of among-group craniometric distances statistically better (P < 0.05). To verify that the findings were consistent across multiple cranial regions, we repeated the analyses for the 3D shape of the whole cranium, as well as the isolated shapes of the face, cranial vault, and the basicranium. We chose to focus on different aspects of cranial shape covariance that could be delineated on the basis of developmental and functional criteria () to account for the possibility that different subsets of cranial shape might reflect different population history signals ().

Fig. 1

Geographic and historical relationship between populations.

(A) Map showing the geographic position of each contemporary population and the Paleoamerican sample from Lagoa Santa (Brazil). Stars denote the waypoints used to calculate more realistic geographic distances between populations. (B) Tree topology representing the hierarchical model of historical divergence among populations. The null history model places the Paleoamericans as a sister group to their nearest geographic neighbor (Chubut, Patagonia).

RESULTS

Craniometric affinity analysis

Multidimensional scaling (MDS) plots of the craniometric distance matrices (Fig. 2) confirm the divergent cranial shape affinities of the Lagoa Santa Paleoamericans relative to other New World populations. In the case of the entire cranium (Fig. 2A), the Paleoamericans lie halfway between a cluster of populations from Africa, Australia, and the Andaman Islands and a cluster comprising New World and East Asian populations. Paleoamericans share affinity with Inuit populations on the second axis. A similar pattern is observed for the cranial vault (Fig. 2B). For the basicranium (Fig. 2C), the affinity between Paleoamericans and Inuit populations is most obvious with a clear separation from other New World and Asian populations on the first dimension. The affinity patterns for the face (Fig. 2D) are somewhat different, whereby the Paleoamericans are distinct from all New World and Asian populations on the second dimension. The MDS results confirm previous observations regarding the generalized affinities of the Lagoa Santa crania and their differences from East Asian and other Native American populations (, , , , , ). Mantel tests () confirmed that the overall among-population affinity patterns displayed by all four cranial data sets were significantly and positively correlated (P = 0.001).

Fig. 2

Population affinities based on craniometric data.

2D nonmetric MDS plots of population morphological affinities based on average Procrustes distance. (A) Cranium. (B) Vault. (C) Basicranium. (D) Face. Mantel tests were used to assess the congruence between affinity patterns based on different cranial data sets. The cranium was strongly correlated with the vault (r = 0.873), face (r = 0.786), and the basicranium (r = 0.704). The vault was also relatively strongly correlated with the face (r = 0.624) and the basicranium (r = 0.609), and the face and basicranium were moderately correlated (r = 0.327). All Mantel tests were significant at the α = 0.0083 level following Bonferroni adjustment.

Assessing alternative history models

Figure 3 and Table 1 summarize the results of the multiple-effects apportionment approach to assessing alternative history models for each of the four cranial data sets. In all cases, the null model suggesting that Paleoamericans share a MRCA with contemporary South Americans was rejected. The basicranium showed the overall weakest fit with the alternative models tested (Table 1), consistent with its overall low levels of among-population differentiation and conservative anatomy (–). The only alternative history model that produced a significant result for the basicranium was one of recent common ancestry between Paleoamericans and Greenland Inuit. In contrast, several alternative history models were found to be significant (P < 0.05) for the whole cranium, vault, and face (Fig. 3 and table S1). For the vault, the best-fit model suggested recent common ancestry between Paleoamericans and Inuit populations but not with the branch leading to subarctic American populations. The same was the case for the entire cranium, although it is worth noting that alternative models of recent common ancestry between Paleoamericans and Australians, Andamanese, and Mongolians were also statistically significant. The best-fit history model for the face data set suggested recent common ancestry between Paleoamericans and the ancestor of all New World and East Asian populations. Alternative models of recent common ancestry with Inuit groups and Australasians were also supported.

Fig. 3

Results of multiple-effects model comparisons.

History models found to be significantly better than the null model show the position of the Paleoamerican population (black diamond) in green (α = 0.05) and red (α = 0.01), with the absolute best model in black for each of the four cranial modules.

Table 1

Summary statistics for null, worst, and best history models tested for each cranial region.

Values represent the overall amount of among-group morphological distance explained (R2), with P values (α = 0.05) in parentheses.

	Worst model	Null model	Best model
Cranium	0.211 (0.841)	0.251 (0.272)	0.348 (<0.0001)
Vault	0.215 (0.931)	0.256 (0.396)	0.389 (<0.0001)
Face	0.291 (0.886)	0.341 (0.298)	0.438 (<0.0001)
Basicranium	0.155 (0.836)	0.205 (0.279)	0.265 (0.015)

DISCUSSION

The null hypothesis that Paleoamericans from Lagoa Santa share a MRCA with contemporary South Americans was rejected, irrespective of which cranial region was used. Moreover, alternative histories where Paleoamericans share a MRCA with subarctic Native Americans were also not statistically supported. In contrast, the best-fit history models suggest last common ancestry between Paleoamericans and other Native Americans outside the New World (Fig. 3). It is also worth noting that, although our results are consistent with the high levels of within-continent diversity noted previously, Lagoa Santa crania were not found to be outliers to contemporary modern human cranial variation. That is, their morphological variability falls within that observed among modern human populations, yet their overall morphology cannot be accounted for by a null hypothesis of shared common ancestry with all subarctic Native Americans.

The best-fit model consistently suggested by the cranium, vault, and basicranium data sets is one where Lagoa Santa Paleoamericans and New World populations share a MRCA in northeast Asia with a population that later gave rise to the populations that colonized the American arctic and Greenland. Moreover, the face data set suggests a more ancient common ancestral link between Paleoamericans and the ancestors of all East Asians and New World populations. This is not inconsistent with the results from the vault and basicranium but does suggest that Paleoamerican facial shape is more generalized and plesiomorphic than vault or basicranial shape. This is also underlined by the fact that some significant alternative histories for the facial data set suggest common ancestry between Paleoamericans and Australasians, which would suggest the retention of a more plesiomorphic, generalized facial shape from the common ancestor of all Australasians and New World populations.

Morphological similarities between Greenland Inuit and Siberian populations have been noted previously (, , ) and have largely been interpreted as evidence for circumarctic gene flow across the Bering Strait. Galland and Friess () suggest that the observed affinity may be due to cold adaptation in these circumarctic populations. Although it is entirely plausible that ancestral Siberian cranial shape was subject to natural selection in response to cold stress (, ) with continuing selective pressure in descendant arctic New World populations, cold adaptation is unlikely to explain Paleoamerican morphology, given the warm climates encountered in eastern South America. It is also possible that the results for the face reflect some similarities in masticatory stress () shared between the Paleoamericans and other hunter-gatherer populations in the data set. However, a functional explanation cannot fully account for the facial results, given the strong population history signal found in facial shape data (see Fig. 2D) () and the fact that the Lagoa Santa Paleoamericans do not share a MRCA with the Chubut, who are also a foraging population. Therefore, our results support a model of most recent shared ancestry between Siberian, Inuit, and Paleoamerican populations as a more parsimonious explanation of the data. This explanation is also consistent with paleogenomic studies that demonstrate shared ancestry between Native Americans and Upper Paleolithic Siberians () that is not found in contemporary East Asians. This suggests that at least some Native American lineages diverged from ancestral Siberians and entered the Americas before the diversification of East Asian populations.

Our results are also in accordance with the recent genomic link found between Amazonian populations and Australasians (, ). Skoglund et al. () suggest than an ancient Native American lineage, named “Population Y,” could have resulted from a highly substructured ancestral northeast Asian population that shared strong genetic affinities with the ancestors of modern Australasians. There is mounting genetic and morphological evidence for at least two major waves of dispersal into Asia from Africa, with Australomelanesians representing modern descendants of the earlier migration (, ). There is also genomic evidence that northeast Asia was continually occupied throughout the Last Glacial Maximum (~21 thousand years ago) (). These spatiotemporal dynamics would have provided ample opportunity for population substructure to emerge within Siberia and Beringia and thus provide several distinct sources of Native American ancestry from the same geographic location through time. Earlier (Paleo)siberian populations would have shared greater genetic affinity with Australasians further south as an outcome of their shared out-of-Africa dispersal history. However, as time progressed, further dispersal from Africa along with differentiation and gene flow within Asia would have altered the genetic signature of the northeast Asian source populations that gave rise to later Paleoeskimo and (possibly) other Native American populations. If such a spatiotemporal scenario is borne out by future genomic studies of northeast Asians and Native Americans, then it would solve the apparent conundrum of how a single “source” population in Siberia/Beringia could have given rise to several distinct genetic and phenotypic lineages within the Americas over time.

CONCLUSION

The nature and timing of the peopling of the Americas is a subject of continuing debate. It is now clear that people entered the Americas from northeast Asia via Beringia by at least 15,000 years ago (, ), and dispersal into South America proceeded quickly, most likely following a coastal Pacific route (, ). A genetic distinction between western and eastern South American populations has been noted (, ), particularly with respect to a signal of Australomelanesian ancestry in eastern (Amazonian) populations (, ). This is consistent with our cranial findings, which suggest a layered population history for South America, with at least two major sources of biological variation from Asia. The earliest (Paleoamerican) migrants were morphologically distinct from later groups, although structured gene flow among the descendants of Paleoamericans and later populations may have contributed to their assimilation in the late Holocene. Finally, our results underline the potential of cranial shape data for drawing powerful inferences regarding past population history when analyzed in an explicit model-bound microevolutionary context. There are still many questions in human prehistory for which relevant genomic data are not, or may never be, available. Therefore, we advocate that neutrally evolving morphological data be viewed as a useful adjunct to, and not an inferior alternative to, genomic data for analyzing our evolutionary history.

MATERIALS AND METHODS

Experimental design

The objective of this study was to test whether there were any alternative population history models for the Lagoa Santa Paleoamericans that provided a statistically better fit to the craniometric data than the null model, which assumed a MRCA linking the Paleoamericans and South American Chubut. The experimental design was as follows: For the same 17 extant human populations, pairwise population distance matrices reflecting “history” (that is, the phylogenetic pattern of shared common ancestry), “geography” (that is, the geographic distances between groups), and craniometric affinity were constructed. The proportion of craniometric distance explained by history, geography, and the interaction between history and geography was partitioned using a multiple-effects model. The null hypothesis for the phylogenetic position of the Paleoamericans was that they share a recent common ancestor with the South American Chubut. The entire craniometric data (including the Paleoamericans) were bootstrapped 1000 times by resampling (with substitution) the individuals within each series, keeping the original sample sizes constant, to generate 1000 craniometric distance matrices that could be compared with the null history model. Thereafter, each alternative history model (that is, where the Paleoamericans were placed as a sister taxon to every terminal branch and every internal node) was tested against the initial craniometric distance matrix. In each iteration of the multiple-effects model, geography was kept constant, so any increased fit of the data to the model was entirely due to the inputted history model. Any alternative history model that explained the craniometric distances better than the 95% confidence interval generated from the bootstrapped analyses was deemed statistically significant.

Craniometric data

The craniometric data set comprised 3D landmark configurations for samples of 17 global human populations and a sample of Paleoamerican specimens from the Lagoa Santa region of Brazil (Table 2 and Fig. 1). All specimens were adult (with fully fused sphenooccipital synchondroses) and were sexed using standard osteological techniques (). The total configuration comprised 135 landmarks and was subdivided to quantify three functional-developmental cranial modules: the vault, face, and basicranium (). Full anatomical descriptions of all landmarks can be found in table S2. All landmark data were collected by a single observer (N.v.C.-T.) using a MicroScribe 3DX digitizer to avoid interobserver bias. Intraobserver error was tested for all landmarks () and was found to be <1 mm using the partial superimposition method ().

Table 2

Human population craniometric samples used.

Sample sizes for Paleoamericans are 18 for whole cranium, 45 for vault, 32 for basicranium, and 27 for face. NHM, Natural History Museum (London, U.K.); MH, Museé de l’Homme (Paris, France); AMNH, American Museum of Natural History (New York, NY, USA); NHMW, Das Naturhistorische Museum, Wien (Vienna, Austria); DC, Duckworth Collections (Cambridge, U.K.); SNMNH, Smithsonian National Museum of Natural History (Washington, DC, USA); MLP, Museo de La Plata (La Plata, Argentina); ZMD, Zoological Museum, University of Copenhagen (Copenhagen, Denmark); RIO, Federal University of Rio National Museum (Rio de Janeiro, Brazil); BH, Museu de História Natural, Federal University of Minas Gerais (Belo Horizonte, Brazil); USP, University of São Paulo (São Paulo, Brazil).

Population	n	Latitude, longitude	Museum
San	31	−21.0, 20.0	NHM, MH, AMNH, NHMW, DC
Biaka	21	4.0, 17.0	NHM, MH
Ibo	30	7.5, 5.0	NHM
Zulu	30	−28.0, 31.0	NHM
Berber	30	32.0, 3.0	MH
Italian	30	46.0, 10.0	NHMW
Basque	30	43.0, 0.0	MH
Russian	30	61.0, 40.0	NHMW
Australian	30	−22.0, 126.0	DC
Andaman	28	12.4, 92.8	NHM
Mongolian	30	45.0, 111.0	MH
Chinese	30	32.5, 114.0	NHMW
Japanese	30	38.0, 138.0	MH
Alaskan	30	69.0, −158.0	AMNH
Greenland	30	70.5, −53.0	SNMNH
Hawikuh	30	33.5, −109.0	SNMNH
Chubut	30	−43.7, −68.7	MLP
Paleoamerican		−19.4, −44.0	ZMD, RIO, BH, USP

Missing data estimation

The Lagoa Santa crania were in variable states of completeness and preservation. In all cases, specimens were required to have at least 70% of the landmarks digitized to be included in an analysis. Missing landmark positions were estimated using reflected relabeling for bilateral points and thin-plate spline interpolation () using the entire contemporary reference data set in R 3.1.3 (estimate.missing function from the geomorph package).

Geometric morphometrics

Each cranial module landmark configuration was subject to generalized Procrustes analysis and tangent space projection in MorphoJ 1.06 to remove variation among specimens due to isometric scaling, rotation, and translation. The resultant Procrustes variables were used to calculate the average pairwise Procrustes distances among all 18 population samples (resulting in four craniometric Procrustes distance matrices). Average between-population craniometric distance patterns were visually inspected via 2D nonmetric MDS analysis ().

Geographic distances

A matrix of between-population great-circle geographic distances (in kilometers) was calculated using the Geographic Distance Matrix Generator version 1.2.3 () based on the geographic coordinates shown in Table 2. The following waypoints were used to calculate more realistic estimates of between-population geographic distances: Cairo, Egypt (30.0, 31.0); Istanbul, Turkey (41.0, 28.0); Phnom Penh, Cambodia (11.0, 104.0); Anadyr, Russia (64.0, 177.0); Prince Rupert, Canada (54.0, −130.0); and Panama City, Panama (9.1, −79.4) (shown as stars in Fig. 1).

Constructing the history matrix

A distance matrix reflecting the nested hierarchical pattern of population divergence (history) was constructed on the basis of the branching topology of the consensus neutral genetic tree of population relatedness (, ). Most of the cranial populations were represented in the data sets analyzed by Pemberton et al. (), so we used their consensus neighbor-joining tree topology based on 246 neutral microsatellite loci as our main framework for the topology of the tree (Fig. 1). The relative branching relationship between the four New World populations (Alaskan, Greenland, Hawikuh, and Chubut) was inferred from the neighbor-joining analysis of over 360,000 single-nucleotide polymorphisms presented by Reich et al. (), who showed that the two Inuit groups share a MRCA [see also the study of Raghavan et al. ()] relative to the Chubut (Patagonians) and the Hawikuh (Central Amerinds). Australians and Andamanese were inferred to share a more recent common ancestor relative to all East Asian populations based on the findings of Rasmussen et al. () and Reyes-Centeno et al. (). The resultant history matrix reflects the hierarchy of relatedness among all populations and was coded such that each node represented one unit. Therefore, the distance between two populations sharing a single node was quantified as 1, whereas two populations separated by six nodes were given a value of 6, etc. We subjected the resultant history matrix to a neighbor-joining analysis () to confirm that the history distance matrix reflected the branching relationship shown in Fig. 1.

Statistical analysis: Multiple-effects model

We adopted a modified version of the approach taken by de Campos Telles and Diniz-Filho (), designed to separate out the unique contribution of two different geographically mediated processes (historical divergence and gene flow) on global human cranial affinity patterns. As shown in Fig. 1, the history model inferred from neutral genetic data is largely geographically structured, because it reflects the dispersal of modern humans out of Africa and the subsequent colonization of major continental regions (). Hence, the history matrix reflects the historical bifurcations among populations as a nested series of hierarchically arranged common ancestors. However, geographic distance has also mediated historical gene flow patterns, because contiguous populations are more likely to experience bilateral gene flow over time (). Following de Campos Telles and Diniz-Filho (), the multiple-effects model used here works by apportioning the variation in among-population craniometric (Procrustes) distance due to (i) pure history, (ii) pure geography, (iii) interaction between history and geography, and (iv) residual variance. We manipulated the multiple-effects model to keep the geographic distance constant but altered the phylogenetic position of the Lagoa Santa sample within the history matrix with each iteration. Therefore, alternative history models that explain a larger proportion of the among-group Procrustes distances would reduce the amount of unexplained variance in the overall model. Our null hypothesis was that Paleoamericans (Lagoa Santa) are the phylogenetic sister group of their nearest geographic neighbor (Chubut, Patagonia). Thereafter, we reran the multiple-effects model, with Paleoamericans iteratively placed as a sister taxon to every possible terminal branch and internal node in the tree shown in Fig. 1. To test whether any alternative history model explained a significantly greater proportion of the overall variance than the null model, we used a bootstrapping procedure. For each of the four cranial modules, we created 1000 bootstrap replicates of the Procrustes coordinate data set by replacing, with substitution, the individuals within each population series, and then we ran the 1000 resultant Procrustes distance matrices through the multiple-effects model using the null hypothetical history model (where Paleoamericans are the sister taxon to Chubut) to generate null distributions of model-fit expectation. Thereafter, any alternative history model that explained more overall craniometric distance variance than the 95% confidence interval for the appropriate null distribution was deemed significantly better than the null history model (one-tailed P < 0.05). Procedures for performing data bootstrapping and the multiple-effects model were undertaken in R. 3.3.0 (R Core Team, 2016) using code written by M.H. and complemented by functions from packages phytools (), vegan (), and MASS (). All data matrices used in the analyses are available from the lead author upon request.