The origin of the Mexicas has been the focus of substantial debate among archaeologists, linguists, geneticists and physical anthropologists in the past. After the downfall of the Toltec Empire (AD 900–1170), the population occupying the capital city of Tula dispersed, enabling immigrants to occupy Central Mexico (Townsend 1992). The classical foundation tale tells that the first Mexicas would have started a pilgrimage from their homeland, the mythological city of Aztlán (the heron's city, located somewhere in the north) towards their Promised Land. Lead by their hero-god Huitzilopochtli, they finally reached the Central Valley and, in 1325, Tenochtitlán, the central city of the empire, was founded on one of the islands of Lake Texcoco (Townsend 1992). His sister city, Tlatelolco, was founded some years later, around 1337 (Matos Moctezuma 1989). This story, as revealed by the Boturini Codex, as well as many archaeological steles, constitutes the modern Mexican legend of the Central Valley people (Orozco y Berra 1944; Corona Núñez 1968; Davies 1973). As recently noted by Kemp et al. (2005), it is uncertain whether this is a true reflection of their history or a fictionalized account, for the Toltecs also claimed an origin to the north. Since the Mexicas valued their ancestors to a great extent, copying them stylistically and architecturally, the Mexica foundation tale might in fact be an adoption and modification of the Toltec story (Townsend 1992; Kemp et al. 2005).

The Mexicas have been extensively studied by linguists, since they belong to the nahua sub-family which is nested into the Uto-Aztecan, one of the most diverse and geographically dispersed linguistic families (Campbell 1997). Uto-Aztecan speaking populations are distributed from the American Southwest to Panama (Miller 1983; Hill 2001). Taking into account the biogeography and reconstruction of Proto-Uto-Aztecan words, many linguists support the Northern Origin model. This model states that Mesoamerican nahua speakers, like the Mexicas, are the final result of a migration from the north (Dyen 1956; Fowler 1983; Miller 1984). However, Hill (2001) suggested that the Proto-Uto-Aztecan speech community was located in Mesoamerica and spread northward into the current range forced by the demographic pressure associated with maize cultivation. Finally, some authors indicate that the source population may well be in Mesoamerica, but not strictly in the Central Valley. Thus, the regions of Michoacán, Guanajuato, Nayarit or even the Gulf of Mexico coast, near Veracruz, should also be considered (Townsend 1992).

But what was the demographic landscape prior to the origin of the Mexica Empire? The archaeological evidence indicates that at the Mexicas' hypothetical arrival, the main settlement prevailing in Mexico Central Valley was Azcapotzalco (Matos Moctezuma 1989; Townsend 1992). This city was ruled by the Tepanecos, and constituted the main economic, military and cultural centre since the Early Postclassic, four centuries earlier. During the long warfare that predated the definitive organization of the Triple Alliance (Tlacopán, Tenochtitlán and Texcoco), Mexicas would have banished the Tepanecos a century after the foundation of the Great Tenochtitlán (Matos Moctezuma 1989; Townsend 1992). Again, it remains uncertain whether the initial belligerent scenario promoted some population exchange and admixture among inhabitants of the basin, or else was just a political conflict without demographic consequences.

Archaeological, linguistic and historical evidence must be complemented with biological data in order to achieve an integrative interpretation of the settlement process. In a recent work, Kemp et al. (2005) analysed the mtDNA of Postclassic Mexicas from Tlatelolco. The haplogroup frequencies exhibited by the Mexicas resemble those of other present-day central and southern Mexican, as well as Central American populations, but differ from American Southwest samples. These results do not support Hill's (2001) hypothesis of a Central Mexican origin of Uto-Aztecan. Instead, they suggest that populations from Mesoamerica had little maternal influence on the genetic structure of groups residing in the American Southwest. However, if the hypothetical expansion of the Proto-Uto-Aztecan from Central Mexico encompassed some degree of genetic homogenization by admixture with northern groups, or else was exclusively a case of cultural diffusion, remains unsolved.

Most issues concerning the Mexica settlement of the Central Valley cannot be directly addressed by archaeology, linguistics or history, since cultural shifts are not always congruent with population replacement, extinction, growth, migration, etc. In addition, mtDNA analyses are usually based on interpretation of past processes departing from the analysis of modern molecular variation. Nevertheless, in some densely populated, historically and demographically complex areas like Mesoamerica, the reconstruction of past events departing from modern groups cannot be accurate enough to identify continuity or replacement across a long period of time.

Diachronic analyses performed to detect biological discontinuities among different periods are of crucial importance to achieve a thorough comprehension of the biological dimension of the numerous cultural transitions that occurred in Mesoamerica. To solve this, analyses of craniofacial size and shape variation across diachronic sequences are extensively used to reconstruct the structure and history of extinct populations. Although these markers are to some extent affected by non-genetic factors, they present some advantages over other kinds of data when dealing with reconstruction of historical processes. For instance, when samples from different regions and periods are available, it is possible to apply different models specifically designed to test for the action of particular microevolutionary agents (Relethford & Blangero 1990) or the existence of discontinuities across diachronic (Konigsberg 1990; Relethford 1991; Steadman 2001; Martínez-Abadías et al. 2005; Stojanowski 2005; Brace et al. 2006) or geographical sequences (González-José et al. 2002, 2003; Brace et al. 2006). In addition, these kind of tests are necessary to understand to what extent great scale technological and cultural transitions are correlated with changes in the gene pool of populations under scrutiny (Brace et al. 2006).

Here, we analyse Mexicas' facial morphology in relation to Preclassic, Classic, Early Postclassic, Colonial groups from the Central Valley and Mexican northern groups in order to asses whether the Late Postclassic facial pattern is compatible with either a scenario of biological continuity in the Central Valley or with population replacement led by northern migrants. More specifically, we test whether facial traits observed in Mexica samples were already present in the Central Valley of Mexico during the Early Postclassic and Classic periods or, on the contrary, are distinctive features of northern groups.

The study includes 331 adult individuals from 11 archaeologically different populations covering all main periods of Mexico Central Valley, as well as main regions from Northern Mexico and one Early Colonial, post-contact group. All skulls are stored at the Dirección de Antropología Física in the Museo Nacional de Antropología, México. The sample composition is provided in table 1. Great care was taken for specimen selection to avoid sample bias, but the availability of material in museum collections did not allow us to gather data from age and sex-matched samples. Three-dimensional landmark coordinates were collected using a Microscribe G2X digitizer. Each individual is represented by 13 homologous landmarks in three dimensions: prosthion; right and left ectomolare; right and left nariale; right and left alare; nasion; right and left zygoorbitale; right and left frontomalare-orbitale; and posterior nasal spine (see figure 1 in the electronic supplementary material). Since artificial cranial deformation is ubiquitous in Mesoamerican cultures, and in some populations it is present in almost all individuals (e.g. Teotihuacan and Tlatilco), we were compelled to include deformed individuals taking care to use a landmark configuration not affected by the effects of deformation. To explore such effects, we performed several preliminary analyses considering just non-deformed skulls and a landmark configuration covering the cranial face, base and vault (prosthion, ectomolare, nariale, alare, nasion, asterion, zygoorbitale, fronto-malare anterior, stephanion, glabella, bregma, lambda, inion, basion, hormion, posterior nasal spine, jugale, superior point on the zygotemporal suture and inferior point on the zygotemporal suture). Since the subset of 13 facial landmarks taken in all deformed and non-deformed specimens provided the same pattern of distances as the overall cranial configuration, we opted for the facial configuration in order to include large-sized samples (despite the high proportion of deformed specimens).

Table 1 Description of the data used in the present study. F, females; M, males.

Landmark configurations were superimposed by generalized Procrustes analysis (Rohlf & Slice 1990; Goodall 1991) using the Morphologika software (O'Higgins & Jones 2004, freely available from http://www.york.ac.uk/res/fme/resources/software.htm). This procedure translates skull configurations to a common origin, scales them to unit centroid size and rotates them to best fit by using a least-squares criterion. The fitted coordinate configurations resulting from these procedures are then placed in the denominated Kendall's shape space (Rohlf 1996). As this shape space is non-Euclidean, further statistical analyses were performed by projecting the coordinates into a linear tangent space (Dryden & Mardia 1998).

Data used in this study were collected by two of us (R.G.-J. and N.M.-A.) using the same landmark collection protocol. The effects of intra- and interobserver error were evaluated by collecting 10 observations of a single specimen, a male African cranium. After the Procrustes superimposition, the Euclidean distance of each landmark to its respective centroid was computed and, for each observer, landmark deviations were calculated relative to the observer landmark mean. Following Singleton (2002), mean deviations were calculated for individual landmarks and subsequently averaged to give a mean deviation for each observer across all landmarks. One-way analysis of variance was performed for each landmark by observer and the root mean squares were evaluated. The root of the within-groups mean squares (root mean square error) is an estimation of the intraobserver error (Sokal & Rohlf 1995), while the root of between-groups mean squares denotes interobserver error. Mean deviation was 0.275 and 0.370mm for N.M.-A. and R.G.-J., respectively, and interobserver error averages 1.22mm. Considering the relatively large size of the faces studied here, these margins of error were considered acceptable.

To model the magnitude and direction of gene flow across regions and periods, shape variation was analysed using three independent methods. First, a canonical variates analysis (CVA) was performed on the Procrustes coordinates, in order to maximize the separation between groups. Shape variance patterns present in the Procrustes shape space coordinates are distorted when the shape space is stretched and compressed by techniques that, like CVA, attempt to maximize group discrimination. Recently, Klingenberg & Monteiro (2005) remarked that distinctness of groups depends not only on the distance between sample means, but also on the direction of the mean difference relative to the directionality of variation within groups (Klingenberg & Monteiro 2005, p. 680). Since in most cases variation is non-isotropic, that is landmarks differ in their amounts and direction of variation, the CVA was implemented by scaling the multidimensional space by the inverse of the within-group variation of the Procrustes coordinates. The importance of such transformation has already been discussed by Albrecht (1980), Carroll et al. (1997), Klingenberg (2003) and Klingenberg et al. (2003,2004). Further details of the calculation can be found in Klingenberg & Monteiro (2005). The distances are obtained after the transformed data measure between-group differences relative to the within-group variation. Thus, the matrix of Mahalanobis distances was obtained after transformation of data. In addition, sex effects were tested by plotting the frequency histograms of both sexes along the first canonical variates (results not shown). Since male and female distributions were fully overlapped (Mahalanobis distance among sexes was not significant), we assume that inter-group differences are much more important than sexual dimorphism. We used a pre-release version of the MorphoJ software (C. P. Klingenberg 2006, unpublished) to perform the transformation and to obtain the Mahalanobis distances. The resulting pattern of distances was visualized by a Neighbour joining tree (Saitou & Nei 1987).

Second, detection of differential gene flow was estimated following the model proposed by Harpending & Ward (1982) and extended to quantitative traits by Relethford & Blangero (1990). To do so, the fitted coordinate configurations of the specimens were analysed using principal components (PCs) analysis, in order to explore how variation is partitioned within and among the samples, as well as to achieve data reduction (Vidarsdottir et al. 2002). PC scores were standardized using a within-sex z-score transformation (Williams-Blangero & Blangero 1989) and the sexes were pooled for all subsequent analyses. Visualization of shape differences along the PCs was obtained by warping the triangulated surface of the mean shape to represent shapes at any position within the plot, using the loadings of original landmark coordinates on these PCs (O'Higgins & Jones 1998; O'Higgins & Vidarsdottir 1999). PC scores were obtained using the Morphologika package (O'Higgins & Jones 2004). The Relethford & Blangero (1990) model states that, under certain assumptions, observed and expected genetic or phenotypic variances can be compared to assess the effects of differential external gene flow. Thus, greater than average gene flow from an external source will result in a population with greater within-group phenotypic variance than expected. Conversely, less than average gene flow from an external source will result in a population with less within-group phenotypic variance than expected (Relethford & Blangero 1990). Estimations of minimum genetic differentiation (Fst), observed, expected and residual within-group phenotypic variance, were obtained after an R matrix computed on the PC's scores. Data analysis and parameter computations were performed using the RMet for Windows v. 5.0 software (Relethford 2004).

Finally, we carried out an analysis of biological, geographical and chronological distance correlation (Konigsberg 1990). According to this model of geographical and temporal isolation, if little gene flow occurred between regions then there will be a significant positive correlation between geographical and biological distances when temporal distance is controlled (Konigsberg 1990). The geographical distance matrix resumes the distances in kilometres among archaeological sites from where skull samples were drawn; chronological distances were estimated as absolute differences in the dating of sites in years before present and biological distances correspond to the Mahalanobis distances matrix obtained above.

The logical basis for the combination of the three procedures is grounded on the underlying principles of population genetics. For instance, if the origin of the Mexica population during the Late Postclassic is the final result of a demographic replacement, greater distances between TLA and the remaining, pre-Late Postclassic Central basin groups (AZC, TUL, TEO and TCO; see table 1 for abbreviation definitions) are expected. In parallel, significant population movement into the Central Valley at the transition between Early and Late Postclassic will be detected as perturbations in the local population structure, since new alleles were introduced after the putative replacement. Finally, if the migrants were coming from somewhere in the north, then biological distances between north and Mexica populations will be smaller than among temporally earlier populations and their within-group heterogeneity would then increase.

Inter-population shape differences were examined using Mahalanobis' distances computed from the Procrustes-aligned coordinates. The Neighbour joining tree obtained after the Mahalanobis' distances is presented in figure 1. Previous analyses made on a sample of exclusively non-deformed skulls studied after facial, basal and vault landmarks resulted in a tree showing the same topology (results not shown). Thus, choice of facial landmarks enables using greater sample sizes, an important detail when dealing with estimations of internal variance (see below). In addition, it supports the robustness of the Neighbour-joining algorithm as a representation of the sample's affinities in this particular dataset. Figure 1 also shows facial shape variations as deformation polygons of each population (targets) onto the consensus (reference). On one hand, Mexicas from Tlatelolco are connected to the Early Postclassic group of Azcapotzalco, as well as to their putative descendants from the early Colonial, and to the Tarahumaras from northeastern Mexico (further details on the chronology of samples and codes are presented in table 1). On the other hand, Classic, Preclassic and northern groups are placed in the remaining branch of the tree showing no geographical or chronological structure. When compared to this non-Postclassic cluster, TLA, AZC, SOL and TAR are characterized by short faces due to a relative reduction in length of the posterior face. Furthermore, these groups present an expansion in width of the posterior face. This change in such a complex morphological module is evidenced by the distance between the right-to-left ectomolars' line and the posterior nasal spine, and by the width of the palate, measured as the distance among right and left ectomolars. In addition, the samples belonging to this cluster present relatively higher orbits and nasal apertures, and relatively low subnasal prognathism. In general terms, all these characters contribute to a greater facial flattening when compared with northerners and Classic groups. Considering the nasal aperture, there is also a distinctive feature characterizing Postclassic, Colonials and Tarahumaras: nasal aperture tends to present a wider basis and the maximum width is found in a lower position in the aperture. The subcluster including Baja California Sur, La Candelaria, Paila and Tlatilco is also defined by an oblique position of the orbits, supported by the more sagittal position of the orbitale.

Figure 1 Neighbour joining tree based on the Mahalanobis distances obtained after the PC scores. Polygons are lateral and frontal visual displays of group differences: the thick, solid lines depict the Procrustes mean of each group and the thin lines represent (more ...)

Even though changes involved are quite subtle, facial morphology evidences a morphological discontinuity across the chronological sequence of the Central Valley. This disruption seems to occur at the Classic–Postclassic transition, around AD 900, rather than during the Premexica–Mexica transition (AD 1200). Continuity across the Classic and the Epiclassic is well supported by the similarity observed among samples from Teotihuacan and Tula. This result is also congruent with archaeological hypothesis pointing to a military, economic and demographic influence of Teotihuacan over Tula at the end of the Classic (Matos Moctezuma 1989; Rattray Childs 2001).

Remarkably, facial attributes defining the whole Preclassic, Classic, Epiclassic and a wide spectrum of northern populations are also found in the Pericú from Baja California, which were recently reported as showing the most resembling phenotypic pattern with Paleoamericans, the earliest inhabitants of the New World (González-José et al. 2003). In other words, the facial pattern present before the rise of the Postclassic world could be widely defined as ‘ancestral’ in terms of New World facial variability, and the Mexica one as ‘derived’, at least in the context of Mesoamerican facial variation.

A first R matrix analysis including all samples shows that the minimum genetic differentiation (Fst) among all populations is 0.052. This value is less than half the minimum genetic differentiation for diachronic populations within the Central Valley, that is, the computation made on the sub-sample including TCO, TEO, TUL, AZC, TLA (Fst=0.125). Furthermore, it is similar to the Fst values obtained considering just the non-Central Valley groups (CAN, PAI, BCS, SON, TAR; Fst=0.039; see table 1 for abbreviation definitions). This indicates that there is greater variation among periods within the Central Valley than between regions when all samples are considered.

The Relethford–Blangero analysis presented in table 2 shows that both Postclassic groups have greater than average external gene flow, but the Azcapotzalco group received significantly more gene flow. Given the relatively large biological distances between Postclassic and Classic plus northern groups (figure 1), it is unlikely that Tepanecos from Azcapotzalco and Mexicas were receiving large amounts of migrants from their Central Valley's predecessors or populations from the far north. Furthermore, the negative residuals in most of the Classic and northern groups suggest that migration in these periods and regions was more restricted than gene flow after the Postclassic period. Thus, our results indicate relatively little gene flow among Classic and Postclassic groups as well as modest interregional contact among the north and Tepanecos or Mexicas.

Table 2 Minimum genetic distances from the centroid (rii), observed mean variances, expected mean variances and residual variances estimated using a heritability of h2=1, and computed according to Relethford & Blangero (1990). (The lower part of table (more ...)

To help clarify the patterns and magnitude of differential gene flow, the Relethford–Blangero analysis was repeated, grouping the Central Valley populations into two main cultural periods: Preclassic plus Classic (TCO, TEO, TUL) and Postclassic (AZC, TLA). Again, the Postclassic sample shows a great deal of heterogeneity, contrasted with the almost average differentiation of the Classic group (table 2).

Konigsberg's (1990) model predicts that if temporal separation is mathematically controlled, then biological differentiation must be associated with spatial separation in the case of limited or null gene flow among regions. The matrix correlation analysis performed to evaluate this model shows that biological and geographical distances (with the temporal distance matrix held constant) were not significantly correlated (partial correlation coefficient r=0.390, p=0.085). This suggests that distance among the Central Valley and the north (via the Sierra Madre inhabited by the Tarahumaras) was not an effective barrier to some level of gene flow. However, a stronger correlation is expected since most populations come from a single delimited area, the Central Valley. This was not the case and, as observed in figure 1, a main cause for this disruption is the divergence among Classic and Postclassic groups: whereas both populations occupied the same region, they presented the higher levels of facial divergence. This result is also compatible with the detection of new polymorphisms coming from an external source suggested by the Relethford–Blangero analysis (table 2).

A possible influence of an external population source during this transition is supported by the Relethford–Blangero analysis and the partial correlation test. However, and considering the Mahalanobis distances, a migration-replacement led by northern groups is unlikely to be responsible for this abrupt change. If so, northern groups must have resembled the facial morphology of Postclassic samples and these groups must have presented average internal variation. This was not the case observed here. Results suggest that some population not considered in our study might be responsible for the origin not just of the Mexicas, but of the Postclassic Central Valley's societies in general. Although we have no direct evidence to recognize the ancestors of the first Postclassic groups (e.g. Tepanecos from Azcapotzalco), most likely candidates can be found in the groups who achieved some degree of demographic, socio-economic and political stability during and after the fall of the Classic centres, as occurred in Nayarit, Jalisco, Guerrero and Michoacán. In the case of Michoacán, settlement and land-use patterns are well documented on the Lake Pátzcuaro Basin archaeological record and paleoenvironmental reconstructions (Fisher et al. 2003). During the Classic, for instance, ceremonial centres appear throughout the region (Pollard 1997; Pollard & Cahue 1999). This pattern continues through the Early/Middle Postclassic with increased population growth. According to Fisher et al. (2003), during the Late Postclassic the Basin served as the geopolitical core of the Tarascan Empire, with centralized economic and political systems, and represented the largest population of the Prehispanic period (Pollard & Gorenstein 1980; Fisher et al. 2003; Pollard 2003).

Thus, even though the collapse of the Classic centres was a vast and ubiquitous phenomenon, some places such as the Lake Pátzcuaro Basin in Michoacán possibly achieved positive population growth rates that enabled migration to empty or poorly settled areas in the Central Valley.

Linguistic, genetic, archaeological and historic–ethnographic hypotheses compete to explain the biological origin of Mexicas. Some linguists support a Uto-Aztecan unity spreading from Mesoamerica to the American Southwest (Hill 2001), while some geneticists (Kemp et al. 2005) discard a regular population contact among Mexicas and northerners. However, we must note that there is no contradiction between both statements, since the expansion of the Uto-Aztecan, if linked to the expansion of crop production, would have been a fairly earlier event. Furthermore, considering that the cultural practices that led to Zea domestication probably occurred as early as 6250 years ago (Piperno & Flannery 2001) and that the complete process of artificial selection was not completed until 2000 years ago (Jaenicke-Despres et al. 2003), it is very likely that the putative expansion of the Proto-Uto-Aztecan was in fact led by Mesoamerican groups much before the raise of the Mexicas.

In this sense, our work is in agreement with the results of Kemp et al. (2005). Both studies disregard intense contact among Mexicas and northerners. However, our analyses identify that the stronger biological disruption occurred previous to the formation of the Mexica Empire. Hence, they support the idea that the expansion of the Uto-Aztecan was a phenomenon that largely antecedes the Postclassic period.

In conclusion, our results show that Classic and northern populations have had little or null effect over the facial shape and maybe the gene pool of Mexicas. Conversely, the most remarkable phenotypic disruption along the chronological sequence of the Central Valley is observed in the Classic–Postclassic transition. This phenotypic disruption is in agreement with the idea that the Classic collapse encompassed important demographic transformations affecting the biology of human populations. It is also in accordance with the paleoenvironmental evidence reflecting climatic changes at this epoch. Tree ring and lake sediment records indicate that some of the most severe and prolonged droughts to impact North America and Mesoamerica in the past 1000–4000 years occurred between AD 650 and 1000, particularly during the eighth and ninth centuries, a period of time that coincides with the Terminal Classic Period (Acuña-Soto et al. 2005). Based on the similarities of the climatic (severe drought) and demographic (massive population loss) events in Mesoamerica during the sixteenth century, drought-associated epidemics of hemorrhagic fever may have contributed to the massive population loss during the end of the Classic period (Acuña-Soto et al. 2005).

The fact that the Mexica shape pattern was already present in some groups of the Early Postclassic suggests that this facial phenotype arose in some Mesoamerican group that successfully survived the challenge of the transition to the Postclassic epoch. When observed in the light of further archaeological, paleoenvironmental, epidemiological and linguistic evidences, our analysis supports the notion that the fall of the Classic world was a wide-impact phenomenon, which not only concerned political, religious and social transformations, but also influenced the demography and biology of Central Valley populations.

We thank the authorities of the Dirección de Antropología Física in the Museo Nacional de Antropología, México, for allowing access to skulls in their care. We also thank Dr Christian Klingenberg and two anonymous reviewers who provided thoughtful insights into the methodological aspects of the paper. This research was supported by the Spanish Ministerio de Educación y Ciencia (CGL2004-00903/BTE) and FEDER (Fondos Europeos de Desarrollo Regional), and the Universitat de Barcelona, Beca de Recerca i Docència (to N.M.-A.).