Introduction

To understand the diversity pattern of human cranial form, many studies have aimed to identify primary factors associated with natural selection in the evolutionary process1. The majority of research has attributed this diversity to environmental factors, particularly climate adaptation2,3,4,5,6,7,8 or masticatory stress caused by dietary conditions6,9,10,11,12,13,14, as well as neutral genetic mutation15,16,17,18,19,20,21,22,23,24.

In terms of climate factors, the globular shape of the spherical neurocranium is explained by adaptation to an extremely cold climate by Allen’s rule25, postulating that mammals save heat loss by minimizing body surface area against volume3,5,17,18,26. Some arguments adopt Bergmann’s rule27, which relates variations in cranial size to temperature4,6,17,26,28 and proposes that body size in cold regions increases to protect against heat loss. Alternatively, mechanical and/or functional effects due to masticatory stress on the cranial form have been postulated10,12,13,29 in light of dietary conditions caused by differences in cooking culture or subsistence, such as those between agriculturalists and hunter-gatherers. The general interpretation is that more release from masticatory stress reduces the volume of masticatory muscles and cranial robustness. Aside from climate adaptation and masticatory stress, some global studies have explained the primary factor of phenotypic diversity of cranial form by neutral evolution22,30,31,32,33. Another explanation for cranial diversity is attributed to the allometric pattern34,35,36,37. For instance, our previous study37 demonstrated an allometric trend in facial proportion correlating to cranial size, where larger crania were associated with narrower upper facial parts.

For these approaches, earlier studies were based on linear measurement data, while more recent studies have adopted spatial three-dimensional geometric morphometric landmark methods. With significant improvements in computer performance, we employed an advanced three-dimensional homologous modelling technique to analyze geographic variations in the human cranial form37. This method uses 3D homologous models created using template fitting technique. Through template fitting, all scan data are represented as mesh models with the same number of data points and the same topology, sometimes referred to as homologous models38. Using this 3D homologous modelling technique, we demonstrated or reconfirmed diversities in cranial form among 148 ethnic groups worldwide, including variations in overall size and specific regional features in facial and neurocranial components.

The large datasets of homologous 3D cranial surface models generated in our previous study are useful for assessing phenotypic morphological variations associated with environmental factors; however, we did not deal with data of external variables represented by climate or dietary conditions, which are highly likely to affect cranial morphology. Therefore, our current study aimed to explore the relationships between external factors and cranial morphology diversity, as quantitatively described through analysis of 3D homologous models.

In conducting this study, we considered the following issues that can arise when preparing datasets for such analyses. Previous studies exploring phenotypical adaptations to different environments, in terms of microevolutionary perspectives of human populations, have not accounted for the issue of global-scale migration. This occurred after anatomically modern Homo sapiens first settled regions out of Africa during the late Pleistocene. Throughout the subsequent Holocene population history, the Neolithic period was revolutionary, with agro-pastoralists from the Indo-European language family spreading from the Levant region to Europe and India about 8 kya. Genome data suggest that most present-day Europeans descended from these later arrivals, who mixed with the indigenous hunter-gatherer Eurasian populations already present, alongside farming dispersals from the Near East39,40,41,42,43,44,45.

Meanwhile, such agricultural-driven population dispersal in West Eurasia is compatible with that in East/Southeast Asia, dating back to at least 8 kya. Archaeological investigations indicate that rice and millet farming first emerged in the Yellow and Yangtze River areas of China, leading to diverse outcomes across East Eurasia, Southeast Asia, the Japanese archipelago, and the Pacific Islands, starting around 4 kya46,47,48,49,50,51. During this period, populations speaking Austroasiatic, Austronesian, and Altaic languages spread across the region. Our “two-layer” model in East Eurasia contributes to discussions on global-scale population dispersals and interactions52,53,54.

These scenarios are based on population history, which suggests that after Homo sapiens spread from Africa to non-African regions about 50 kya, the first migrants had not necessarily continued to settle regions where they initially occupied. Later Holocene agricultural-driven migratory populations, which may have replaced by or interbred with pre-existing hunter-gatherers, settled for several thousands of years at most.

In examining climate adaptations among early- to mid-Holocene agricultural-driven immigrants, it is crucial to consider the rapid development of human civilizations followed by the agrarian revolution. This development significantly reduced the selection pressure on physical and physiological adaptations to natural climate stresses55. Advances in civilization-associated conveniences, such as housing, clothing, heating, and cooling technologies, have improved survival rates without the need for physical adaptations to harsh climates. Consequently, the morphological characteristics of migrants may not fully reflect local climate conditions. When studying climate adaptations in human morphology, it is essential to account for the effects of migration after the agricultural revolution, which introduces potential biases in population samples. Despite this, few studies have addressed these effects. This study aimed to analyze the impact of bioclimate and masticatory factors on the regional diversity of cranial forms using 3D homologous models from worldwide ethnic groups, taking into account the large-scale migrations during the early- to mid-Holocene Neolithic period and beyond.

Materials and methods

Cranial specimens

This study included 405 male cranial specimens selected from 150 populations across nine geographic units (Table 1). We recorded the 3D surface data of the specimens using a high-resolution 3D scanner (EinScan Pro produced by Shining 3D Co. Ltd., Shenzhen, minimum resolution: 0.5 mm, https://www.shining3d.com/) and then generated mesh wire models consisting of approximately 300,000 vertices. The 3D scan data from most of these specimens were reused when they were collected in our previous study37. The current study added new data collected from the collection of the American Museum of National History and the University Museum of the University of Tokyo (Supplementary Table S1). These collections are permitted to use for biological anthropology through a process of each museum’s or institution’s internal approval ethical assessment. Many of these cranial specimens were selected from the cranial measurements database after Martin’s craniometric definitions provided by Tsunehiko Hanihara26,31,32. In the same way as in our former study, we selected representative adult male skulls from various ethnic groups worldwide, most of which were collected between the 18th and early 20th centuries. Although the exact ages of many skulls were unknown, we focused on mature specimens, approximately 20 years old or older, excluding those with advanced age-related changes, such as complete alveolar bone and tooth loss. Consequently, the age distribution within and between groups was random. We assumed that intergenerational differences in cranial form were negligible over the short term based on a reference study56. To identify the individuals representing each population, we calculated the Euclidean distance from the population mean to each individual based on 37 cranial measurements, selecting the specimens with the smallest distance from the mean.

For statistical comparisons, the samples were grouped according to the major geographical units, as shown in Table 1. The “Africa” group consisted of samples from the sub-Saharan regions only. The samples from North Africa were included in the “Middle East” group, together with the environmentally similar West Asians. The “Northeast Asia” group included individuals of non-European origin only, and the “America” group consisted of American Natives only. The samples from Micronesia were combined with the “Polynesian” group. Most of the groups were geographically native specimens, whereas some groups (listed in italics) in Africa, Northeast/Southeast Asia, and America were defined by ethnicity.

Table 1 Cranial specimens (males) used for this study and combination patterns of sample sets (Italic: ethnic group).
Full size table

Template fitting to generate a homologous model

This study used a homologous modelling by template fitting to generate statistically analyzable 3D models. The modelling concept and technical method were described in detail in our previous study37. Below, we briefly summarize the procedure used in this study. For this purpose. we used HBM-Rugle software, which was developed through a collaboration between Medic Engineering, Kyoto, and the Digital Human Research Center, National Institute of Advanced Industrial Science and Technology, Tokyo (http://www.rugle.co.jp/).

We used the same template mesh model of the cranium (4,485 vertices, 8,728 polygons) with 56 landmarks described in a head-oriented coordinate system as in the previous paper37. X-axis represents the transverse axis, the Y-axis represents the frontal axis, and the Z-axis represents the sagittal axis, respectively.

Using HBM-Rugle, the template model was fitted onto the target scan data using 56 landmarks. The deformed template mesh model was then subdivided by the loop subdivision surface algorithm to create a finer mesh model consisting of 17,709 vertices (34,928 polygons). After final fitting, the mean distance b the mean distance between the other data points of the homologous model and the target scan data was 0.322 mm.

Morphometric variables

To detect major factors in cranial size and shape, given the high intercorrelation among bone segments, we performed principal component analysis (PCA) using the 3D coordinates of all vertices (53,127 XYZ coordinates) across homologous cranial models. We applied PCA to both non-normalized and centroid size-normalized datasets. The PCA based on non-normalized data provided a clearer characterization of cranial forms by geographic regions and allowed for easier interpretation of components. This is because the components with a contribution rate greater than 1.0% were less fragmented and had larger eigenvalues compared to those from the normalized dataset. Therefore, we chose to use PCA based on non-normalized data. Additionally, size is a critical factor related to environmental conditions, as discussed in the introduction and in line with Bergman’s rule. The software used for PCA was HBS composed in HBM-Rugle.

Bioclimate variables

Bioclimate data (BIO) were cited from the WorldClim website (version 2.1, https://www.worldclim.org/)57 that is freely available for academic use. This database provides 19 variables for biologically useful data analysis of worldwide regions. The variables are the average for the years 1970–2000. Although the cranial specimens and modern climate data were not from the same generation, we assumed that regional climate conditions, such as the contrast between worm/cold or wet/dry, remained relatively parallel. Our study focused on 6 of the 19 variables related to temperature and precipitation for their straightforward interpretation in analyzing responses to cranial forms: BIO 1 = annual mean temperature; BIO 10 = mean temperature of the warmest quarter; BIO 11 = mean temperature of the coldest quarter; BIO 12 = annual precipitation; BIO 16 = precipitation of the wettest quarter; BIO 17 = precipitation of the driest quarter. We selected the geographic coordinates for these bioclimate data based on the central ___location (town) or capital city, if applicable, of each population sample’s settlement area. This selection was made using bioclimate maps generated by DIVA-GIS software (https://www.diva-gis.org/) (Supplementary Table S2).

Masticatory variable

Masticatory function is directly affected by nutritional conditions, including food resources and cooking techniques, which determine the need for chewing power. Examining the impact of dietary practices on cranial morphology poses challenges due to the difficulty in collecting quantified food data from worldwide ethnic groups. Nonetheless, it is known that dietary conditions affect the development of masticatory organs, notably the jaw, and peripheral bone segments, such as the zygomaxillary portion. The growth of masticatory muscles, crucial for cranial form, is closely linked to dietary factors and includes muscles such as the temporalis, masseter, and lateral and medial pterygoids. Among the bone segments associated with these mastication muscles, the temporal fossa serves as a direct indicator of temporalis muscle volume, which is essential for mandibular movement. Additionally, the temporalis muscles have been the focus of discussions on how extensional stress affects neurocranial form, as reviewed in relation to Allen’s rule in the introduction. However, estimating the development of other masticatory muscles from the skull alone, without the mandible, is challenging. For this reason, we focused exclusively on the temporalis muscles. To measure this, we calculated the area within the infratemporal fossa on a horizontal section, delineated by the zygomatic arch, the external lateral walls of the temporal and sphenoid bones, and the posterior external wall of the maxillary sinus, projected onto the Frankfurt horizontal plane (Fig. 1) using homologous cranial mesh models. The masticatory variable, referred to as the InfraTemPoral Fossa (ITPF), was determined by the summed area of the left and right sides. Although this variable is a partial component of the overall cranial form, we treated it as a pseudo-external variable, similar to a climatic variable.

Fig. 1
figure 1

An example of homologous fitting model of 3D-scanned crania, and the measured projected area of infratemporal fossa (ITPF) on the Frankfurt horizontal plane (data consist of the summed area from both sides).

Full size image

Correlation analysis

To assess the response of the cranial form to climate conditions and mastication factors, correlation coefficients between the scores of principal components (PCs) and external variables comprising the six bioclimate variables (BIO) and the masticatory variable (ITPF) were calculated using the Bell Curve for Excel (version 4.06) software. Typically, multiple regression analysis is used to identify the primary factors impacting a target variable; however, in this study, this method was deemed inappropriate due to the high multi-collinearity among the external variables, as demonstrated later.

Correlation coefficients between those variables were calculated on the basis of the following three combination patterns of the datasets (Table 1), considering the sample bias due to past population movement during the late Pleistocene and/or early- to mid-Holocene periods:

Pattern A: 150 ethnic groups (405 individuals) from all regional units.

Pattern B: 90 ethnic groups (236 individuals), excluding Europeans (EUR), South Asians (SAS), Southeast Asians (SEA), and Polynesians (POL). However, Dravidians, Veddas, and Andamans (SAS) and Negritos (SEA) were maintained in this data pool, while Japanese and Koreans (NEA) were excluded. Most of these excluded regional samples were assumed to be descendants of early- to mid-Holocene Neolithic agricultural migrants or those who interbred with preexisting hunter-gatherers.

Pattern C: 63 ethnic groups (160 individuals), excluding American Natives from the dataset used in Pattern B. These excluded samples were assumed to have experienced more complex environmental adaptations during the late Pleistocene, as they migrated from Northeast Asia across Beringia and entered North America, with some groups continuing through Central America into South America.

Through above computing correlation coefficients, in case that both the bioclimate (BIO) and masticatory variables (ITPF) are highly significantly covariate with PC scores, we calculate partial correlation coefficients, in order to evaluate the relationship with each pair of two variables after removing the effect of one variable. In order to assess the significance of these coefficients, upward p-values adjusted using the Holm–Bonferroni method58 for multiple testing of R-package were provided in addition to normal probabilities.

Results

The eigenvalues and contribution rates of PCA applied to the non-normalized datasets of the 405 homologous cranial models and PC scores computed for the 405 cranial samples are given in Supplementary Tables S3 and S4. Consequently, 14 primary PCs were detected with contribution rates greater than 1% of the total variance, with the cumulative proportion of variances being 83.09%. The results of PCA analyses based on Pattern B samples and Pattern C samples were very similar to those based on Pattern A samples. Therefore, results of Pattern A samples were used in the following analyses.

Next, we calculated the correlation coefficients between these PC scores and the six bioclimate variables (BIO) (Supplementary Tables S2) and area of the masticatory variable (ITPF) (Supplementary Tables S4). The top of Table 2 gives those correlation coefficients greater than 0.3 with a highly significance level under 1% through both normal and multiple probability tests, which were calculated using the dataset of Pattern A of Table 1. Consequently, at this correlation level threshold, only three PCs were extracted at a weak or moderate level: PC1, PC3, and PC4. These probability values are given in Supplementary Tables S5 together with the upward p-values obtained through multiple testing.

Likewise, as shown in Table 2 in the middle, correlation coefficients with BIO and ITPF were computed again using the Pattern B sample set, excluding assumable descendants from agricultural-driven Neolithic or later migrants. This consisted of samples from Africa (AFR), American Natives (AMR), Melanesia (MEL), Middle East (MID), most of Northeast Asia (NEA), and hunter-gatherers in South Asia (SAS) and Southeast Asia (SEA). PC1, PC3, and PC4 were again detected in this analysis at a higher level than those in Pattern A. Furthermore, PC6 and PC12 were detected as significantly but weakly correlated with external bioclimate or masticatory variables. The coefficients resulting from the dataset in Pattern C, which excludes American Natives (AMR) from the dataset in Pattern B, are provided at the bottom. It is noteworthy that PC2 emerged as a factor that moderately correlated with the bioclimate variables. Among the BIO variables, the mean temperature of the coldest quarter (BIO 11) tended to be moderately correlated with these PC scores. The masticatory variable (ITPF) showed a higher correlation with PC1 scores in any pattern, as compared with other component scores.

Table 2 Correlation coefficients between PCA scores and bioclimate variables (BIO) and absolute projected area of space in the infratemporal fossae on the Frankfurt horizontal plane (ITPF). This table indicates only correlations greater than 0.3 and 0.4 (bold cells: greatest values in BIO) at the significant levels greater than 0.01. Partial correlation coefficients with BIO (Bold) and ITPF are given respectively in case that PCA scores are correlated to both kinds of variables.
Full size table

Since these some PCs are correlating with both the climate and masticatory variables, this table gives partial correlation coefficients which were calculated without effect one variables in each other of pair. In any pattern, The PC1 is more strongly correlating to IPTF than BIO, the attribution to either PC3 and PC4 is slightly higher by BIO than ITPF.

Figures 2, 3, 4, 5 and 6 represent the regression lines given in the above three datasets. Green lines are based on the dataset of Pattern A, including all samples; blue lines are drawn in the Pattern B sample set (see box symbols), excluding assumable populations affected by early- to mid-Holocene agricultural migrants; and the red lines are generated on the basis of Pattern C, excluding American samples (box symbols, except green circles). On the top of each graph, color maps of cranial surface are illustrated, which are highly correlated vertices to each PC. On the opposite side of the graph, virtual shapes with PC scores of ± 3 SD are shown. These morphing images are provided in Supplementary Videos S1-S5.

Fig. 2
figure 2

The regression lines between PC1 scores and bioclimate and/or masticatory variables: the area of infratemporal fossa (ITPF) and the mean temperature of the coldest quarter (BIO 11). The color map of the cranial surface at the top of the graph indicates vertices highly correlated with PC1. Virtual shapes with PC1 scores of ± 3 SD are also shown. The scale is given as a green colored sphere with a diameter of 50 mm. Abbreviations: see Tables 1 and 2.

Full size image

Figure 2 visually confirms that PC1 represents overall cranial size due to the large part of the cranial surface correlating with XYZ coordinates (white) as well as the correlation with ITPF and BIO 11 (mean temperature of the coldest quarter). The larger crania were associated with larger ITPF, reflecting the size of the infratemporal fossa; however, many Northeast Asians and Americans are plotted above the regressing lines, indicating that those ITPF sizes were above the proportionate largeness. The correlation coefficients with BIO 11 tended to be higher in the datasets of Pattern B and Pattern C, indicating that cranial size increases in colder regions, such as Northeast Asia and some American Natives, but decreases in South Asia and a proportion of Africans and Melanesians.

Figure 3 shows the chart plotting of PC2 which is highly correlating with both the Z (sagittal)-axis (dark blue: cranial length) and X (transverse)-axis (red: cranial breadth), clearly exhibiting a variation in cranial length/breadth proportion. The highly significant response to climate variables was detected only in the dataset of Pattern C, from which agricultural migrants and American Natives were excluded. BIO 11 and BIO 16 (precipitation of the wettest quarter) were moderately correlated with the cranial proportion represented by PC2. The Africans and Melanesians settled warmer and wetter environmental regions, possessing longer neurocrania. In contrast, Northeast Asians, Americans, and Europeans settled in colder and dryer environments, tending to have shorter and more globular cranial shapes. In terms of ITPF, no significant correlations were detected in any dataset.

Fig. 3
figure 3

The regression lines between PC2 scores and bioclimate variables: the mean temperature of the coldest quarter (BIO 11) and mean precipitation of the wettest quarter (BIO 16). The color map of the cranial surface at the top of the graph indicates vertices highly correlated with PC2. Virtual shapes with PC2 scores of ± 3 SD are also shown. The scale is given as a green colored sphere with a diameter of 50 mm. Abbreviations: see Tables 1 and 2.

Full size image

Figure 4 depicts another obvious cranial variation represented by PC3, showing differences in cheek forms relating to ITPF size. The greatest contrast was seen in broad and high cheek forms, associated with a higher face (enlarged on the XY axes: yellow, and the Y axis: green). This feature was observed in cold-climate occupants, such as Northeast Asians and Americans, who had a larger infratemporal fossa. In contrast, warmer-climate populations, such as Middle East Asians and parts of the African and Melanesian groups, exhibited narrower cheeks directly associated with a smaller infratemporal fossa.

Fig. 4
figure 4

The regression lines between PC3 scores and bioclimate and/or masticatory variables: the area of infratemporal fossa (ITPF) and the mean temperature of the coldest quarter (BIO 11). The color map of the cranial surface at the top of the graph indicates vertices highly correlated with PC3. Virtual shapes with PC3 scores of ± 3 SD are also shown. The scale is given as a green colored sphere with a diameter of 50 mm. Abbreviations: see Tables 1 and 2.

Full size image

Additionally, PC4 (Fig. 5) was clearly related to cheek form and position, particularly concerning the zygomatic bones and maxilla. Northeast Asians and American Natives (depicted on the XZ axes: pink) exhibited a large infratemporal fossa (broad zygomatic arch, which is enlarged along the X axis: red), contributing to a forward projection of the face. In contrast, Middle East Asians and Europeans tended to have a zygomaxilla region that was more posteriorly oriented. This difference in cheek and maxilla positions resulted in flatter facial profiles in the former groups compared to the latter groups, who are characterized by noticeably protruding nasal bones. Furthermore, substantial differences were observed in the neurocranium form. Northeast Asians and Americans typically showed a more inclined forehead, whereas Middle East Asians and Europeans had a lower orientation of the occipital inion.

Fig. 5
figure 5

The regression lines between PC4 scores and bioclimate and/or masticatory variables: the area of infratemporal fossa (ITPF) and the annual mean temperature (BIO 1). The color map of the cranial surface at the top of the graph indicates vertices highly correlated with PC4. Virtual shapes with PC4 scores of ± 3 SD are also shown. The scale is given as a green colored sphere with a diameter of 50 mm. Abbreviations: see Tables 1 and 2.

Full size image

PC6 characterized the breadths in the frontal and temporal proportions, which significantly correlated with ITPF at a weak level (Supplementary Fig S2). As a whole, this variability seems to be independent of any specific pattern of region-assemblage.

PC12 (Fig. 6) was related to the breadth of the orbital and nasal (piriformis) apertures, which respond to warm/cold climatic conditions. In colder regions, such as those inhabited by Northeast Asians and some Americans, these apertures tended to be narrower than those of tropical-region occupants, represented by Africans. This variation in the tightness of orbital and nasal openings influenced the overall frontal to mid-facial profile, where the former groups exhibited a higher and narrower structure, contrasting with the lower and wider profiles of the latter groups.

Fig. 6
figure 6

The regression lines between PC12 scores and bioclimate variables: the mean temperature of the coldest quarter (BIO 11). The color map of the cranial surface at the top of the graph indicates vertices highly correlated with PC12. Virtual shapes with PC12 scores of ± 3 SD are also shown. The scale is given as a green colored sphere with a diameter of 50 mm. Abbreviations: see Tables 1 and 2.

Full size image

In addition to the individual-level scatter diagram of the PCA scores above, we presented scatter plots of group averages for broader comparisons (Supplementary Fig. S1 and scores in Table S6). These plots were based on the average cranial homologous models derived from vertex datasets of individual homologous models across 150 ethnic groups. They offer a clearer view of cranial affinities across major regional distributions, with less overlap compared to individual plotting, which facilitates the interpretation of the depicted patterns.

Discussion

Through our analysis, it is evident that masticatory and bioclimate factors interact intricately with cranial forms, making it challenging to find distinctive attribute in some segments. Additionally, the most affected bioclimate variable was the mean temperature of the coldest quarter (BIO 11), rather than the annual mean temperature or temperature of the warmest season. Another notable finding was the increase in the number of PCs factors with higher correlation coefficients with climate factors when sample bias affected by hypothesized replacement or admixture due to early- to mid-Holocene agricultural-driven demic population expansion and migration was removed, or when samples such as Americans who experienced multiple environmental adaptations through long-distance migrations during the late Pleistocene were excluded. From these perspectives, environmental influences on cranial forms will be discussed, focusing on several key points.

The response of the overall cranial size (represented in PC1) to climate highlights the contrast between smaller crania in tropical, warmer regions of South Asia and some parts of Africa, as well as among near-equatorial islanders in Southeast Asia and Melanesia, and the larger crania in cold regions of Northeast Asia and America. This corresponds to Bergmann’s or the Island rule59 based on ecogeographical conditions4,6,17,26,28,60. The former relates to temperature, while the latter depends on the available space in the niche and food resources. It is worth noting that when using a sample dataset selected without consideration of migration, the correlation with temperature is weaker than when using a sample set excluding such biases in population movement. In addition, in PC1 and other PCs, the most affected bioclimate variable was the mean temperature of the coldest quarter rather than that of the warmest season. When calculating the standard deviation of these temperatures (Supplementary Table S2), regional variance is considerably larger in the cold season (SD = 16.4) than in the warm season (SD = 6.5), indicating greater geographic differences in thermal stress during the cold season. This larger geographic fluctuation in temperature during the cold season may contribute to a higher covariance with the PC scores. Furthermore, the cold climate likely posed significant challenges to human survival under increased selective pressures compared to warm environments61,62.

PC2 assigns the length/breadth ratio of the cranial vault, which describes the variation from the proportionally elongated crania of Australo-Melanesians and Africans to the globular crania of other regions, including Northeast and Middle East Asians, Europeans, and American Natives. These features have been described in numerous previous somatometric and osteometric studies63,64,65. This characteristic, known as brachycephalization of non-Africans, has long been argued with various hypotheses. The primary hypothesis is linked to adaptation to cold climates by reducing body surface area, in accordance with the Allen’s rule17,18,26, assuming that a shorter, globular head minimizes the surface area for the same volume compared to a long, elliptical form. Another well-known hypothesis suggests a reduction of masticatory stress, represented by gracilization of the temporal muscles due to decreased chewing power from softening food, which reduces stress on the lateral sides of the head, thereby suppressing transverse growth of the neurocranium6,10,11,12,13,14. The size of the infratemporal fossa, which well correlates with overall cranial size within the same regional group. However, between different regional groups, the size of the infratemporal fossa can be proportionally and absolutely discordant even if cranial sizes are similar (PC1: Fig. 2). For example, the infratemporal fossa of Northeast Asians and American Natives was markedly larger than that of Europeans and Middle East Asians, despite similar overall cranial sizes and short, globular cranial forms (PC2: Fig. 3). This may explain the weak correlation between masticatory variables and the length/breadth proportion of the cranial vault. Furthermore, this long/short head component was not considerably correlated with other facial proportions. Therefore, our analysis does not support the hypothesis that cranial length/breadth proportion is substantially affected by the development of masticating organs.

Regarding correlation with climatic conditions, more PCs were highly correlated with BIO 11 when using the Pattern C subset data excluding American Natives and some Eurasian samples that hypothetically include early- to mid-Holocene agricultural-driven demic migrants from other regions. It is well known that the population history of American Natives assumes Northeast Asian origins, migrating across the Bering Strait land bridge (Beringia) during the late Pleistocene at the oldest, irrespective of multiple migrations66. This scenario suggests that the ancestors of American Natives, once settled in Northeast Asia, including the Siberian region, experienced selective pressures from a cold geo-environment to a certain extent before crossing Beringia and dispersed across a broad area of North and South America in different environments. Therefore, many American indigenous peoples, except those in northern high-latitude areas, experienced secondary or multiple adaptations (e.g., colder in Northeast Asia > warmer in Central America > colder South America) during their migration history. As a result, no significant correlation of neurocranium shape was detected with any climate variable when including such groups in the dataset due to the complexity of their population history. When using a dataset that excludes samples from multiple environments, our analysis supports interpretations based on Allen’s rule to explain short, rounded crania in cold-climate regions.

The trait combination of cheekbone orientation, nasal projection, and facial height/breadth profile represents another distinctive facial characteristic among regional phenotypic variability, described through PC3 and PC4 (Figs. 4 and 5). Each component showed significant covariation with both climate and masticatory variables. However, the challenge in interpreting these results lies in the fact that both variables themselves correlate almost equally with each other (for instance, the correlation coefficient between BIO 1 and ITPF is 0.448 in Pattern B; Supplementary Table S7, Fig. S3), making it difficult to ignore their significance level. Therefore, even if multiple regression analysis is adopted, it would be challenging to narrow down the primary factor based on these coefficient variables. Nevertheless, the response to external factors attributing cranial variation could be assessed by focusing on several aspects. The wide, protruding cheekbones of Northeast Asians and American Natives reflect their large infratemporal fossa, necessitating strong chewing power due to dietary conditions in cold-temperature niches. In particular, high-latitude ethnic groups in severely cold Siberia and the Arctic region show well-developed shoveling and double shoveling, related to cold climate adaptation and powerful biting67, such as frozen meet, and using teeth as a tool for tanning leather by increasing the structural strength of incisors68. This trait is genetically controlled by EDAR polymorphism which had experienced strong positive selection69. Additionally, the large teeth, palatine, and mandibular torus—which are bony protuberances on the palatine and lingual surfaces of the mandible of American Natives are well-known characteristics frequently present among Arctic populations63,68,70,71,72, corresponding to enhanced mastication function with large infratemporal muscles. Although the partial correlation coefficients given in Table 2 indicates slightly higher correlation with BIO than ITPF, considering these morpho-functional interpretations, dietary conditions in cold regions may significantly affect facial morphology else than climate conditions alone. On the other hand, their flat face, represented by a flat nose in forward cheeks, is often linked to Allen’s rule, in which animals reduce protruding parts like limbs and ears to conserve heat loss. The facial characteristic, along with a globular neurocranium, superficially fits Allen’s rule among Northeast Asians and some American Natives, including Arctic groups. However, as mentioned earlier, the forward cheek position also appears to be related to the large infratemporal fossa too. In this skeletal facial orientation, nasal bones seem buried in forward-projected cheek maxillary bones, enhancing facial profile flatness. This contrasts with Middle East Asians and Europeans, whose protruding nasal bones remain in front of posterior-set cheeks, leading to a smaller space for the infratemporal fossa surrounded by a zygomatic arch. Their reduced temporal muscle space and masticatory stress are associated with the emergence of agriculture14,29. The impact of dietary changes on masticatory stress due to cultural shifts has been well documented, as in the study of global mandibular shape variation73. The study explained that the transition to systematic food production led agricultural societies to consume softer diets, which in turn significantly influenced the form of the human mandible. Given this interpretation, although the facial flatness trait superficially appears to follow Allen’s rule because the covariances with temperature are equal or slightly greater than that with masticatory variables, enlarged infratemporal fossa responding to masticatory stress boosts facial flatness as another primary factor.

The facial profiles described are highly correlated with the calvaria outline, particularly the degree of frontal and occipital inclines. Europeans are linked to a reduction in frontal incline and a flat, narrow occipital bone with a lower inion. In contrast, Northeast Asians tend to have an inclined forehead with a lifted occipital region. As discussed in our previous study37, the mechanism remains entirely uninterpretable with respect to either geo-environmental or mastication factors. Various possible interactive covariation factors, other than climate conditions, such as cranial weight balance, development of brain lobes including cerebellum74, muscle development of the neck, etc., remain unanswered questions.

Finally, variations in the shape of orbital and nasal apparatuses were described by PC12 (Fig. 6) in response to temperature, reflecting wider mid-facial and frontal bones among warm-climate occupants compared to narrower ones among settlers in cold regions. In the virtual shapes illustrated in the scatter diagram of PC12, from which size factors were likely removed, the facial features reflected a proportionally wider, lower face versus a higher, longer face. The primary climate attribution of these characteristics is not easily interpretable. One possibility is that increased maxillary sinus volume leads to the high narrowness of the nasal cavity and apparatus, as stated by several studies8,18,26,75,76,77,78, addressing the functional attribution of accommodating nasal apparatus with heat loss protection in cold dry climates. Although this phenotypical appearance in cold climate regions does not contradict our results, further analysis, including internal nasal form, such as maxillary sinus shape and volume, will be crucial in proving this argument. It is worth noting that the narrow, high nasal apparatus correlates with the narrowness of the orbital opening and forehead. This finding cannot be interpreted solely from the functional adaptation of the nasal cavity. In fact, a similar pattern of facial profiles was characterized through the analysis of allometric trends in our previous study37. The detected allometric pattern indicates that as cranial size increases, the face becomes higher, while the forehead and center of the mid-face (including orbital and nasal openings but excluding the zygomatic arch) tend to be proportionally narrower, characterizing not only Northeast Asians and some American Natives but also other large cranial samples from other regions. There is a possibility that such allometric trends boosted climate adaptation in the nasal portion among occupants of high-latitude cold and dry regions.

Through our global-scale analysis of human cranial forms based on 3D homologous models, it is evident that signatures of climate environmental stress fluctuate according to sample sets. Ultimately, the covariation between cranial forms and climate tends to be higher when removing assumed migration bias represented by early- to mid-Holocene agricultural-driven migration and/or removing population samples that experienced multiple climate adaptations during the late Pleistocene. A comparison of climate adaptation between humans and non-human primates78 suggests that humans exhibit fewer cranial adaptations to climatic stresses than primates. This implies that technological development by humans has mitigated stress from natural environments, reducing selection pressure on physical characteristics caused by climate conditions. Considering these studies, human cranial morphology may be more affected by cultural and behavioral adaptations, such as dietary conditions, food resources, and cooking techniques, than by natural environmental stress alone. Our results correspond to Bergmann’s and Allen’s rules to explain climate adaptation in mammals to a certain extent. However, our analysis in consideration with agricultural migratory bias for sample set demonstrated that such climate stress for selective pressure on the cranial forms was weaker during only several millennium span in the Holocene, supporting the view of reduction of selective pressure on human’s physical characteristics caused by natural environmental stress due to the development of civilization. Besides, our cranial analysis demonstrated the complexity of responses to environmental factors, particularly regarding facial forms, necessitating comprehensive interpretations including masticatory stress and/or allometric trends.

The aim of this study was to explore the relationships between external factors and cranial morphometric diversity using PCA. As a result, some PCs (such as PC3 and PC4) showed significant correlations with both masticatory and climate factors; however, interpretation was complicated by the high inter-correlation between these external factors. To address this issue, we adopted partial correlation analysis to assess the primary influence on cranial form and concluded that both factors contributed almost equally. Although we used PCA to evaluate morphometric variables representing major morphological variations similar to our previous study37, it is possible that PCA may not capture all biologically meaningful variation driven by environmental adaptation. To resolve this issue, it may be useful to consider utilizing the multivariate shape data as the response variable and adopting partial least squares regression to address multicollinearity. Testing such alternative methods could help confirm our findings or provide a more comprehensive analysis.