Homo erectus is the first hominin species with a truly cosmopolitan distribution and resembles recent humans in its broad spatial distribution. The microevolutionary events associated with dispersal and local adaptation may have produced similar population structure in both species. Understanding the evolutionary population dynamics of H. erectus has larger implications for the emergence of later Homo lineages in the Middle Pleistocene. Quantitative genetics models provide a means of interrogating aspects of long-standing H. erectus population history narratives. For the current study, cranial fossils were sorted into six major palaeodemes from sites across Africa and Asia spanning 1.8–0.1 Ma. Three-dimensional shape data from the occipital and frontal bones were used to compare intraspecific variation and test evolutionary hypotheses. Results indicate that H. erectus had higher individual and group variation than Homo sapiens, probably reflecting different levels of genetic diversity and population history in these spatially disperse species. This study also revealed distinct evolutionary histories for frontal and occipital bone shape in H. erectus, with a larger role for natural selection in the former. One scenario consistent with these findings is climate-driven facial adaptation in H. erectus, which is reflected in the frontal bone through integration with the orbits.

1. Introduction

Homo erectus occupies a central position in human evolution. Advances in dating species divergences from ancient DNA position H. erectus (or some portion of the lineage) as a potential last common ancestor of the modern human (Homo sapiens) and Neanderthal–Denisovan lineages in the Middle Pleistocene. If the latter two groups split prior to 500–700 ka [1–3], then this postdates the vast majority of fossils assigned to the other candidate for this position, Homo heidelbergensis s.l. Recent H. sapiens and H. erectus are also the most geographically expansive species of hominins. The current study provides novel insights into relative genetic diversity and the microevolutionary processes that accompanied population differentiation of H. erectus by applying methods grounded in population history theory to 27 H. erectus cranial fossils and a comparative sample of approximately 300 recent H. sapiens individuals.

Homo erectus shares with H. sapiens some key features of evolutionary history, including migration out of Africa to occupy a range of more seasonal and temperate habitats across Eurasia and Asia. Homo sapiens therefore serves as a useful model for thinking about H. erectus population history. The genetic landscape of H. sapiens was shaped by a major out-of-Africa population bottleneck [4], serial founding effects [5], limited gene flow across their range [6,7] and local adaptation [8,9]. These same population history events left an imprint on the species' cranial phenotype [5,10–12].

Homo erectus, broadly defined, spans 1.8 Myr and traverses latitudes from 25 °S to 40 °N in Africa, Eurasia and Asia [13–15]. Conventional wisdom holds that H. erectus originated around 1.9 Ma in East Africa [16], where there are some of the oldest H. erectus sites and abundant evidence of more ancient Homo species (figure 1). From there, populations of African H. erectus dispersed east into Asia, a migration that included the ephemeral occupation of W. Asia (Caucasus) and southern China by 1.8 Ma [17–21] (but see [22–24]). East and Southeast Asian groups from China and Java, respectively, are typically viewed as the product of a single migration event (but see [25]). Many scenarios imply periods of regional isolation that facilitated in situ phyletic transformation through genetic drift and/or environmental adaptation but with sufficient gene flow to maintain species cohesion [15,23,25–34]. Some workers also discern a north-south cline of phenotypic variability in the Asian part of the range [35]. Most workers view these spatial and temporal variations as population- or subspecies- rather than species-level phenomena [13,23,26,36,37] (but see [38–41]).

Figure 1. Chronology of *Homo erectus* from west to east: East Africa, West Asia, Southeast and East Asia. Black rectangles indicate time ranges of samples included in each palaeodeme, while thinner grey rectangles indicate the full time range of fossils assigned to *H. erectus* in those regions. Key *H. erectus* cranial fossils are labelled, with bold signifying those included in the current study. Abbreviations as in table 1 or as follows: OH: Olduvai Hominid; BSN: Busidima North; DAN: Dana Aoule North. The question mark indicates some uncertainty regarding the attribution of Daka (BOU-VP-2/66) to *H. erectus*. Inset: landmarks and semilandmarks from the frontal bone (yellow) and occipital bone (green) superimposed on the Sambungmacan 3 crania. Landmarks are indicated by larger spheres than semilandmarks. (Online version in colour.)

Table 1. *Homo erectus* fossils included in analysis, organized by palaeodemes.
Collapse
population	sites	sample^a	age of sample
West Asia (WAS)	Dmanisi	Both: D2280, D2700, D3444	1.83–1.77 Ma [17,68]
Early East Africa (EAF)	E. Turkana	Both: ER 3733, ER 3883O: ER 730	1.63–<1.5 Ma [69–71]
East Asia (EAS)	Zhoukoudian	Both: Zkd 3, Zkd 10, Zkd 11, Zkd 12	0.75–0.8 Ma [33]
Early Southeast Asia (ESA)	Sangiran	Both: S17, Skull IX (Tjg-1999.05)O: S3, S10, S12F: S2	0.9–0.79 Ma [72,73]
Late Southeast Asia-S (LSA-S)	Sambungmacan and Ngawi	Both: Sm 1, Sm 3, Sm 4, Ngawi	Middle Pleistocene [74]
Late Southeast Asia-N (LSA-N)	Ngandong	Both: Ng 6, Ng 11, Ng 12O: Ng 1, Ng 7F: Ng 5, Ng 10	0.12–0.11 Ma [75,76]

^aBold indicates that a fossil was included in the occipital and frontal bone analyses. O and F indicate inclusion only in the occipital bone or frontal bone analysis, respectively. Abbreviations as follows: D: Dmanisi; ER: East Rudolf; Zkd: Zhoukoudian; S: Sangiran; Sm: Sambungmacan; Ng: Ngandong.

However, there are long-standing debates surrounding the alpha taxonomy of H. erectus (reviewed in [26,42]). At the extremes, there are small minorities of workers that advocate subsuming several early Homo species, including H. erectus, into a single species [19] or splitting the traditional H. erectus sample into many distinct species [43]. Most workers, however, view H. erectus as either a broadly distributed and polytypic species [13,14,23,44] or a specialized Asian lineage, in which case the earlier African and possibly Georgian fossils are assigned to a separate species, H. ergaster [16,45,46]. The degree of morphological variation in H. erectus is at the high end of the range documented in extant species, but is relatively low considering the geographical and temporal range of this species [26] and, moreover, its neurocranial shape is distinct from both earlier and later archaic Homo species [15]. Additionally, dentognathic anatomy in particular has hinted at a more complex set of relationships among regional groups through time than is captured by a simple east-west dichotomy [24,47,48]. For these reasons, H. erectus is defined here as including fossils from Africa, Eurasia and East/Southeast Asia, with the recognition that there are other plausible taxonomic hypotheses.

The extent to which recent H. sapiens and H. erectus converged in their population history has not been explored previously. Statistical models grounded in population history theory provide a powerful set of tools to detect past signals of microevolutionary processes that figure prominently in these population history scenarios. Quantitative genetics is concerned with the evolution of complex traits, including phenotypic traits determined by polygenic inheritance, such as cranial shape. Evolutionary morphological analyses adapted from quantitative genetics are already yielding valuable insights into human evolution, including recognition of a greater role for neutral processes (e.g. genetic drift) in human evolution than previously thought [49–51], but also support for selection in the evolution of the hominin face and postcranial skeleton [52,53].

The current study assesses the hypothesis that similar, primarily neutral, microevolutionary factors related to shared evolutionary history shaped cranial morphology in H. erectus and recent H. sapiens. This hypothesis is an oversimplification; the details of population size, gene flow, genetic drift, mutation and selection certainly varied between the species. However, the predictions are generalized rather than granular and will reflect the most common pattern over time rather than nuanced details of population history.

—	Prediction 1: H. erectus exhibits a similar magnitude of intraspecific cranial variation as recent H. sapiens. This prediction follows from the observation that intraspecific cranial variation is correlated with neutral genetic variation across hominoids, including H. sapiens [54], indicating an important role for population history (e.g. genetic drift, gene flow) in shaping cranial variation.
—	Prediction 2: H. erectus populations are more divergent than recent H. sapiens populations. This is the expectation under mutation-drift equilibrium if population divergence in both species is due primarily to neutral processes (e.g. genetic drift) given the greater time since divergence among H. erectus populations [55].
—	Prediction 3: neutral evolutionary processes (e.g. genetic drift) can account for phenotypic divergence among groups of H. erectus and H. sapiens. This prediction is consistent with long-standing ideas about geographical and genetic isolation of regional H. erectus populations [32] and empirical evidence for neutral cranial diversification in other Homo species, including H. sapiens [56].

The hypotheses were assessed separately for the frontal and occipital bones both to maximal fossil sample sizes and to assess mosaic cranial evolution. Differential selection on the face and brain could yield divergent evolutionary histories in the frontal and occipital bones, respectively, given modularity between braincase and face and/or generally low levels of cranial integration observed for H. sapiens [57–59]. These analyses represent the first use of quantitative genetics to shed light on the evolutionary processes that shaped H. erectus cranial variation.

2. Methods

Three-dimensional landmarks and semilandmarks acquired from the occipital and frontal bones (see figure 1 and electronic supplementary material, Supplemental Methods) in H. erectus and recent H. sapiens are the raw data used in subsequent analyses. Standard missing landmark estimation procedures were used for each bone, including reflected relabelling and TPS based on species-specific means [60] (electronic supplementary material, Supplemental Methods and tables S1 and S2). Semilandmarks were slid to maximize geometric homology by minimizing the bending energy matrix [61] and all landmarks/semilandmarks were superimposed by generalized Procrustes analysis separately for each bone [62]. All analyses were performed in the R statistical environment, relying heavily on the ‘Morpho’ [63], ‘geomorph’ [64], ‘sampling’ [65] and ‘ggplot2’ [66] packages.

The H. erectus sample consisted of n = 23 for the occipital analysis and n = 22 for the frontal analysis (n = 27 total) subdivided into six spatially and temporally circumscribed palaeodemes (sensu [67]) of 2–6 individuals each (table 1 and figure 1). Some well-preserved fossils were excluded to constrain the temporal depth of each sample to less than or equal to 130 kyr (see [77] for a similar approach). The late Indonesian samples from Ngandong and Sambungmacan/Ngawi were treated as distinct demes due to subtle differences in discrete traits and neurocranial shape between these groups [30,78]. Small fossil samples introduce uncertainty into estimates of palaeodeme averages, but relatively strong morphological homogeneity within palaeodemes [15] suggests that we are sampling near the group centroids in most cases. Given the controversy surrounding H. erectus species composition, statistical analysis was re-run for H. erectus s.s. from Asia only and presented in the electronic supplementary material.

(a) Prediction 1

Recent H. sapiens from six populations (average of n = 32 each, total n ∼ 193) representing a similar geographical distribution as H. erectus were used to evaluate Prediction 1 (electronic supplementary material, table S3). Several populations were insular groups that were isolated for periods of their history, perhaps mirroring the geographical isolation of H. erectus demes.

Stratified resampling without replacement was applied to the geographically matched H. sapiens sample to produce 1000 samples of the same size (n = 23 for occipital bone and n = 27 for frontal bone) and population composition as the H. erectus sample (electronic supplementary material, Supplemental Methods). The sum of eigenvalues (SEV) was calculated from both the original H. erectus sample and the 1000 H. sapiens samples to test Prediction 1. The SEV is equivalent to the sum of all trait variances and is used here to measure intraspecific shape variation. This measure scales linearly with neutral genetic variation across hominoids [54]. The H. erectus value is compared to the empirical distribution of modern H. sapiens values to determine the probability that H. erectus exhibits similar intraspecific variation as H. sapiens.

(b) Prediction 2

Centroids for H. erectus palaeodemes and H. sapiens populations were calculated using the same samples described above (electronic supplementary material, table S3), again using resampling to generate H. sapiens samples of equal size as the H. erectus demes. Between-deme variation in H. erectus was calculated using SEV calculated from the covariance matrix of palaeodeme averages. This H. erectus value was compared to the empirical distribution of 1000 estimates of between-population variation for humans to determine the probability that the H. erectus value exceeded that of H. sapiens.

(c) Prediction 3

The quantitative genetics tests for Prediction 3 require an estimate of the G (or P) matrix, which contains the within-population additive genetic (or phenotypic) variances and covariances of traits and is essential for formulating predictions regarding the magnitude and direction of between-population variation under certain microevolutionary conditions [79]. The small sample sizes of H. erectus preclude direct estimates of either the G or P matrix. Therefore, the recent H. sapiens pooled within-population P matrices for the occipital and frontal bones were used as proxies for H. erectus. H. sapiens was chosen because it is the closest extant phylogenetic match for H. erectus [49,52,80] and the human P matrix is a good estimate of its G matrix [81–83]. In total, 145 recent H. sapiens from three populations with n ∼ 50 each were used to calculate the pooled within-population covariance matrices that were used as proxies for H. erectus in testing Prediction 3 (electronic supplementary material, table S3). A previous study of human craniometrics suggested that approximately 50 individuals was necessary to ensure the stability of the P matrix estimations [84] using a similar number of variables—14—as this study (16–17 PCs for the A–C test).

Two tests were used to evaluate Prediction 3. First, p_max (also known as the line of least evolutionary resistance [85]) was calculated as PC 1 of the pooled-within population covariance matrix for recent H. sapiens. Populations are predicted to diverge primarily along p_max if their differentiation occurred via neutral evolutionary processes [85–87]. This was assessed by calculating the angle as the arccosine of the dot product of the normalized vectors for p_max and PC 1 of a between-group PCA of the H. erectus palaeodeme means or H. sapiens population means. This angle measures divergence between the primary axes of within- and between-group variation. This is a conservative test as p_max (PC 1) is likely to be conserved across species to a greater extent than more minor PC axes. However, evolutionary change can follow an axis of high intrapopulation variability that is not perfectly aligned with p_max if the variance is spread across several dimensions [88].

Therefore, the Ackerman–Cheverud (A–C) test was used to evaluate whether divergence among populations (B) is proportional to within-population variation (P) across multiple dimensions [49]. First, the centroids of all six H. erectus palaeodemes (and all six H. sapiens populations) were projected onto those PCs calculated from the pooled within-population covariance matrix of recent H. sapiens accounting for greater than or equal to 1% of the variation. Next, the natural logarithm of between-deme (or -population) variance on each PC was regressed of the natural logarithm of within-population variance [49]. Under a pure drift model, the regression coefficient (β) is expected to be 1, and deviations from this expectation can indicate non-random processes. Both of these tests evaluate an average signal which may mask some episodes of adaptive evolution. The A–C test is a relatively conservative test and its power is reduced when dealing with a small number of groups, so any significant results are likely to be meaningful [49]. Both p_max and the A–C test assume that the P matrices of H. erectus and H. sapiens are proportional. This conjecture finds some support in the generally conserved nature of the G matrix within Homo [89,90], although some minor differences in integration between archaic and modern Homo have been documented [89,91].

3. Results

(a) Shape variation

The primary axes of variation for frontal and occipital shape separated the two species (H. erectus and H. sapiens) (figure 2a,b). Shape differences reflect well-documented distinctions between the species: H. erectus has a flatter frontal squama with a taller supraorbital torus and greater constriction posterior to the supraorbital torus, while the occipital bone is relatively wider but more tightly angled in its midline profile (figure 2c,e). The pattern of intraspecific variation in H. erectus differed for the two bones. The greatest contrast in occipital bone shape (PC 1) was between the geologically oldest (WAS) and youngest palaeodemes (LSA-N and LSA-S) (electronic supplementary material, figure S1a). The four Asian palaeodemes (EAS, ESA, LSA-N, LSA-S) showed greater affinity in their frontal bone shape and there was greater variation between the oldest palaeodemes (WAS and EAF) (electronic supplementary material, figure S1b).

(b) Prediction 1 intraspecific variation

The H. erectus frontal bone SEV was 67% higher than the average SEV for recent H. sapiens and exceeded 100% of the resampled H. sapiens values (figure 3c,d). The H. erectus occipital SEV value was 27% larger than the average H. sapiens SEV and was higher than 98% of those values. Restricting the analysis to only H. erectus s.s. from Asia yielded a similar pattern, but the within-deme variation in occipital shape was not significantly higher in H. erectus (electronic supplementary material, Supplemental Results).

(c) Prediction 2: between-group variation

The between-group variation in frontal bone shape was approximately twice as great in H. erectus as in recent H. sapiens and the H. erectus value exceeded 100% of the resampled H. sapiens values (figure 2c,d). The between-group variation was 50% higher for the occipital bone in H. erectus and the H. erectus value exceeded 99% of the recent H. sapiens values. Restricting the analysis to only H. erectus s.s. from Asia produced similar results (electronic supplementary material, Supplemental Results).

(d) Prediction 3: neutral evolution of cranial shape during population divergence

Homo erectus populations diverged along the p_max calculated from the human sample, but also exhibited significant spread along other axes (figure 3). The shape vectors representing p_max and the direction of maximal between-deme variation in H. erectus are moderately aligned for occipital shape (angle = 51.9°; Pearson's r = 0.60), but more orthogonal for frontal bone shape (angle = 75.0°; Pearson's r = 0.22). The angle between p_max and the axis of maximal between-population variation in recent H. sapiens was 61° (r = 0.47) for the occipital bone but only 46° (r = 0.70) for the frontal bone.

Figure 3. Ordination of the first two principal components of within-population variation in (a) occipital and (b) frontal bone shape. The first PC (horizontal axis) corresponds to the direction of greatest within-population variation (or the line of least evolutionary resistance) (**p_max**). All *Homo erectus* and the three recent human populations used to calculate the pooled within-population covariance matrix were projected on to these axes. Colours for *H. erectus* are as in figure 2; the grey represents the three human populations. Abbreviations for labelled fossils are from table 1. Intraspecific variation in *H. erectus* (red dashed line) compared to resampled values for recent *H. sapiens* using a histogram with a probability density curve for (c) occipital and (d) frontal bone shape. (Online version in colour.)

Less than a quarter of the total variation in occipital and frontal shape is concentrated onto their respective p_max vectors. The A–C test was applied to 17 (occipital) or 16 (frontal) PC axes, accounting for greater than 90% of total shape variation in each bone. The slopes for the occipital and frontal bones are β = 0.97 (R² = 0.65) and β = 0.59 (R² = 0.50), respectively, for H. erectus (electronic supplementary material, figure S3). The slope was not significantly different from 1.0 for the occipital bone (p = 0.89), but was for the frontal bone (p = 0.02). A slope less than 1.0 implies stabilizing selection in directions of high within-population variance and/or directional selection in directions of limited within-population variance. A slope of 1.0 was not rejected for the occipital (β = 1.11, p = 0.54) or frontal bone (β = 1.19, p = 0.53) in recent H. sapiens. These same tests of neutrality yielded similar results when applied to Asian H. erectus only, although the H. erectus frontal bone β did not differ significantly from 1.0 (electronic supplementary material, Supplemental Results).

4. Discussion

The current study examined population history signals in the frontal and occipital bone for two species, H. erectus and recent H. sapiens. These species are unique among hominins in experiencing a major migration out of Africa to occupy more diverse habitats across Eurasia and East and Southeast Asia. The species differ in that H. erectus has a much deeper fossil record of nearly 2 Myr and was more constrained in its latitudinal distribution across this range. Results indicate limited support for shared population history in H. erectus and recent H. sapiens and mosaic evolution of the cranium in H. erectus.

The results highlight distinct evolutionary histories for the frontal and occipital bones in H. erectus. This is supported by the stronger signal of natural selection in frontal bone diversification among populations compared to the occipital bone, but also the discordant patterns of shape differentiation among demes for these two bones (electronic supplementary material, figure S1). Homo erectus and modern humans themselves diverged more or less along the axis of maximal within-population variation—or p_max—for frontal bone shape, but nearly orthogonal to p_max for occipital bone shape, suggesting different species histories for these regions as well. Modularity between the cranial vault and the face may have facilitated these different evolutionary paths [92].

It is tempting to attribute the higher intraspecific variation in H. erectus (see also [26,93]) to the greater time depth of the H. erectus sample, but empirical evidence suggests that this had a modest effect on the total intraspecific variation, even on large timescales [94]. Higher cranial variation in H. erectus compared to recent H. sapiens implies correspondingly higher genetic diversity given the strong association of cranial and neutral genetic variation across hominoids (also measured by SEV in that study [54]). Low genetic variation in recent non-African H. sapiens is attributed to a major population bottleneck during the recent human diaspora out of Africa [4,11], so higher genetic variation in H. erectus implies that this species did not experience the same drastic population bottleneck during its migration.

Higher population differentiation variation in H. erectus is consistent with a hypothesis of neutral evolution for both species due to the greater time depth in H. erectus (Prediction 2), but is insufficiently nuanced to rule out other microevolutionary processes, such as reduced gene flow among populations or stronger local adaptation H. erectus [15,23,25–30]. The tests of neutral evolution provide some clarity. Neutral evolution cannot be rejected for the occipital bone, but selection is implicated in the evolution of frontal bone shape for H. erectus s.l. (but not when H. erectus is restricted to only the Asian lineages). Recent H. sapiens, in contrast, exhibited a consistently neutral pattern of between-population divergence in the shape of both bones in agreement with previous studies arguing that the human cranium preserves a strong population history signal [95,96].

These results point in potentially interesting directions. For example, given the strong relationship between brain morphogenesis and neurocranial form, the selection on the brain could translate to changes in the frontal bone. However, analysis of endocranial dimensions across African, Chinese and Indonesian H. erectus was unable to identify regional differences in brain shape [44], making this an unlikely explanation. Directional selection in modern human cranial shape, and particularly facial shape, is typically linked to very cold climates or a change in diet [96–101]. Frontal bone morphology, and particularly the form of the browridge, reflects the integration of the neural and bony structures of the upper face and anterior neurocranium [102,103]. Thus, local selection in H. erectus, perhaps related to climatic adaptation in the face, could have produced both the overall higher variation, greater population differentiation and the signal of selection in the analyses presented here. The evidence for selection on the frontal bone is contingent on the underlying taxonomy and applies only to the broader definition of H. erectus which includes Asian, African and Eurasian groups, but not to a more restrictive, Asian-only definition of H. erectus.

5. Conclusion

Homo erectus and H. sapiens share important population history characteristics, most notably a major dispersal across the Old World which probably involved serial founder effects, local adaptation and variable gene flow. Yet, quantitative genetic tests reveal distinct evolutionary histories for these species. The magnitude of cranial variation among populations is higher in H. erectus than recent H. sapiens occupying a similar geographical range. This may imply a more prominent role of natural selection in H. erectus. Indeed, tests of neutral evolution revealed a role for natural selection in frontal bone diversification among H. erectus, but not H. sapiens populations. This contrasts with the primarily neutral signal revealed for the occipital bone, implying distinct evolutionary trajectories for two bones. This pattern may have been facilitated by modularity between the face and braincase such that frontal bone evolution is a correlated response to selection on the face via integration of the anterior frontal bone and face.

Data accessibility

Data are available from the Dryad Digital Repository: https://dx.doi.org/10.5061/dryad.5qfttdz32 [104].

Competing interests

I declare I have no competing interests.

Funding

I received no funding for this study.

Acknowledgements

I thank Monica Castro for help with data processing. I appreciate the various institutions and individuals that provided access to fossil and comparative materials, including Bandung Institute of Technology, Geological Museum (Bandung), Gada Madjah University, Lembaga Ilmu Pengetahuan Indonesia (LIPI), National Museum of Kenya, Senckenberg Museum, American Museum of Natural History, Natural History Museum (London), Duckworth Laboratory (Cambridge University), Iwan Kurniawan, Yahdi Zaim and Yousuke Kaifu.