Subarctic climate for the earliest Homo sapiens in Europe.

The expansion of Homo sapiens across Eurasia marked a major milestone in human evolution that would eventually lead to our species being found across every continent. Current models propose that these expansions occurred only during episodes of warm climate, based on age correlations between archaeological and climatic records. Here, we obtain direct evidence for the temperatures faced by some of these humans through the oxygen isotope analysis of faunal remains from Bacho Kiro Cave, Bulgaria, the earliest clear record of H. sapiens in Europe. The results indicate that humans ∼45,000 years ago experienced subarctic climates with far colder climatic conditions than previously suggested. This demonstrates that the early presence of H. sapiens in Europe was not contingent on warm climates. Our results necessitate the revision of key models of human expansion and highlight the need for a less deterministic role of climate in the study of our evolutionary history.

INTRODUCTION

Models of the expansion of Homo sapiens posit that their dispersals into Europe and central Asia during the Late Pleistocene largely occurred during millennial-scale warm phases [Greenland Interstadials (GIs)] of the Last Glacial period (–). According to these models, the early wave of this expansion into Europe from southwestern Asia followed the cold phase of Heinrich event 5 [Greenland stadial (GS) 13] or GS 12 that potentially triggered a Neanderthal (Homo neanderthalensis) depopulation (–). Such models are critical as they make key contributions toward understanding the processes whereby H. sapiens spread across diverse climate zones and replaced Neanderthals in a few millennia. However, rather than relying on paleoclimate evidence coming directly from the archaeological record itself, these modeled scenarios predominantly rest on correlating the chronometric ages of archaeological finds with climatic phases documented in archives, such as ice cores or speleothems. Here, we use direct evidence of temperatures faced by humans during the Initial Upper Paleolithic (IUP) of Bacho Kiro Cave. We show that, in contrast to existing models proposing that warm climates were necessary for range expansion, at least some dispersals or the early continued presence of H. sapiens in Europe occurred during cold conditions.

Bacho Kiro Cave, located near the town of Dryanovo in central northern Bulgaria (fig. S1), has yielded one of the richest archaeological records of early H. sapiens in Europe here associated with the IUP. The rich IUP deposits include H. sapiens fossil remains and are currently thought to represent one of the earliest occurrences of H. sapiens in Europe (, ). The IUP occupations uncovered during recent re-excavations now date to 45,040 to 43,280 cal B.P. (calibrated years before 1950; 95.4% probability) in Layer N1-I & I and likely begin as early as 45,990 cal B.P. (95.4% probability) in Layer N1-J () [14C dates recalibrated using IntCal20 (); see site background in supplementary text S1 and fig. S2]. Analyses of the mtDNA of the Bacho Kiro Cave H. sapiens remains show that the IUP population represented at the cave did not contribute to the genome of modern-day Europeans (), suggesting that the Bacho Kiro IUP humans were part of a local population extinction, as previously proposed for other early H. sapiens finds (, ). Here, we apply oxygen and strontium stable isotope analysis to equid (Equus ferus and Equus hydruntinus) and aurochs/bison (Bos/Bison) tooth enamel from Layers N1-I & I and N1-J to characterize local seasonal paleotemperatures experienced by H. sapiens that produced the IUP archaeological record. In addition, we generate comparative data from the underlying Layer N1-K that was deposited between 61 ± 6 ka (thousand years) ago and 51 ka B.P. (infinite radiocarbon age), and which is attributed to the Middle Paleolithic (MP; supplementary text S2 and table S1). Strontium isotope analyses are used to confirm a lack of long distance migratory behavior and, therefore, suitability for local climatic reconstruction for the analyzed animals. Our sample combines teeth recovered in the 2015–2019 excavations in the Niche 1 sector (marked by “N1” prefix) and in the Main sector (supplementary text S1). IUP layers from these two sectors are clearly correlated (supplementary text S1), so we treat them as the same archaeological unit and use the designation of N1-I & I to denote the combined samples. We also use material recovered at the contact between two layers, such as N1-I/J. As the IUP faunal record in Layer N1-I & I is predominantly accumulated through human activity (supplementary text S3 and fig. S3), paleotemperature estimates generated by stable isotope analysis of faunal tooth enamel are directly representative of climatic conditions during human presence at the site (see methodological background in supplementary texts S4 to S8).

The fauna from the IUP layers comprises a mixture of taxa that range from temperate forest–adapted species to cold temperature–adapted taxa and species that can thrive in a large range of climates (). The predominant taxa are cave bear, cervids (especially Cervus elaphus), Bos/Bison, caprines (especially Capra ibex), and equids (). Layer N1-J & J also contains a few specimens of indisputably cold-adapted taxa, such as woolly mammoth, reindeer, giant deer, and wolverine, whose presence is very unusual for marine isotope stage (MIS) 3 in southeastern Europe (, ). Preliminary results from the study of micromammals suggest the presence of open habitats with a cooler climate than today (, ). However, given the ecological flexibility of many—especially larger—animals and the role of prey choice in the accumulation of faunal remains in the cave, more quantitative and high-resolution data on paleoclimatic conditions experienced by humans at Bacho Kiro Cave are needed.

RESULTS

Sequential oxygen isotope measurements of enamel bioapatite phosphate (δ18Ophos) form complete or partial sinusoidal curves with summer (high temperature) peaks and winter (low temperature) troughs for most studied teeth (Nteeth = 13; supplementary text S6 and fig. S4). Three Bos/Bison teeth (AA7-141, AA8-334, and AA7-2017) do not record a clearly visible peak or trough (as they are more heavily worn and preserve only a short period of isotopic input) and are excluded from further analysis. In contrast, a number of Equus sp. teeth record up to two complete annual cycles, enabling the extraction of several peak and trough data points. Last, 87Sr/86Sr values do not show substantial differences in the inhabited geology between summer and winter in any analyzed teeth and a low range of values with all measurements identical to the third decimal place (0.7090 to 0.7097, mean = 0.7093 ± 0.0002, 1 SD; fig. S5). This indicates that neither horses nor aurochs/bison were undertaking long distance migrations across different lithologies during this time. As no plant or bedrock 87Sr/86Sr baseline data are available for the vicinity of the site, determining the range of movement consistent with these data is difficult and reliant on inferences from similarity of lithology. Bacho Kiro Cave is located in an area of carbonate lithology, and similar lithology is present in a narrow band along the northern edge of the Balkan mountains (). This means that animals could potentially move relatively far within the same lithology in an east-west direction but only a few tens of kilometers in a north-south direction. However, such a scenario where at least three different animal species followed the same specific migratory pattern is much less parsimonious than a local origin within a few tens of kilometers of the site. In addition, east-west gradients of δ18Oprecip are very small, making such a migratory behavior extremely unlikely to substantially bias paleoclimatic reconstructions. We therefore conclude that δ18O values can be used to reconstruct local paleotemperatures without variability introduced by geographical changes in δ18O of precipitation (δ18Oprecip).

Seasonal δ18Ophos values are low compared to other Late Pleistocene faunal samples in all archaeological layers and show a small to moderate seasonal (summer to winter) amplitude in Layers N1-K, N1-I & I, and N1-I/J, with more pronounced seasonal δ18Ophos differences in Layer N1-J (Fig. 1). Overall, δ18Ophos values range from 11.3 to 16.4 ‰ in summer, from 10.9 to 14.9 ‰ for annual midpoint values, and from 10.0 to 13.8 ‰ in winter. There is a moderate shift to generally lower values from the Middle Paleolithic (meansummer = 14.1 ± 0.5 ‰ 1 SD, N = 4; meanmean annual = 13.6 ± 0.5 ‰ 1 SD, N = 5; meanwinter = 13.0 ± 0.6 ‰ 1 SD, N = 5) to the IUP as well as throughout the IUP from Layer N1-J to the top of Layer N1-I & I (Layer N1-J: meansummer = 14.2 ± 2.1 ‰ 1 SD, N = 3; meanmean annual = 13.4 ± 1.5 ‰ 1 SD, N = 4; meanwinter = 11.5 ± 1.1 ‰ 1 SD, N = 3; Layer N1-I & I: meansummer = 12.6 ± 0.7 ‰ 1 SD, N = 5; meanmean annual = 12.2 ± 0.8 ‰ 1 SD, N = 7; meanwinter = 11.5 ± 1.2 ‰ 1 SD, N = 6). Analysis of variance of δ18Ophos values shows significant differences between layers for summer, winter, and annual midpoint values, respectively (psummer = 0.015, pwinter = 0.022, pmean annual = 0.0003; normality and equivalence of variance confirmed by qqnorm plots and Levene’s test). A Tukey test reveals that significant differences can be detected between the Middle Paleolithic (Layer N1-K) and the IUP occupation of Layer N1-I & I for annual midpoint and winter δ18Ophos values (pmean annual = 0.00034, pwinter = 0.027). In addition, statistically significant differences can be seen between the lower IUP material (Layer N1-J) and the upper IUP layers (N1-I & I) for annual midpoint and summer δ18Ophos values (pmean annual = 0.018, psummer = 0.012). However, it should be noted that the sample size in Layer N1-J is small. At the same time, we observe on average a stronger seasonal amplitude (difference between summer and winter values) of 2.6 ± 1.2 ‰ in Layer N1-J at the start of the IUP than in either Layer N1-K or Layer N1-I & I, where the seasonal amplitude is smaller (Layer N1-K meanampl = 1.0 ± 0.7 ‰, 1 SD; Layer N1-I & I meanampl = 1.5 ± 0.8 ‰, 1 SD). However, these diachronic changes occur within a framework of generally low δ18Ophos values in all layers in comparison to many other studies of MIS 3 fauna [e.g., ()].

Fig. 1.

Seasonal oxygen isotope data for IUP and Middle Paleolithic (MP) layers of Bacho Kiro Cave.

Temperature-driven summer (pink), annual midpoint (gray), and winter (blue) oxygen isotope values (δ18Ophos) extracted from sequentially sampled Equus sp. (diamonds, right) and Bos/Bison (circles, left; only Layer N1-I & I) tooth enamel show overall low values, particularly in the IUP Layer N1-I & I. Some higher summer δ18Ophos values occurred at the beginning of the IUP in Layer N1-J and at the contact between N1-I and N1-J (N1-I/J), but the sample size for these layers is small. Plotted points represent δ18Ophos summer peaks, winter troughs, and annual means of individual years represented in sinusoidal δ18Ophos time series obtained from sequential samples taken along each analyzed tooth (see supplementary text S4). Summer peak and winter trough values were obtained by visual inspection of each δ18Ophos measurement series and are marked individually in fig. S4. Annual midpoint values represent the mean of the summer and winter values. A comparison of annual midpoint values and full annual averages can be found in supplementary text S8 and figs. S6 and S7. Means for summer, annual midpoint, and winter records, respectively, are connected by shaded lines, while shaded ribbons visualize the maximal spread of the data. For summer and winter values, error bars represent measurement uncertainty (1 SD) as determined by replicate measurements of each sample. For annual midpoint values, error bars represent the uncertainty around the mean derived by error propagation of the measurement uncertainty.

Last, reconstructed drinking water δ18O (δ18Odw) values of Bos/Bison and Equus sp. from Layer N1-I & I are systematically below modern-day estimates of δ18Oprecip for the area of the cave (fig. S8), underlining that the values obtained here are generally much lower than expected for a warmer phase climate, such as the modern-day climate of the Balkans. Bos/Bison δ18Odw values fall inside the range of variation of Equus sp. δ18Odw values within error, although they commonly fall on the upper end of the Equus sp. range (fig. S8). Considering the uncertainty introduced by converting to δ18Odw, the two taxa are in good agreement. Therefore, we combine the estimated δ18Odw values of both taxa to generate a more precise and robust temperature estimate for Layer N1-I & I (Fig. 2).

Fig. 2.

Reconstructed summer, winter, and mean annual paleotemperatures.

Air temperatures (°C) reconstructed for the MP (Layer N1-K, green) and IUP (Layer N1-J and following, blue) at Bacho Kiro Cave fall substantially below modern-day temperatures [horizontal lines (mean) and shaded areas (SD); 2012–2020 data obtained from Gorna Oryahovitsa ()] for summer (pink), mean annual (gray), and winter (blue) temperatures. In particular, the IUP occupation of Layer N1-I & I shows especially cold conditions with mean annual temperatures below freezing. Reconstructed temperatures are similarly low at the contact between N1-I and N1-J (N1-I/J); however, the small sample size for this stratigraphic unit means that this result is less secure. A small temperature decline can be seen from the MP to the IUP, but sample numbers are low for the lower IUP layers [noted in gray as number of data points (Ndata) and number of teeth (Nteeth)]. Plotted points represent calibrated temperatures reconstructed for each layer derived from Equus sp. (diamonds) or a combination of Equus sp. and Bos/Bison (triangles) oxygen isotope measurements. Error bars indicate compound error around each temperature reconstruction. Summer and winter temperatures were converted from maxima and minima of the inverse model described in supplementary text S7, while mean annual temperatures were converted from the annual midpoint of unmodeled oxygen isotope values, as described in supplementary text S8.

Paleotemperature estimates fall substantially below local modern-day conditions (Fig. 2). In Layer N1-I & I, reconstructed temperatures are on average approximately 14°C lower than today, with the largest difference in summer and the smallest in winter. In Layer N1-K, temperatures are approximately 11°C below modern-day conditions, again with the smallest difference in winter. Diachronic trends mirror those noted in δ18Ophos, but the error ranges of temperature estimates from different layers overlap owing to the uncertainty introduced by conversion to temperature estimates. It should be noted that seasonal temperature differences reconstructed here most likely represent minimum seasonal amplitudes because of the specifics of the inverse model that was used to remove the effect of tooth enamel mineralization time averaging on the δ18Ophos seasonal amplitude (see supplementary text S7). For this reason, summer temperatures may have been higher and winter temperatures may have been lower than estimated here.

DISCUSSION

Using oxygen isotope values of Equus sp. and Bos/Bison from the IUP, we provide evidence that at least a portion of the human use of Bacho Kiro Cave in the context of the IUP took place during pronounced cold conditions consistent with a GS. This contrasts with models that posit H. sapiens preferences for warm environments during their expansions into Europe in the Upper Pleistocene. On the basis of seasonally homogeneous 87Sr/86Sr values, it is unlikely that these oxygen isotope results are affected by any migratory behavior of the analyzed animals, which renders them a faithful proxy for paleotemperature in the area around the site.

While directly comparable δ18Ophos values for the Upper Pleistocene are sparse, values from Layer N1-I & I most closely resemble those coming from glacial phases (i.e., MIS 4 or GSs) from sites at higher latitudes than Bacho Kiro Cave. This is the case for all analyzed IUP and Middle Paleolithic deposits of the cave [and while some diachronic differences can be seen (Figs. 1 and 2), these took place in a generally cold framework]. For example, annual midpoint δ18Ophos values in Equus sp. from Bacho Kiro Cave Layer N1-I & I (meanmean annual = 12.2 ± 0.8 ‰ 1 SD) are similar to δ18Ophos annual midpoint values of 12.5 ‰ (n = 4; converted from carbonate oxygen isotope values, δ18Ocarb) and 12.3 ‰ (n = 1, converted from δ18Ocarb) of equids from glacial phase Weichselian sites Bocksteinhöhle and Vogelherdhöhle in southwest Germany (), as well as annual midpoint values of 12.4 and 13.2 ‰ obtained from equids recovered from the sites of Boncourt Grand Combe (n = 2; MIS 3) and Courtedoux-Va Tche Tcha (n = 11; MIS 5a) in Switzerland (). This indicates that climatic conditions of the Bacho Kiro Cave IUP are consistent with conditions of a GS. This is supported by comparisons of converted δ18Odw values of other faunal species from the Upper Pleistocene, where Layer N1-I & I results match well with those obtained from woolly mammoths from sites in the Baltics dating to between 47.5 and 42.5 ka (). Low δ18Ophos values in equids and Bos/Bison could theoretically also be caused by substantial consumption of drinking water from high altitudes, snow or glacier meltwater, or deep groundwater. However, we argue that the hydrological situation in the region does not support a scenario where local rivers had substantially lower δ18O values (supplementary text S4), and the seasonal variability seen in δ18Ophos time series suggests that animals did not regularly drink from strongly buffered water sources such as groundwater (supplementary text S4).

Comparing both reconstructed temperatures and δ18Odw values to modern-day temperatures and δ18Oprecip in Eurasia suggests that climatic conditions during the IUP occupations of Bacho Kiro Cave are most closely comparable to current conditions in Scandinavia and Russia (–) (see supplementary texts S4 and S8 for assumptions underlying paleotemperature reconstruction). Furthermore, these reconstructions demonstrate that the animals analyzed here—and the humans that hunted them—lived during conditions that are consistent with the more intense millennial-scale glacials observed in other climatic records from southeastern Europe (e.g., Heinrich events). Air temperatures of ∼10° to 15°C below modern-day conditions are thought to have occurred in the region during millennial-scale stadials, such as the Heinrich events, while GIs are commonly reconstructed much closer to modern-day conditions (–).

A tentative trend of decreasing temperatures and temperature seasonality (summer-winter temperature difference) can be observed from the earlier IUP in Layer N1-J (∼46 to 44 ka cal B.P.) toward the later IUP in Layer N1-I & I (∼45 to 43 ka cal B.P.; Fig. 2), but it should be noted that the sample size for Layer N1-J is small and that this trend occurs within the framework of generally low temperatures in all analyzed layers. While still consistent with a stadial phase, the climatic conditions during the earlier IUP at Bacho Kiro Cave show substantially higher summer temperatures paired with quite harsh winter conditions—more similar to continental climates found in Russia or central Asia today—while the summer-winter difference in the later IUP is much smaller, closer to modern-day northern Scandinavia. A limited amount of preliminary data on the season of death of equids, Bos/Bison, and ursids indicates that human presence at Bacho Kiro Cave during the IUP was not restricted to the warmer summer season (supplementary text S3) despite the low teeratures reconstructed here for the colder seasons. Given the limited nature of the season of death data, however, it is currently not possible to determine the frequency or intensity of site use between different seasons. Because of the anthropogenic nature of the IUP (Layer N1-I & I and upper N1-J) faunal record, the climatic conditions reconstructed from faunal remains are directly tied to incidences of human activity at the site for these phases. In contrast, this is not necessarily the case for the fauna from the Middle Paleolithic deposits (Layer N1-K), where rates of carnivore modifications are more substantial (supplementary text S3), making it less clear whether faunal remains, and our samples, from Layer N1-K were accumulated as prey of humans or carnivores. The stable isotope results from these layers, therefore, are not tied specifically to instances of human occupation of the cave. Instead, climatic data for these layers need to be interpreted as generally representative of the time of layer formation, which occurred between 61 ± 6 ka ago and >51 ka B.P. (infinite radiocarbon age).

Cold conditions during the IUP occupations of Bacho Kiro Cave could perhaps explain the unusual presence of woolly mammoth, reindeer, giant deer, and wolverine in the faunal record for that time period at the site. The presence of species adapted to more temperate conditions, such as red deer (C. elaphus) might, however, suggest that there was some climatic variability within the IUP occupation (). It does seem unlikely though that a substantial portion of this record formed during temperate phases, as no signals of milder climate are captured in the isotopic data in either horses or bison/aurochs—species with a large climatic tolerance that are commonly found in both warm and cold phases. Geographical heterogeneity in local climates and habitats of the site area may also explain the presence of a mix of species with different climatic preferences. Results of cold conditions during the IUP occupations of Bacho Kiro Cave are additionally in broad agreement with colder and more open habitats reconstructed from preliminary analysis of the micromammals of the same deposits ().

A stadial phase IUP occupation of Bacho Kiro Cave is in contrast to previously proposed climatic conditions of H. sapiens dispersals into Europe immediately before the Upper Paleolithic. Currently, both the IUP and the subsequent Aurignacian expansions of H. sapiens into higher latitudes are thought to have occurred during GIs and potentially after a (proposed) climatically driven decline of Neanderthal populations in Europe during Heinrich event 5 [(–, ), but see () and ()]. These models are mostly founded on the warm phase assignment of a number of archaeological deposits. The beginnings of the IUP in central Europe have been correlated to GI 13/14 at the site of Bohunice by comparing it to the chronometric dates of the North Greenland Ice Core Project (NGRIP) record (, ). Chronometric dates of the IUP in Central Asia suggest the appearance of the IUP in the Altai during GI 12 (, ), which is supported by pedogenesis indicators in the Tolbor 16 stratigraphic sequence (). To support such data, a general point about a preference of early H. sapiens in Europe for warmer environments has been made on the basis of comparisons of Aurignacian site densities with climate simulations (, ). Furthermore, correlations of archaeologically sterile layers with the GS 12 cold phase shown in speleothem records and the NGRIP have been used to propose that the expansion of the Aurignacian took place after a decimation of the Neanderthal population by this stadial (). Some studies, such as the investigations at Willendorf, have, in contrast, yielded indications of cold phase dispersal of H. sapiens speculatively tied to an expansion of steppe landscapes (), but these remain the exception in an overall “warm phase” model.

A cold phase IUP occupation of Bacho Kiro Cave, however, suggests a different and likely more varied relationship between climatic conditions and the peopling of higher latitudes by our species during the Upper Pleistocene, as well as the importance of nonclimatic factors. Given the lack of direct climatic evidence for other European IUP sites, it is challenging to robustly support alternative models such as a steppe corridor dispersal (). However, the newly generated Bacho Kiro Cave isotope data show that currently favored models need to be revised to account for potentially diverse dispersal mechanisms and cold stage presence of H. sapiens in Europe at the start of the Upper Paleolithic. Our results indicate that humans associated with the IUP were better adapted to harsher climate conditions and were more environmentally flexible than previously thought. In addition, a relatively early presence of H. sapiens in southeastern Europe during a GS suggests that their spread at this time was unlikely to have been dependent on an immediately preceding climatically induced decline of the Neanderthal population. As the earliest H. sapiens groups in Europe appear to have been present during a GS with a duration of only a few millennia, it appears unlikely that Neanderthal populations could have been substantially decimated within the same millennial-scale cold event but before the arrival of H. sapiens or would not have recovered if such a depopulation had taken place in an earlier stadial. However, the demographic and environmental mechanisms for the subsequent dispersal of H. sapiens, commonly related to the material record of the Early Upper Paleolithic (Protoaurignacian, Aurignacian), may be different from those behind the expansion during the time of the IUP.

At the same time, the contrast between our results and established models underlines the importance of obtaining climatic evidence directly from archaeological deposits to supplement age correlations with long-term climate archives. On the basis of comparisons of the radiocarbon chronology with speleothem records from Romania (, ), Layers N1-I & I and N1-J & J were originally considered to most likely fall into the warm phase around GI 12 (, , ). Similarly, a comparison with the Tenaghi Philippon pollen record would also support this assignment (). However, the direct evidence of cold climatic conditions presented here shows that an assignment of the IUP of Bacho Kiro to a GI is inappropriate. We believe that this apparent disagreement is caused by a combination of the large combined uncertainty of the chronometric dates of archaeological deposits and climatic records, spatial differences in the timing and expression of climatic events, and the different climatic and environmental features represented by each climatic proxy rather than methodological inconsistencies in any of the climatic reconstructions. Combined uncertainty of ice-core layer counting, 14C, and U-series dates often span several thousand years, easily covering several climatic oscillations that can additionally vary in timing or intensity across the globe. Moreover, the calendar ages obtained from radiocarbon dating are dependent on the calibration curve used, and recalibration of the Bacho Kiro Cave 14C dates, for instance, has shifted the dates for the IUP occupation by up to 950 years, creating a larger overlap with GS 12. The case of Bacho Kiro Cave again highlights that chronometric correlations between archaeological deposits and spatially distant climate records offer climatic context of a resolution that is not always appropriate for the questions asked. Such correlations should therefore be strengthened by direct climatic evidence where possible.

MATERIALS AND METHODS

Tooth sample selection and sequential sampling of tooth enamel

Oxygen isotope analyses were applied to sequential samples of Bos/Bison and Equus sp. tooth enamel from 13 individual teeth (nBos/Bison = 5, nEquus sp. = 8; Table 1) from the faunal collection of the 2015–2019 re-excavation of the site. Samples of Equus sp. and Bos/Bison teeth were taken from the IUP Layers N1-I & I and N1-J and the Middle Paleolithic Layer N1-K, with most samples coming from the Niche 1 sector. The IUP layers of the two excavation sectors (Main sector and Niche 1) can be clearly correlated with each other but are named separately with N1 prefixes denominating layers in the Niche 1 sector. Samples from Layers I and N1-I are treated as a combined unit for the purposes of this study and termed N1-I & I. Artifacts found at the contact between Layer H and Layer I or N1-H and N1-I (labeled as H/I or N1-H/I, respectively) appear to be reworked from the surface of Layer N1-I & I () and are therefore grouped with the Layer N1-I & I unit here. One sample was obtained from the contact of Layers N1-J and N1-I and is labeled as N1-I/J. All samples from Layer N1-J were obtained from the Niche 1 sector, with none coming from the Main sector, where the layer is designated as Layer J. Samples were only obtained from the upper parts of this stratigraphic unit, which exhibits radiocarbon dates ranging from 46 to 44.5 ka cal B.P. (95.4% probability; not including dates for Layer J; see supplementary text S1). This section of the layer is securely attributed to the IUP and technologically consistent with the overlying Layer N1-I & I. Artifacts are spare in the lower parts of Layer N1-J preserved in the Niche 1 sector, and the layer exhibits a gradual contact with the underlying Middle Paleolithic in Layer N1-K, as exemplified by radiocarbon dates >51 ka B.P. obtained from the lowest parts of Layer N1-J. No animal remains from these lower layer portions were included here, and all results generated for Layer N1-J are assigned to an early phase of the IUP of Bacho Kiro Cave. A series of sequential enamel samples was obtained from each tooth, covering the full length of the tooth. A total of 179 sequential samples were processed for this study. Sequential tooth enamel samples were obtained as strip slices from loph sections of each tooth (see full details in supplementary text S5). Tooth enamel sections were cleaned of adhering dentine and calculus and sectioned along the growth axis into a series of ∼1-mm-wide sequential sample strips. Usually, every third strip sample was further processed for isotope analysis, resulting in a series of sequential samples with ∼2-mm gaps between samples. Each strip that was chosen for isotope analysis was split in two parts, one part for oxygen stable isotope analysis and the other half reserved for strontium isotope analysis. From a subset of individuals, two sequential samples each representing the summer and winter season input were chosen to be analyzed for 87Sr/86Sr. Only individuals with both summer peaks and winter troughs were analyzed for 87Sr/86Sr to facilitate detection of seasonal migratory behavior by comparing summer and winter 87Sr/86Sr values.

Table 1.

Context and sample information for teeth sampled for oxygen and strontium stable isotope analysis.

The period of tooth mineralization is given as the age in months of the start of enamel mineralization to the age when enamel mineralization is complete. Tooth mineralization data are based on studies of modern horses () and bison () (who are very similar to cattle in tooth formation).

Findnumber	Layer	Taxon	Tooth position	No. of seq. samples	Mineralization period
AA7-141	N1-H/I	Bos/Bison	Left mandibular M3	9	9–24 months
AA7-121	N1-H/I	Equus ferus	Left maxillary M2	16	7–37 months
AA7-52	N1-H/I	Equus ferus	Left mandibular M3	22	21–55 months
AA7-2017	N1-I	Bos/Bison	Right mandibular M3	6	9–24 months
AA8-334	N1-I	Bos/Bison	Left maxillary M2	7	2–14 months
CC8-18	N1-I	Bos/Bison	Right mandibular M3	13	9–24 months
CC7-2813	N1-I	Equus ferus	Left maxillary P2	20	13–31 months
A7-534	I	Bos/Bison	Left mandibular M2	9	2–14 months
CC7-2397	N1-I/J	Equus hydruntinus	Right maxillary M3	13	21–55 months
CC7-2605	N1-J	Equus ferus	Left maxillary P4	12	19–51 months
CC7-2478	N1-J	Equus ferus	Right maxillary M3	18	21–55 months
CC8-2419	N1-K	Equus ferus	Right maxillary M3	17	21–55 months
CC7-3018	N1-K	Equus ferus	Right maxillary P3	17	14–36 months

Oxygen stable isotope analysis

Tooth enamel strip samples were ground to powder using a clean agate mortar and pestle. Approximately 5 mg of each powder sample was converted to silver phosphate for oxygen stable isotope analysis of bioapatite phosphate using digestion with hydrofluoric acid followed by crash precipitation of silver phosphates following an adapted version of the protocol developed in () and modified in () as described in () (see full details in supplementary text S5).

Oxygen isotope ratios of Ag3PO4 were analyzed using a High-Temperature Conversion Elemental Analyzer (TC/EA) coupled to a Delta V isotope ratio mass spectrometer via a Conflo IV interface (Thermo Fisher Scientific, Bremen, Germany) at the Max Planck Institute for Evolutionary Anthropology (MPI-EVA). Details of the analytical setup can be found in supplementary text S5. Oxygen isotope delta values were two-point scale normalized to the Vienna Standard Mean Ocean Water (VSMOW) scale using matrix matched standards calibrated to international reference materials, and scale normalization was checked using three separate quality control standards in each run. Scale normalization was conducted using the B2207 silver phosphate standard (δ18O = 21.7 ± 0.3 ‰, 1 SD; Elemental Microanalysis, Okehampton, UK) and an in-house silver phosphate standard (KDHP.N, δ18O = 4.2 ± 0.3 ‰, 1 SD). The accepted value of this in-house standard was determined by two-point calibration using B2207 and IAEA-SO-6 [barium sulfate, δ18O = −11.35 ± 0.3 ‰, 1 SD, as given in ()]. Aliquots of an in-house modern cow enamel standard (BRWE, later replaced by BRWE.2 owing to exhaustion of this material) and the standard material National Institute of Standards and Technology (NIST) Standard Reference Material (SRM) 120c [formerly National Bureau of Standards (NBS) 120c] were precipitated and measured alongside each batch of samples to ensure equal treatment. In addition, a commercially available silver phosphate (AS337382, Sigma-Aldrich, Munich, Germany) was used as a third quality control standard to check across-run consistency of scale normalization independent of silver phosphate precipitation. Measurements of these standards gave δ18O values of 14.9 ± 0.4 ‰ for BRWE (1 SD, n = 53), 21.7 ± 0.5 ‰ for NIST SRM 120c (1 SD, n = 37), and 13.9 ± 0.2 ‰ for AS337382 (1 SD, n = 159). This compares well to the consensus value for NIST SRM 120c of 21.7 ‰ (), as well as the long-term averages for BRWE of 15.2 ± 0.3 ‰ and for AS337382 of 14.0 ± 0.3 ‰. The BRWE.2 standard (used in the later part of the study to replace the BRWE standard) gave a mean of 14.3 ‰ with a between-run reproducibility of 0.3 ‰ (1 SD, n = 17). Samples were usually measured in triplicate and average reproducibility of sample replicate measurements was 0.3 ‰.

Strontium isotope analysis

Sample aliquots reserved for strontium isotope analysis were processed as whole enamel pieces, which were first transferred to a PicoTrace clean laboratory facility at the MPI-EVA and cleaned in the facility before sample preparation. Sample preparation was conducted using a digestion and ion-exchange chromatography protocol following methods outlined in () (see full details in supplementary text S9). All samples were analyzed for 87Sr/86Sr using a Neptune multicollector inductively coupled plasma mass spectrometer (MC-ICPMS, Thermo Fisher Scientific, Bremen, Germany) at the MPI-EVA. Resulting 87Sr/86Sr measurements were normalized for instrumental mass bias to 88Sr/86Sr = 8.375209 (exponential law) and corrected for 87Rb interference. External data normalization was conducted using the NIST SRM 987 reference material [87Sr/86Sr accepted value = 0.710240 (); average of measured values = 0.710283 ± 0.000009, 1 SD, n = 16] using correction offsets of −0.000042 and −0.000043. Measurements of NIST SRM 1486 gave an average value of 0.709296 ± 0.0000056 (1 SD, four measurements of two aliquots), which is very close to the expected value of 0.709299. All samples were measured in duplicate with an average reproducibility of 0.0000068 (1 SD). Two procedural blanks processed alongside samples gave Sr concentrations of ∼0.02% of typical sample concentrations. A lack of correlation between strontium concentration and 87Sr/86Sr values indicates an absence of diagenetic alteration in all samples (fig. S9).

Temperature reconstruction

Air temperature estimates were obtained by a two-step regression process based on the empirically determined relationships between (i) tooth enamel δ18Ophos and drinking water δ18O (δ18Odw) and (ii) δ18O of precipitation (δ18Oprecip) and air temperature. Regression methods and determination of temperature reconstruction uncertainty were applied following (). Details on modern calibration datasets and the conversion procedure are available in supplementary text S8, fig. S10 (map of Global Network of Isotopes in Precipitation (GNIP) stations used for water isotope data), fig. S11 (water isotope to temperature calibration data), and fig. S12 (enamel-drinking water calibration data). An inverse model following () was applied to seasonal sinusoidal curves of δ18Ophos to correct for damping of the seasonal amplitude caused by time averaging from tooth enamel mineralization and the sampling procedure (see supplementary text S7 and fig. S13). The resulting curve minimum and maximum values were converted to air temperature to yield estimates of summer and winter temperature. It should be noted that the inverse model does not take into account the progressive slowing of tooth growth and enamel mineralization that is known to occur particularly in horses. Seasonal amplitudes reconstructed here, therefore, most likely represent minimum amplitudes (see supplementary text S7 for details). Mean annual temperatures were obtained by converting directly from the annual midpoint of unmodeled oxygen isotope means, which were in turn calculated as the average of the summer peak and winter trough value. It should be noted that modern-day observations of mean annual temperatures are calculated as averages of equidistant observations along a full annual cycle rather than as the midpoint between the warmest and the coldest month. Because of the dynamics of tooth enamel mineralization and the sampling geometry, δ18O measurements represent averages over different amounts of time and are not equidistant in time, making it impossible to derive annual averages of δ18O measurements that are conceptually identical to annual temperature means. Analyses discussed in supplementary text S8 show that there is no difference between paleotemperatures obtained from summer/winter midpoint δ18Ophos values compared to averages of all measurements within the complete period of the sinusoidal δ18Ophos curve. We therefore conclude that these methods are equally useful proxies for mean annual temperature and use the summer/winter midpoint to reconstruct such paleotemperature estimates.

Software, code, and data

This article, including code for all data analyses, was written in R version 3.6.2 (), and the manuscript and Supplementary Materials were produced using RMarkdown (). Package information and version details can be found in supplementary text S10. All δ18O and 87Sr/86Sr values for individual samples can be found in tables S2 and S3. Data and code to reproduce the manuscript files, figures, and analyses are available at https://osf.io/tk9dc/.