Human bony labyrinth is an indicator of population history and dispersal from Africa.

The dispersal of modern humans from Africa is now well documented with genetic data that track population history, as well as gene flow between populations. Phenetic skeletal data, such as cranial and pelvic morphologies, also exhibit a dispersal-from-Africa signal, which, however, tends to be blurred by the effects of local adaptation and in vivo phenotypic plasticity, and that is often deteriorated by postmortem damage to skeletal remains. These complexities raise the question of which skeletal structures most effectively track neutral population history. The cavity system of the inner ear (the so-called bony labyrinth) is a good candidate structure for such analyses. It is already fully formed by birth, which minimizes postnatal phenotypic plasticity, and it is generally well preserved in archaeological samples. Here we use morphometric data of the bony labyrinth to show that it is a surprisingly good marker of the global dispersal of modern humans from Africa. Labyrinthine morphology tracks genetic distances and geography in accordance with an isolation-by-distance model with dispersal from Africa. Our data further indicate that the neutral-like pattern of variation is compatible with stabilizing selection on labyrinth morphology. Given the increasingly important role of the petrous bone for ancient DNA recovery from archaeological specimens, we encourage researchers to acquire 3D morphological data of the inner ear structures before any invasive sampling. Such data will constitute an important archive of phenotypic variation in present and past populations, and will permit individual-based genotype-phenotype comparisons.

The Late Pleistocene dispersal of modern humans from Africa (also known as “Out-of-Africa II”) has been well documented with fossil evidence (1–3). The dispersal also left a complex genetic signature in recent worldwide populations, reflecting the effects of random drift, selection, and gene flow between dispersing and local populations (4, 5). While molecular data have proven highly effective in reconstructing population history and identifying local adaptations (6–10), phenotypic data have the potential to provide similar but also complementary information. Indeed, modern human skeletal variation has been shown to reflect both population history and local adaptation (11–15). However, patterns of heritable variation tend to be superimposed with nonheritable variation due to developmental plasticity and differences in lifestyle (16, 17). As a result, different parts of the human skeleton (such as the neurocranium, cranial base, face, mandible, teeth, pelvis, and long bones) provide differential insights into population history, local adaptation, and lifestyle (15, 18–25). Results on the influence of each of these factors on different skeletal regions depend to some extent on the morphometric variables and methods of analysis used. Overall, however, it appears that the external morphology of the temporal bone most clearly reflects human population history as it seems to be relatively unaffected by environmental effects, such as local adaptation and phenotypic plasticity (18, 23, 26–29).

Given the current state of research on human phenetic phylogeography, it is worth scrutinizing a wider range of human hard tissue structures for their potential to track population history. Here we assess the potential of the cavity system of the inner ear—the so-called bony labyrinth—as a source of information about human population history. The labyrinth reaches its final morphology by birth (30). Accordingly, postnatal developmental differences and phenotypic plasticity can largely be excluded as potentially confounding effects. Furthermore, the labyrinth is contained in dense cortical bone such that it is typically well-preserved, even in fragmentary archaeological specimens. Archaeological samples are relevant because they provide a temporal perspective on modern human phenotypic variation. Furthermore, the labyrinth is less liable to taphonomic distortion than other skeletal parts because of the compactness of the petrous bone.

We compare patterns of worldwide variation in modern human labyrinth morphology with patterns of variation in genetic markers (SNPs) known to track human population history and assumed to evolve neutrally. Three-dimensional data of bony labyrinths were acquired with computed tomography (CT) from a sample of 221 skeletal specimens representing 22 worldwide populations [19 historical populations (<3,000 y B.P.), and three prehistoric populations (>3,000 y B.P.) from Central Europe (Neolithic: Vaihingen, Germany), Japan (Mesolithic: Jomon), and Indonesia (Paleometallic: Nusa Tenggara Timur = Eastern Lesser Sunda Islands)] (Table 1 and Fig. S1). Labyrinth morphology was quantified with 26 anatomical 3D landmarks (). Genetic data from modern populations were retrieved from the Human Genome Diversity Project (HGDP) database (). Typically, studies comparing genotypic and phenotypic population-level data use standard measures of neutral genetic differentiation (Fst), and analogous measures of phenetic differentiation (Pst) (13, 31). The interpretation of these measures, however, depends on specific demographic model assumptions, and on estimates of parameters, such as heritability, effective population size, and within-population variation, which are all difficult to validate (7, 32, 33). To limit possible model bias, we compare “raw” shape distances derived from geometric-morphometric data with allele-sharing distances derived from SNP data. Multidimensional scaling (MDS) was used to evaluate specimen coordinates in multivariate phenetic and genetic space, respectively. To permit direct comparisons between morphometric data and molecular data, we defined 10 major geographic regions (Table 1) and pooled the data region-wise.

Table 1.

Sample for labyrinth data

Population	src*	Region	n	Longitude	Latitude	Dd†
South Africa	Z	SAf	19	21	−12	2.8
North Africa	Z	NAf	7	28	20	7.9
Egypt	Z	NAf	10	31	30	8.9
Western Asia	Z	WAs	8	55	41	12.1
Central Europe	Z	Eur	15	9	47	13.3
Vaihingen	L	Eur	16	9	49	13.3
Northern Europe	Z	Eur	6	30	59	14.1
India	Z	WAs	7	79	21	15.9
China	Z	NAs	5	114	38	19.5
Southeast Asia	Z	SAs	8	98	13	19.6
Indonesia	Z	SAs	8	107	−2	22.4
Ainu	T	Jap	10	141	43	23.5
Japan	T	Jap	10	140	36	24.6
Jomon	T	Jap	8	138	35	24.6
Sunda	A	SAs	4	120	−9	25.0
Australia	D	Aus	11	125	−21	27.2
Torres Strait	D	Aus	8	143	−9	28.2
Arctic	Z	NAm	19	−118	64	31.6
North America	Z	NAm	10	−117	39	33.3
Mesoamerica	Z	SAm	7	−91	15	37.7
South America	Z	SAm	15	−69	−9	43.7
Patagonia	Z	SAm	10	−70	−53	49.5

Collection source: A, Airlangga University, Surabaya; D, Duckworth Collection, Cambridge; L, Landesdenkmalamt Baden-Württemberg, Stuttgart; T, Tokyo University Museum, Tokyo; Z, Department of Anthropology, University of Zurich.

Dispersal distance in 1,000 km.

Results

Fig. 1 shows morphological variation of the human bony labyrinth, revealing a fairly isotropic spatial pattern of variation of the sample around the mean shape. Fig. 1 shows patterns of variation in multivariate shape space. The first two shape components (SCs) account for 10.1% and 8.8% of the total sample variance, respectively (see Fig. S2 for variance proportions of all SCs). Along SC1, populations exhibit an approximate geographical ordering from Africa (blue tones in Fig. 1) to the Americas (red tones) (correlation between SC1 and dispersal distance from Africa: r2 = 0.58). Superposition of phenetic and genetic data (Fig. 1) shows closely matching distributions in multivariate space, and also reveals the geographic ordering of genetic data known from earlier studies (8, 34, 35). Within-group variation of labyrinth morphology is wide compared with between-group variation (Fig. 1 and Table 2), and this is also the case for the genetic data (Table 2).

Fig. 1.

Worldwide variation of human bony labyrinth shape in physical space (A) and in multivariate shape space (B and C). (A) Visualization of the mean shape of the bony labyrinth, and patterns of variation around each landmark (data points and 3 SD ellipsoids). (Scale bar, 5 mm.) (B) For each population the location of the mean shape is indicated with a diamond, and the range of variation with a 1 SD ellipsoid; the spectral colors indicate dispersal distance from Africa (see Table 1 for details). (C) Procrustes superposition of genetic data (blue) on phenetic data (red: region-mean locations; gray: population-mean locations, as in A); note that the distribution of populations along SC1 and SC2 approximately reflects their geographic distribution. Arct, Arctic; Aust, Australia; CEur, Central Europe; Chin, China; Egyp, Egypt; Indi, India; Indo, Indonesia; Japa, Japan; Jomo, Jomon; MAme, Mesoamerica; NAfr, North Africa; NAme, North America; NEur, Northern Euope; Pata, Patagonia; SAfr, South Africa; SAme, South America; SEAs, Southeast Asia; Sund, Sunda; Torr, Torres Strait; WAsi, Western Asia; Vaih, Vaihingen.

Table 2.

Within-population and between-population variances

Source	Within population (Vw)*	Between population (Vb)	Vb/Vw
Labyrinth SC1–SC72	0.0053	0.0014	0.2680
HGDP 8192 SNP data	0.0328	0.0156	0.4759

Averaged over all populations.

Patterns of variation in the full set of shape components can be used to assess whether labyrinth morphology follows a neutral-like pattern of evolution and population differentiation. The neutral hypothesis posits that between-group variation evolves as a linear function of within-group variation (36, 37). Indeed, the regression of between-population on within-population variances of shape components does not deviate significantly from linearity (Fig. S3).

The phenetic distance between sub-Saharan African populations and populations outside Africa (Dp) increases with dispersal distance from Africa (Dd), and there is close correspondence between phenetic (Dp) and genetic (Dg) distance functions (Fig. 2). Together, these results indicate that modern human labyrinth morphology bears an isolation-by-distance signature that corresponds to the dispersal-from-Africa signature known from genetic data. The phenotypic signal can be recovered from comparatively small population samples (∼10 individuals per population, 22 populations), and it can be better quantified with individual-based morphometric shape distances than with population-based Pst distances (Fig. 2 and Fig. S4). Additionally, the genotypic signal can be quantified equally well with individual-based allele-sharing distances as with population-based Fst distances (Fig. 2 and Fig. S4).

Fig. 2.

Correlations between phenetic, genetic, and dispersal distance data. (A) Phenetic distance Dp (red) and genetic distance Dg (blue) versus dispersal distance Dd; note that slopes for phenetic and genetic data are statistically similar (P > 0.74). (B) Dp vs. Dg by region (to facilitate comparisons, Dp and Dg are log-mean-centered). See population name abbreviations list in Fig. 1.

The labyrinth morphologies of the prehistoric populations show both commonalities and differences compared with those of modern populations at geographically similar locations (Figs. 1 and 3). The early populations from Europe and Japan group most closely with corresponding modern populations, indicating a likely population continuity since the Neolithic. On the other hand, the early Indonesian sample groups with Austromelanesian populations rather than with today’s predominant (Malay) populations in the same region, indicating some level of population replacement in the recent past, as suggested in an earlier study (38).

Fig. 3.

Consensus tree of phenetic similarities among populations. Numbers indicate branch support as evaluated from 1,000 resamplings on random subsamples of n = 6 specimens per population (note that local population affiliations are better supported than global affiliations; for details, see ). See population name abbreviations list in Fig. 1.

Discussion

The morphology of the bony labyrinth has been shown to be a complex indicator of both functional (i.e., locomotor and auditory) adaptations (39–42; but see ref. 43) and phyletic relationships among primates and other mammals (40, 44). Recent studies show potential to elucidate morphological differentiation at the level of hominoid species (42) and subspecies (45). Morphometric quantification, such as that used in the present study, differentiates labyrinth morphology among great ape species and gibbon (42), and between subspecies of chimpanzees (45), although the contributing factors to such shape differences remain to be elucidated. Patterns of intraspecific variation of the primate labyrinth have only recently been investigated (46), and the observation that slow-moving species exhibit higher levels of intraspecific variation than fast-moving species has been interpreted as evidence for relaxed stabilizing selection in slow-moving species. Compared with great apes and fossil hominids, such as Australopithecus and Paranthropus, the bony labyrinth of the genus Homo is characterized by relatively large superior and posterior semicircular canals, and a relatively small lateral canal. Whereas this pattern has been related to the evolution of obligate terrestrial bipedalism or bipedal running in Homo (39, 47), other factors, such as neurocranial developmental integration not related to vestibular function (30, 42) and neutral genetic drift (43, 44), may have led to differences among closely related taxa. Our results show that, in addition to functional and phyletic information, the labyrinth morphology of modern humans contains significant phylogeographic information; it exhibits a geographically structured pattern of worldwide variation similar to that observed in neutral genetic data and attributed to the dispersal from Africa.

The mixed functional, phyletic, and phylogeographic signal evinced by labyrinth morphology raises the question as to how these factors contribute to the overall pattern of shape variation of the bony labyrinth of modern humans. Labyrinth morphology is constrained by the physics of vestibulo-cochlear signal transduction (48), so it is expected to be under taxon-specific–stabilizing selection that optimizes auditory, motion, and gravity perception. How can the action of stabilizing selection be reconciled with the evidence presented here that labyrinth morphology shows a neutral-like pattern of worldwide variation? This question expresses a well-known paradox of phenotypic evolution that “many quantitative traits show substantial heritable variation and yet appear to be subject to stabilizing selection” (49). One basic assumption is that the constraints exerted by stabilizing selection are similar for all human populations worldwide, and are largely independent of the local environment. In addition, as an effect of in vivo neural plasticity, an individual’s sensory performance tends to be optimized within the constraints imposed by the signal transduction system (50). Under these conditions, the classic model of Lande (51, 52) states that the expected variance of a phenotypic feature under stabilizing selection depends on the width of the adaptive zone w around a single fitness peak (with w2 >> σ2, where σ2 is the within-population phenotypic variance). Under these model assumptions, the patterning of population-mean variation is due to random genetic drift: that is, “stochastic wandering of populations in the vicinity of the fitness peak” (53).

An additional model prediction is that within-population variation will be larger than variation between population means (52) (Table 2). We thus expect that—within the human-specific constraints given by w—the morphology of the bony labyrinth exhibits nondirectional variation according to a neutral-like evolutionary process. The results presented here largely correspond to this pattern, indicating no dominant direction of variation and wide intragroup compared with intergroup variation. Furthermore, our data permit a tentative estimate of the relative width of the adaptive zone, w/σ (). Assuming effective population sizes between 5,000 and 10,000 individuals (54), estimates for w/σ range between 39 and 55 (Table S2). These values are larger than the theoretical value of 10 assumed in Lande (51), but within the range of empirical values reported in the literature (55).

Various models have been proposed to study genetic and phenetic evolution in a population expanding over space and time (52, 56, 57). The neutral expectation is that genetic and phenetic distances between the founder population and subsequent populations increase as a function of time and dispersal distance, and that deviations from this relationship indicate nonneutral effects. Fig. 2 shows a similar increase of genetic and phenetic distance with dispersal distance, providing support for the notion that labyrinth morphology evolves neutral-like. The phenetic data exhibit relatively more variation around the regression line than the genetic data. It remains to be clarified to what extent these deviations reflect factors, such as phenotypic plasticity, nonadditive genetic variation, and local adaptation, and to what extent they merely reflect sampling bias due to the limited number of independent phenetic features or small population-specific samples. An earlier study on worldwide modern human cranial variation showed that the interplay of these factors tends to complicate the discrimination between signatures of “neutral” population history and local adaptation, and that a combination of genomic, phenomic, environmental, and fossil/archaeological data will be required to reconstruct the diversity of evolutionary processes that have shaped modern human cranial variation (15).

Additional evidence for a strong influence of genetic drift, and thus a neutral-like pattern of human labyrinth variation, comes from the observation that the first shape component (SC1 in Fig. 1) exhibits correlation with dispersal distance from sub-Saharan Africa (r2 = 0.55). This is a typical property, revealed by dimension reduction methods, of multivariate datasets exhibiting seriation in time and space (58, 59). Analogous effects have been described for neutral genetic data in which similarity decreases with dispersal distance, be it as an effect of range expansions beyond Africa, or isolation-by-distance within continents (35, 60). The correlation between SC1 and geographic distance primarily reflects isolation by distance, but it also permits inferences on the geographic origin of the dispersal ( and Figs. S5 and S6).

In view of these findings, it is worth reexamining skeletal structures, such as the cranium, for evidence of neutral-like evolution under stabilizing selection. For example, the evolutionary divergence between human and Neanderthal cranial morphologies has been characterized as compatible with neutral expectations, whereas the divergence between cranial morphologies of chimpanzee subspecies is smaller than expected under neutrality, thus indicating stabilizing selection and developmental constraints reducing within-group variation (61, 62). These findings may indicate a wider adaptive zone w/σ in Homo than in Pan, probably an effect of cultural buffering favoring relaxed stabilizing selection (62, 63).

Given the strong population-history signal of labyrinth morphology, one may ask whether the labyrinth of a single individual can be used to infer its population affiliation. This is clearly not the case, as has been shown previously for cranial morphologies (64). Similar limitations apply to inferring population affiliations of isolated fossil hominin labyrinths (65–67). One reason is that within-population variation is larger than between-population variation. The other is that the highly integrated nature of the phenotype results in a small number of independently evolving features (the first few SCs with significant loadings) available to discriminate between populations. As has been shown for genetic data, the analysis of population structure faces a trade-off between the number of available markers and sample size (68). As demonstrated here, evaluating local population-mean phenotypes is an effective means to infer population affiliations and history, even for comparatively small samples.

We have therefore shown that the petrous portion of the human temporal bone is a valuable source of information for “paleo-phenomics.” In recent years it was shown that the petrous bone is a valuable source of information for paleo-genomics as well (69). Indeed, the dense cortical bone surrounding the labyrinth has become a major target of ancient DNA and collagen retrieval from archaeological and fossilized human specimens (69, 70). Typically, destructive sampling is performed before 3D data acquisition with CT, with the effect that unique phenotypic information is irrevocably lost. We thus urge all parties involved in extracting organic material from structures, such as the petrous bones or teeth, to acquire high-resolution CT data beforehand. The X-ray doses used in medical CT as well as in low-dose μCT (<200 Gy) have been reported as not affecting ancient DNA preservation in the sample (see ref. 71, pending confirmation, including blind tests by independent laboratories). Combining phenotypic and genotypic information obtained from the same specimen and the same skeletal structure will thus open new venues for the coinvestigation of genetic and phenetic aspects of human evolutionary processes.

Materials and Methods

Morphometric Analysis of the Bony Labyrinth.

Volumetric data of the bony labyrinths were acquired with medical and μCT from a global sample of modern humans (n = 221; 22 populations; institutional review board approvals were unnecessary). When available, left labyrinths were used, otherwise mirror-imaged right labyrinths. Labyrinth shape was quantified with a set of 26 3D landmarks equally distributed over the semicircular canals, the vestibular region, and the cochlea, as specified in ref. 44 and visualized in Fig. 1. Shape variation in the sample was analyzed with principal components analysis (PCA) of shape (72). Note that this procedure is equivalent to the evaluation of principal coordinates via MDS analysis of the between-individual Procrustes distance matrix. Population-mean shapes were used to evaluate a between-population phenetic (Procrustes) distance matrix. The neighbor-joining method was then used to derive a phenetic similarity tree. To obtain statistical estimates of branch support, 1,000 trees were constructed with n = 6 randomly selected specimens per population (resampling with replacement), and a consensus tree was evaluated [procedures “neighbor” and “consense” of the PHYLIP software (73)] (Fig. 3).

Modern Human Genetic Data.

Genetic data were obtained from the HGDP panel (74). The following filtering procedures were applied: First, we removed individuals that have been labeled as related (75), then we used the pipeline reported in Gazal et al. (76) to exclude individuals having up to third-degree relationships. Next, we ran ADMIXTURE (77) to identify and filter out individuals exhibiting recent population admixture. Using the PLINK software (78), we selected only the autosomal biallelic SNPs with minor allele frequency >0.05 or <0.01% missing data. We applied a Hardy–Weinberg equilibrium threshold by population (HWE P < 10−5) to exclude SNPs deviating from it, and removed SNPs that were fewer than 2,000 bases apart from each other. Following these quality-control procedures, we obtained genotype data on 351K SNPs for 902 individuals from 52 populations (Table S1). Finally, we used the Genomic Evolutionary Rate Profiling method (79) to obtain site-specific estimates of evolutionary constraint. After this step, we retained a subsample of 8,192 (=213) SNPs that are unlikely to be affected by natural selection. The average distance between these SNPs is 0.345 Mb.

These data were analyzed with methods closely matching the methods used for the analysis of the phenetic data. First, we calculated allele-sharing distances (which, in terms of information theory, are Hamming distances) between all individuals, resulting in an n × n interindividual distance matrix (no Jukes–Cantor correction was performed, which corresponds to the assumption that back-mutations or mutation saturation are not present in the sample). MDS was used to transform this matrix into specimen coordinates in multivariate genetic space (principal coordinates). These data were used to evaluate population-mean locations, and between-population Euclidean distances, using the same methods as for the phenetic data. Note that the procedures used here differ in some aspects from traditional genetic PCA as proposed by Patterson et al. (80), where each individual’s position in multivariate genetic space is directly coded relative to a predefined reference (the most frequent alleles). Here, each individual’s position in that space results from MDS on the interindividual allele-sharing distance matrix.

Dispersal Distances from Africa.

Traditionally, geographic distances from Africa are estimated as sums of great-circle segments connected through geographic way-points (8). These estimates implicitly assume directional migration paths rather than dispersal. The estimates used here are based on a nondirectional population growth and diffusion model (the “Fisher–Kolmogorov–Petrovsky–Piskounov” model), which was implemented in an agent-based simulation framework (81, 82). We simulated the global spread of a human population originating in South Africa (Johannesburg) on a spherical grid with 660,492 nodes (internode distance 30–35 km) and a total number of 70,000 1-y time steps. This simulation strategy permits evaluating, for each habitable grid node, the accumulated dispersal distance Dd from Africa (Fig. S1). In a next step, geo-locations for all individuals in the phenetic and genetic samples were evaluated. For the labyrinth sample, geo-locations were estimated from the population affiliation information of each individual, using the present-day geographic center of each population as the reference location. For the HGDP sample, such information was gathered from the metadata associated with each individual. Finally, each individual’s geo-location was used to read out from the simulation map its corresponding dispersal distance.

Comparison of Distance Data.

To permit direct comparisons between genetic and phenetic distances (Fig. 2), populations were further grouped according to larger geographical regions (Fig. 2, Table 1, and Table S1). The group-wise phenotypic and genotypic distance metrics used here differ from analyses relying upon Pst and Fst distances (83, 84). These latter distance measures depend on estimates of parameters, such as effective population size Ne, heritability h2, and within-population variation, while the measures used here represent distances between groups of individuals from geographically confined areas. For phenetic data, an additional advantage of the measurements used here is that possible differences between female and male labyrinth morphologies (85) tend to be averaged out. Nevertheless, to permit comparisons with Pst–Fst-based studies, we analyzed our data with the respective methods (83, 84), and provide the results in Fig. S4. To facilitate scale-independent comparisons between phenetic, genetic, and dispersal distances, all data were log10-transformed (51).

Estimation of the Strength of Stabilizing Selection on Labyrinth Shape.

Under weak stabilizing selection, the expected equilibrium variance σ2mean of the mean phenotype around its optimum is σ2mean = (w2 + σ2w)/2N, where w2 characterizes a Gaussian fitness function for stabilizing selection, σ2w is the within-population variation, and N is the effective population size (51). Resolving by w2 and standardizing by σ2w yields the relative width of the adaptive zone, w/σ = (2Neσ2mean/σ2w − 1)1/2. To obtain estimates of w/σ we use the following approach: σ2w and σ2mean are estimated by the average of the within-population variances Vw, and the between-population variance Vb, respectively. PCA of shape (see above) yields Vw and Vb for each shape component SC, as well as for labyrinth size S. Because SCs are statistically independent of each other (and, in the present case, of size S), each of these variables can be treated as an independent phenotypic feature. Accordingly, w/σ can be estimated for each feature (Table S2). Because variation in SC1 is largely an effect of dispersal distance from Africa, and thus unlikely to be in an equilibrium state (57), we omit this component from our analyses. Published estimates of effective human population size Ne vary widely; here we use a range of Ne = 5,000–10,000 (54). Together, these data yield the estimates of w/σ reported in Table S2.