Respiratory adaptation to climate in Upper Palaeolithic modern humans: the case of Sungir and Mladeč.

When our human ancestors migrated out of Africa, they faced a considerably harsher climate, but the extent to which human cranial morphology has adapted to climate is still debated. In particular, it remains unclear when such facial adaptations arose in human populations. Here, we explore climate-associated features of face shape in a worldwide modern human sample using 3D geometric morphometrics and a novel application of reduced rank regression. Based on these data, we assess climate adaptations in two cruci al Upper Palaeolithic human fossils, Sungir and Mladeč, both of which have been archaeologically associated with a boreal-to-temperate climate. We found several aspects of face shape, especially the relative dimensions of the external and internal nose as well as of the maxillary sinuses, that are strongly associated with temperature and humidity, even after accounting for auto-correlation due to geographical proximity of populations. For these features, both fossils revealed adaptations to a dry environment, with Sungir being strongly associated with cold and Mladeč with warm -to-hot temperatures. As both fossils are dated among the earliest recent modern humans in Europe, our results suggest a relatively fast rate of climate adaptation in human respiratory morphology.


Introduction
The presence and the nature of climate adaptation in modern humans is a highly debated question, and not much is known about the speed with which these adaptations are obtained by our species.
Previous studies demonstrated that the facial morphology of recent modern human groups has likely been influenced by adaptation to cold and dry climates [1][2][3][4][5][6][7][8][9] . Although the age and rate of such adaptations have not been assessed, several lines of evidence indicate that early modern humans faced variable and sometimes harsh environments of the Marine Isotope Stage 3 (MIS3) as they settled in Europe 40,000 years BC 10 . In the present study, we explore whether the facial morphology of two of the earliest Upper Palaeolithic humans from Europe, Mladeč-1 (Czech Republic) and Sungir-1 (Vladimir, Russia), demonstrate signatures of adaptation similar to that of boreal-totemperate-adapted modern human groups, as suggested by their archaeological context.

Modern population evidence of craniofacial diversity
Current patterns of modern human craniofacial diversity evolved by a combination of both neutral and selective processes. Connections between craniofacial morphology, climate, as well as genetic and geographical similarities have long been investigated by various multivariate statistical approaches [1][2][3][4][5][6][7][8][11][12][13] . Several of these studies indicated that the neurocranium and upper face partly track neutral genetic and molecular distances, while midfacial shape also reflects climate 3,11 .
However, Betti et al. 5 claimed that neutral processes were more important than climate in shaping the overall human cranium. In their point of view, a large proportion of the signal for climate-related natural selection is due to the inclusion of populations from extremely cold regions (i.e., Arctic samples). Along these lines, von Cramon-Taubadell 14 reported that all of the cranial "one-function" modules are highly correlated with genetic distances between populations. Finally, and contrary to Harvati and Weaver 3 and Smith 11 , Reyes-Centeno et al. 12 suggested that while overall cranial shape was significantly correlated with neutral genetic distances, the temporal bone and the face exhibited stronger associations with genetic distances than the neurocranium.
Other studies pointed to strong climatic influences on midfacial structure, particularly those of the internal nasal region. Evteev et al. 15 reported a strong association of mid-facial morphology in Asian groups with climatic variables that contrast the temperate climate of East Asians and the very cold-dry climate of North Asians. However, when genetic distances between groups were accounted for, these correlations appeared less significant. In line with Evteev et al. 15 , Maddux et al. 8 demonstrated that the internal nasal cavity exhibited an ecogeographic distribution consistent with climatic adaptation, with crania from colder and/or drier environments displaying internal nasal cavities that are longer, taller, and narrower (especially superiorly) compared to those from hotter and more humid environments. Clinical and experimental studies suggest that tall, narrow noses assist in thermoregulatory processes to warm and humidify the inspired air to protect lung tissues from desiccation see 8 and references therein . Furthermore, cold-dry adapted crania tend to display tall zygomatic-maxillary interfaces, with concomitantly taller and wider maxillary sinuses 16,17 . No strong correlations between maxillary sinus size and climatic conditions have been shown, but current literature suggests that sinuses act as zones of accommodation for ontogenetic and phylogenetic changes of surrounding craniofacial structures, inclusive of the nasal cavity [17][18][19][20] .
Overall, the discussion about the significance that climate might have played in modern human variation and adaptation is still unresolved, partly owing to methodological challenges and misconceptions. For instance, correlations between neutral genetic distances and measures of morphological dissimilarity (morphometric distances) not necessarily indicate a neutral mode of evolution, as is often assumed in the literature reviewed above. Both neutral genetic markers as well as climate are more similar in geographically adjacent populations than in more distant populations; climate adaptation would thus also lead to a correlation between neutral genetic and morphometric distances 21 . Moreover, multivariate morphometric distances pool all measured variables, which may have very different functional and evolutionary roles, and thus lack statistical power to identify signatures of adaptation. Nevertheless, it is becoming clear that the mid-face (particularly nasal morphology) is more strongly associated with climatic variables than the neurocranium.
With this in mind, it can be expected that changes in nasal morphology were particularly important in helping our tropical-based ancestors settle in the colder and drier environments of Europe during the Late Pleistocene. Indeed, several recent empirical studies focusing on cranial morphology among early Europeans and Northeast Asian populations suggest that Neanderthal nasal morphology was compatible with the need for warming and humidifying inhaled air 22,23 .
Evteev et al. 7 also showed that several Upper Palaeolithic Europeans were similar to modern North-East European groups in terms of mid-facial morphology; however, these researchers did not extend their study to directly compare the Upper Palaeolithic European to modern North Asian samples, which may be more reflective of cold-adapted morphologies 6 . Furthermore, it is unclear as to when these nasal adaptations developed in the modern human lineage. Two individuals that can shed light on this question are the Mladeč-1 and Sungir-1 specimens, which represent two well-preserved crania of early Upper Paleolithic individuals living in relatively harsh environments.

Mladeč and Sungir: geographical, temporal, and climatic context
The Mladeč and Sungir sites are well known owing to the numerous individuals whose remains were buried or otherwise deposited at these sites, suggesting a prolonged rather than a fleeting occupation. First excavated in 1955, the Sungir site (56°10'30"N, 40°30'30"E) is located near the city of Vladimir, about 192km northeast of Moscow, Russia 24 ; excavations at this site have unveiled several graves and at least nine individuals 25,26 . The Sungir remains have been recently dated to around 35,000 14 C years BP 25 . Sungir-1, an almost complete skeleton of an older male, was unearthed in 1964 from Grave 1 with numerous grave goods including ochre, ivory beads, and other types of body ornamentation 24,26 .
The Mladeč caves are located in central Bohemian massif, Olomouc Krai of the Czech Republic, and yielded remains of up to 10 individuals 27 . This site was largely excavated under the direction of Szombathy in the late 1880s; excavations continued under various researchers through the mid-1900s. The Mladeč-1 cranium was unearthed from the main cave in 1881 28 . Note that earlier reports assigned the Mladeč-1 cranium as male but it could also belong to a younger female 29 . Direct radiocarbon dates for Mladeč human material cluster around 31,000 14 C uncalibrated years before present 27 .
Roughly contemporaneous, several lines of evidence suggest that the Mladeč and Sungir inhabitants may have also lived in broadly similar climatic conditions. During MIS 3, the neighbouring Mladeč Moravian area was partly covered by woodland dominated by conifers with the accompaniment of some deciduous trees, including oak, beech, and yew 30,31 . At the same time, geological studies of the frost features and molluscs biodiversity suggests a cold subarctic tundra 32 , thus pointing to a variable and changing environment that ranged from the steppe, shrub and steppe, and partially forested landscape 31 . This picture is repeated at the Sungir site, where soil and pollen analysis points to the Bryansk interstadial (32-24 kyBP), a slightly warmer period in the MIS 3, where boreal vegetation was dominated by spruce and pine with the admixture of birch [33][34][35] . Thus, both individuals likely lived in environmental conditions that could be more continental than at present, with sharper changes between warm and cold periods, similar to those of the modern Southern and Eastern Siberian regions [36][37][38] . These environmental conditions offer a unique opportunity to see how relatively early modern humans may have started adapting to colder-drier conditions in terms of craniofacial anatomy.
As such, several studies have assessed how the morphology of the Sungir specimens compares to other peri-glacial hominins, and whether morphological changes indicative of climatic adaptation are evident 7,26,39-41 . Bunak 39 found that Sungir-1 is morphologically most similar to early central Europeans, including Predmostí 3 and other Cro-Magnon-like specimens. In a more recent study, Evteev et al. 7 also found close similarities between Sungir-1 and Late Holocene modern Europeans. While informative, these studies may be considered limited in that they focus on external cranial features, without including internal features of the nasal cavity (such as internal nasal height and breadth) that are more functionally related to climate 8 .

The present study
The purpose of this study is to investigate whether two Upper Palaeolithic individuals, Mladeč-1 and Sungir-1 (hereby referred to simply as Mladeč and Sungir), display mid-facial and respiratory features characteristic of modern groups adapted to the boreal-to-temporal climate, as suggested by the archaeological evidence. Given their close geographical location and temporal proximity, we also expect that the two fossils will have equal manifestation of these climate adaptations.
The majority of methods used in previous studies pooled different cranial features and did not take into account the autocorrelation of morphological and climate variables between human groups due to their geographical proximity. Here, we rectify this problem by introducing a reduced rank regression approach based on generalised regression coefficients weighed by the geographical distances. The data has been collected as three-dimensional landmark coordinates of the external and internal mid-facial morphology from CT-scans, inclusive of 233 Holocene human crania from diverse ecogeographic regions, and analysed by a geometric morphometric approach (see Methods).

Reduced Rank Regression of the Shape of the Mid-Face
We analysed the shape and size of the mid-face separately. The relationship between group average facial shapes, as described by the shape coordinates of the 61 three-dimensional landmarks, and four climate parameters was explored by reduced rank regression (RRR). Similar to the more familiar two-block partial least squares analysis (PLS), RRR represents the multivariate association between two blocks of variables (here, shape and climate) by several pairs of linear combinations, or latent variables, of the measured variables. Unlike PLS, RRR distinguishes between a dependent and an independent block and maximizes the regression slope of one linear combination on the other, i.e., the effect of one unit change of the independent variables on the dependent variables, regardless of the occuring variance. Hence, RRR also captures effects of climatic properties on facial shape even if they vary less than other aspects. In our case, the regression slopes were also corrected for the autocorrelation due to geographic proximity (see Methods for more details).
The climate variables explained 54.5% of the variance in facial shape among recent modern human groups (corrected for geographic distances; see Methods). Of this accounted variation, over 94% was summarized by the first two dimensions of the RRR (Supplementary Tab. S1). The first dimension had a strong negative loading on both temperature parameters and on maximum humidity. The second dimension is mainly loaded on minimum humidity. In other words, the first dimension described how a cold and dry climate was associated with facial shape, and the second dimension described how the humidity during the driest months was associated with the shape variables ( Figure 1).

Figure 1.
Climate loadings for the first two dimensions of reduced rank regression (RRR). The first RRR dimension accounts for 73% of the total association between the four climate variables and facial shape (i.e., 73% of the summed squared partial regression coefficients), and the second dimension 21%. TMN -average minimum temperature of the coldest month; TMX -average maximum temperature of the hottest month; HMN -average minimum humidity (vapour pressure) of the least humid month; HMX -average maximum humidity (vapour pressure) of the most humid month. All climate variables were transformed to Z-scores.

Dimension 1 Dimension 2
As in other geometric morphometric contexts, the loadings of the shape vectors can be visualized as shape deformations ( Figure 2). The first dimension corresponded to the relative width of the external and internal nose and the relative size of the sinuses. In other words, populations living in warmer and more humid climates tend to have relatively wider and more prominent external noses, wider and lower internal nasal structures and overall smaller sinuses, all in combination with a slightly narrower and lower mid-face. A colder and drier climate, by contrast, was associated with an overall taller mid-face, narrower and flatter external nose, tall internal nose and large sinuses ( Figure   2). Along the second dimension, populations from drier climates, as reflected by minimum humidity, tended to have taller faces with "high cheek-bones", a narrow but prominent external nose, and a narrow internal nose with slightly more posterior superior ethmoidale points. In dry climates, the maxillary sinuses are shifted anteriorly compared with their position in more humid climates

Discussion
Our results demonstrate that multiple aspects of internal and external mid-facial shape are associated with climate, even when accounting for spatial autocorrelation among groups.
Temperature and minimum humidity were the main factors in driving climate-related variation in the shape of the upper respiratory tract. Groups from colder and drier climates tended to have an overall taller mid-face, a narrower and flatter external nose, a tall internal nose and large sinuses.
They also tended to have a larger absolute size of the face and the upper respiratory tract. Groups from warmer regions, irrespective of their geographical origin, had a relatively wider and more prominent external nose, wider and lower internal nasal structures and overall smaller sinuses, all in combination with a slightly narrower and lower mid-face. The association of the mid-facial and respiratory morphology with minimum humidity was not as strong as that with temperature and maximum humidity.
Our results correspond with and develop on the previous accounts of climate-related trends. In particular, previous studies suggest that the relative width of the external nose decreases in colder climates, whereas nasal width increases and the height of the nasal bridge decreases with increasing maximum humidity [42][43][44] . The morphology of the internal nose has been previously shown to correlate with humidity and to reflect environmental adaptation better than other parts of the respiratory system 8,45 , in agreement with the clinical and experimental literature 46,47 . Our study confirmed that populations in colder climates tend to have a relatively taller and narrower internal nose (Fig. 2). The dryness of the climate, irrespective of the temperatures, was reflected in a narrower superior part of the internal nasal cavity with a posterior positioning of the left and right superior ethmoidale points also see 48 . Furthermore, the absolute, but not relative, size of the internal nose separated groups from drier and more humid climates, in agreement with results by Maddux et al 8 .
The absolute size of the maxillary sinuses did not clearly relate to climate (Fig. 4), which is in agreement with the observations by Buck 49 and Butaric 50 . However, the relative sinus dimensions as represented by the RRR shape score were larger in groups from cold climates (Fig. 2: Dimension 1).. 17 suggested that the form of the internal nasal cavity and sinuses are largely determined by the shape and size of the maxillary and zygomatic bones. Hence, future studies on maxillary sinus morphology should consider this structure in conjunction with the surrounding internal and external anatomy.

Maddux and Butaric
Among the four climate variables, minimum temperature had the strongest effect on the absolute size of the mid-face and the respiratory organs (Fig. 4, Supplementary Tab. S2). In contrast to Noback et al. 1 , we found that -unlike for overall face shape -the relative sizes of the separate respiratory organs did not consistently relate to climate, (which is in disagreement with 6 . This inconsistency between absolute and relative sizes may owe to differential constraints imposed by heat retention versus those by the demands for warming and humidification of the inhaled air 46,47,52 .

Climate-associated features in Sungir & Mladeč midfacial anatomy
Some previous studies aligned Mladeč and Sungir by facial morphology with European groups 7,39,48 , but this is not supported by our analyses. This discrepancy is due likely to the different statistical and morphometric methods. Morphological similarities among populations are shaped by many developmental and evolutionary factors, including evolutionary drift and population history, which often dominate morphometric distances and ordinations (such as the CVA in Fig. S1). To circumvent this problem, we first identified the features of face shape that are strongest affected by climate in the worldwide sample, and then calculated individual scores for these "adaptive" features ( Fig. 3). In other words, we quantified morphological similarities of Mladeč and Sungir only for the supposedly climate-adapted features, not overall shape. At the same time, these multivariate scores are expected to trace adaptive signatures more reliably than univariate traits because climate adaptation likely affects each single trait only little, but summed over multiple traits that are all affected by climate the adaptive signature should stand out 53 . Additionally, spatial autocorrelation can induce spurious associations between morphology and climate but, on the other hand, climate is geographically patterned. As we used a weighted least squares approach to correct for spatial autocorrelation, the presented results are very unlikely to reflect geographical patterns or population history only.
Based on this approach, Sungir showed adaptation to cold due to a taller mid-face, narrower and relatively flatter external nose, tall internal nose and relatively large maxillary sinuses and choanal regions (Dimension 1 , Fig 2 and 3). The absolute size of the mid-face and respiratory features, aligned Sungir with Mongolian-Buryats, a group from a colder and drier climate, which is in agreement with studies describing this fossil's morphology as '…currently found only between Siberian mongoloids, Inuits and some American Indians' 54:148 .
However, it is important to note that thermoregulatory processes are not the sole function of the nose. Nasal size-particularly in terms of internal height dimensions and the choanal region-is largely driven by oxygen demand 8,55-57 . As metabolic rates increase in colder climates 23,58,59 , concomitant increases in oxygen demand may result in the relatively larger size of the internal nose and choanae 23,56,58,59 .
Previous studies note that the Sungir individuals, particularly Sungir-1, were relatively tall and large Mladeč showed features of adaptation to a warmer environment along the first RRR dimension (relatively wide internal and external nose, relatively small maxillary sinuses; Fig. 2 and 3).
Additionally, Mladeč showed adaptation to a dry climate with "high cheek-bones," a prominent external nose, and an internal nose with a slightly more posterior position of superior ethmoidale points, similar to African groups from dry climates and to Inuits (RRR dimension 2, Fig. 2 and 3). The such as a woolly mammoth, to wolverines and ground squirrels who prefer more wooded environments 30,63,64 . In addition, there is evidence of several sharp climate oscillations during the middle-interglacial that have been traced throughout Western and Eastern Europe 65,66 . In Dolni Vestoniče, the climatic oscillations were evident from the palaeosoil stratigraphy and associated malacofauna 67 . In the Mid-Russian European Plain the evidence comes from the palaeosoils and associated vegetation 64 . During these oscillations, colder and drier tundra conditions were followed by warmer periods with slightly more humid environments that encouraged the development of broad-leaf tree cover and associated malacofauna 68,69 . The warmer time periods may account for the adaptation to a warmer climate in Mladeč. However, it is also possible that Sungir's ancestors had longer time in the periglacial environment and therefore more generations to adapt to the local climate than did the ancestors of Mladeč.
The association of Sungir with modern groups from cold-dry environments and of Mladeč with the groups from a warmer and dry climate regarding their functional nose shape, irrespective of the overall morphological similarity between the two fossils ( Supplementary Fig. 1), incites the question about the speed of facial and respiratory adaptation to the changing climate in the Upper Palaeolithic. To answer this question, further investigation with a wider set of modern and fossil humans are required-particularly to sort out morphological differences among cold-dry versus hotdry conditions and to analyze corollary effects, such as the relationship between metabolic rate, temperature, and upper respiratory morphology. Such studies will provide a better understanding of how our ancestors could have adapted to different environments, and how researchers can best infer this from the fossil record.

Materials and Methods
To assess both external and internal craniofacial morphology, we used computed  * Location coordinates are transferred into a decimal system.
Cranial and maxillary sinuses were digitally segmented and rendered in Amira 5.6 70 using previously described protocols 17,50,71 . Visual inspection of the fossil scans indicated the presence of foreign materials, including matrix in the nasal cavity, ethmoidal air spaces, and several paranasal sinuses. Processing these scans to remove matrix and segment the maxillary sinuses followed standard semi-automated and manual techniques see 72,73,74 for reviews on working with fossil specimens . Additional details regarding the processing of the fossil scans is provided elsewhere 71 . Supplementary Fig. S4 provides the external and internal views of the segmented maxillary sinuses for the Sungir and Mladeč specimens.
The rendered cranial and maxillary sinus models were imported into the freeware program 3D-Slicer   Table S3 for abbreviations and external facial landmarks not illustrated here. Weather stations were chosen based on proximity to samples. Using the Monthly field observations database, the climatic values were obtained following protocols of Maddux et al. 8 . Specific geographic and climatic variables for each modern human sample can be found in Table 1.
Geographical distances between groups were obtained following previously published protocols 17 and accounted for wayward points along the possible migration routes avoiding migration across large bodies of water and additional environmental barriers 77 .

Statistical analyses
Data preparation and statistical analyses were carried out in R 78 using Geomorph 79 , Morpho 80 and vcvComp 81 packages. The landmark configurations were registered by Procrustes superimposition 82 and symmetrised by averaging original and reflected sets of landmarks for each individual. Because of the unbalanced sex composition of groups, we corrected for sexual dimorphism by subtracting from each individual the sex-specific mean shape and than adding the overall mean following 83 . The same treatment was applied to correct for sexual dimorphism in centroid size values. Both fossils were assumed to be male. But we also repeated the analyses without correction, which led to qualitatively similar results.
For each specimen, including the two fossils, we computed the centroid size (CS) of the fullface landmark configuration along with the centroid sizes of the five anatomical regions (face, external nose, internal nose, choanae and sinuses). Additionally, we calculated relative size values for the five regions by dividing their CS through that of the full configuration of 61 landmarks.
Population-specific size distributions were compared in boxplots.
Reduced rank regression (RRR) is a method to decompose the multivariate dependence of one set of variables onto another set of variables into some pairs (or "dimensions") of linear combinations (latent variables) that best represent this statistical dependence [84][85][86] . Unlike the more familiar two-block partial least squares analysis (PLS), RRR is not symmetric because the pairs of linear combinations have maximal regression slope (not maximal covariance as in PLS). This difference matters if the independent variables are highly anisotropic (differ in variances and covariances), which can be the case even for variance-standardized variables. For instance, if a combination of independent variables (say, the difference between maximum and minimum humidity) varies considerably less than other combinations (sum of minimum and maximum humidity), it will contribute little to the net covariances and not show up in the first PLS dimension(s). By contrast, RRR maximizes the regression slope, i.e., the effect of one unit change of the independent variable on the dependent variable, regardless of the variances. E.g., if one unit change in the difference in maximum and minimum humidity had a stronger effect on face shape than one unit change in average humidity, the first dimension of RRR would represent this effect of the differences in max. and min. humidities, even if it varies little across the populations. In other words, we seek the features of face shape (linear combinations of shape variables) that are most responsive to climatic differences, which are not necessarily those features that vary most across populations. As in standard regression models, RRR can be corrected for spatial dependencies by a weighted least squares approach.
More formally, for a sample of size n, let X be an n x p matrix of the p independent variables (the four climate variables in our application) and Y be the n x q matrix of dependent variables (the 183 shape coordinates). The symmetric n x n weight matrix W represents non-independences among the cases resulting from geographic proximity. Multiple methods have been proposed to estimate or model these spatial autocorrelations (e.g., Dormann et al. 87 ). We simply subtracted the geographic distance between each pair of the populations from the maximum of these pairwise distances and divided the result by the same maximum distance. The entries of the resulting weight matrix thus ranged from 0 to 1, with 1s corresponding to pairs of individuals in the same location. Compute the p x q matrix of partial regression coefficients, B=(X'W -1 X) -1 X'W -1 Y, and decompose B into its singular vectors and singular values: UDV'=B.
The first column, u1, of the p x p matrix U is the first loading vector for the climate variables and the first column, v1, of the q x q matrix V the first loading vector for shape. The second pair of vectors are the loading vectors of the second dimension, and so on. The first singular value, d11, equals the regression slope of the linear combination Yv1 on Xu1, and similarly for the subsequent dimensions.