Subarctic climate for the earliest Homo sapiens in Europe
Abstract
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 (1–3). 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 (1–4). 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 (5, 6). 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 (6) [14C dates recalibrated using IntCal20 (7); 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 (5), suggesting that the Bacho Kiro IUP humans were part of a local population extinction, as previously proposed for other early H. sapiens finds (8, 9). 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 (5). The predominant taxa are cave bear, cervids (especially Cervus elaphus), Bos/Bison, caprines (especially Capra ibex), and equids (5). 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 (5, 10). Preliminary results from the study of micromammals suggest the presence of open habitats with a cooler climate than today (5, 11). 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 (12). 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., (13)].
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 (41)] 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 (14), 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 (13). 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 (15). 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 (16–18) (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 (19–22).
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 (5). 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 (5).
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 [(1–3, 23), but see (24) and (25)]. 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 (2, 26). Chronometric dates of the IUP in Central Asia suggest the appearance of the IUP in the Altai during GI 12 (23, 27), which is supported by pedogenesis indicators in the Tolbor 16 stratigraphic sequence (23). 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 (2, 28). 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 (1). 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 (25), 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 (25). 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 (1, 29), Layers N1-I & I and N1-J & J were originally considered to most likely fall into the warm phase around GI 12 (1, 6, 29). Similarly, a comparison with the Tenaghi Philippon pollen record would also support this assignment (2). 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 (5) 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.
| 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 |
Table 1. Context and sample information for teeth sampled for oxygen and strontium stable isotope analysis.
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 (30) and modified in (31) as described in (32) (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 (33)]. 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 ‰ (34), 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 (35) (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 (36); 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 (37). 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 (38) 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 (39), and the manuscript and Supplementary Materials were produced using RMarkdown (40). 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/.
Acknowledgments
The re-excavation of Bacho Kiro Cave was jointly conducted by the National Institute of Archaeology with Museum, Bulgarian Academy of Sciences, Sofia and the Department of Human Evolution at the MPI-EVA. We would like to thank the National Museum of Natural History (Sofia), the Archaeology Department at the New Bulgarian University (Sofia), the Regional Museum of History in Gabrovo, and the History Museum in Dryanovo for assistance on this project and the opportunity to study the Bacho Kiro Cave faunal material. We would like to thank M. Trost, S. Hesse, M. Kaniecki, and P. Dittmann (MPI-EVA) for technical assistance during stable isotope sample preparation. S. Steinbrenner is thanked for technical assistance with TC/EA-IRMS maintenance. Thanks are also due to H. Temming and U. Schwarz (MPI-EVA) for the production of microCT scans and replicas of the sample materials. We would also like to acknowledge the assistance of to D. Veres with taking OSL samples. Last but not least we would like to thank the handling editor, S. Ortman, as well as three anonymous reviewers for their thoughtful comments that greatly improved this manuscript. Funding: The field work was financed by the Max Planck Society. The stable isotope work was funded by the Max Planck Society as part of S.P.’s doctoral project. S.P. was supported by the Max Planck Society and the University of Aberdeen. K.B. was supported by a Philip Leverhulme Prize from The Leverhulme Trust (PLP-2019-284). N.B.’s work was supported as part of a grant by the German Research Foundation (“PALÄODIET” Project 378496604). V.A. was supported by a grant from the Foundation for Science and Technology, Portugal (IF/01157/2015/CP1308/CT0002). Author contributions: The study was devised by S.P., K.B., S.P.M., J.-J.H., and T.T. Archaeological excavation was undertaken by N.S. and T.T. in collaboration with Z.R. and S.P.M. who all contributed contextual information. V.A. collected sedimentological data at the site and untertook micromorphological investigations that provided information on site formation for this study. Zooarchaeological and paleontological analyses were performed by G.M.S. and R.S. OSL dating was carried out by T.L. Radiocarbon dating and recalibration of radiocarbon dates were conducted by H.F. MC-ICPMS analysis was conducted by N.B. and S.P. Sampling, sample processing for oxygen and strontium stable isotope analysis, and TC/EA-IRMS analysis were carried out by S.P. Code and data analyses were written and conducted by S.P. N.-H.T. consulted on statistical analysis and coding. S.P. wrote the paper with input from all authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. This paper and the Supplementary Materials were generated completely using R and Rmarkdown, making all analyses and figures reproducible via the RMarkdown documents. These documents containing all analysis code, as well as all underlying data files, can be found in an OSF repository at https://osf.io/tk9dc/. All stable isotope data generated here are additionally available in tables S2 and S3 of this article. The faunal remains that were sampled for stable isotope analysis have been returned to the History Museum in Dryanovo and are available to any researcher to be studied.
Supplementary Materials
This PDF file includes:
Supplementary Texts S1 to S10
Figs. S1 to S13
Tables S1 to S3
References
REFERENCES AND NOTES
1
M. Staubwasser, V. Draguşin, B. P. Onac, S. Assonov, V. Ersek, D. L. Hoffmann, D. Veres, Impact of climate change on the transition of Neanderthals to modern humans in Europe. Proc. Natl. Acad. Sci. 115, 9116–9121 (2018).
2
U. C. Müller, J. Pross, P. C. Tzedakis, C. Gamble, U. Kotthoff, G. Schmiedl, S. Wulf, K. Christanis, The role of climate in the spread of modern humans into Europe. Quat. Sci. Rev. 30, 273–279 (2011).
3
J. J. Hublin, The modern human colonization of western Eurasia: When and where? Quat. Sci. Rev. 118, 194–210 (2015).
4
S. L. Kuhn, N. Zwyns, Rethinking the initial Upper Paleolithic. Quat. Int. 347, 29–38 (2014).
5
J.-J. Hublin, N. Sirakov, V. Aldeias, S. Bailey, E. Bard, V. Delvigne, E. Endarova, Y. Fagault, H. Fewlass, M. Hajdinjak, B. Kromer, I. Krumov, J. Marreiros, N. L. Martisius, L. Paskulin, V. Sinet-Mathiot, M. Meyer, S. Pääbo, V. Popov, Z. Rezek, S. Sirakova, M. M. Skinner, G. M. Smith, R. Spasov, S. Talamo, T. Tuna, L. Wacker, F. Welker, A. Wilcke, N. Zahariev, S. P. McPherron, T. Tsanova, Initial upper palaeolithic Homo sapiens from Bacho Kiro Cave, Bulgaria. Nature 581, 299–302 (2020).
6
H. Fewlass, S. Talamo, L. Wacker, B. Kromer, T. Tuna, Y. Fagault, E. Bard, S. P. McPherron, V. Aldeias, R. Maria, N. L. Martisius, L. Paskulin, Z. Rezek, V. Sinet-Mathiot, S. Sirakova, G. M. Smith, R. Spasov, F. Welker, N. Sirakov, T. Tsanova, J. J. Hublin, A 14C chronology for the Middle to Upper Palaeolithic transition at Bacho Kiro Cave, Bulgaria. Nat. Ecol. Evol. 4, 794–801 (2020).
7
P. J. Reimer, W. E. N. Austin, E. Bard, A. Bayliss, P. G. Blackwell, C. B. Ramsey, M. Butzin, H. Cheng, R. L. Edwards, M. Friedrich, P. M. Grootes, T. P. Guilderson, I. Hajdas, T. J. Heaton, A. G. Hogg, K. A. Hughen, B. Kromer, S. W. Manning, R. Muscheler, J. G. Palmer, C. Pearson, J. van der Plicht, R. W. Reimer, D. A. Richards, E. M. Scott, J. R. Southon, C. S. M. Turney, L. Wacker, F. Adolphi, U. Büntgen, M. Capano, S. M. Fahrni, A. Fogtmann-Schulz, R. Friedrich, P. Köhler, S. Kudsk, F. Miyake, J. Olsen, The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0-55 cal kBP). Radiocarbon 62, 725–757 (2020).
8
Q. Fu, M. Hajdinjak, O. T. Moldovan, S. Constantin, S. Mallick, P. Skoglund, N. Patterson, N. Rohland, I. Lazaridis, B. Nickel, B. Viola, K. Prüfer, M. Meyer, J. Kelso, D. Reich, S. Pääbo, An early modern human from Romania with a recent Neanderthal ancestor. Nature 524, 216–219 (2015).
9
Q. Fu, C. Posth, M. Hajdinjak, M. Petr, S. Mallick, D. Fernandes, A. Furtwängler, W. Haak, M. Meyer, A. Mittnik, B. Nickel, A. Peltzer, N. Rohland, V. Slon, S. Talamo, I. Lazaridis, M. Lipson, I. Mathieson, S. Schiffels, P. Skoglund, A. P. Derevianko, N. Drozdov, V. Slavinsky, A. Tsybankov, R. G. Cremonesi, F. Mallegni, B. Gély, E. Vacca, M. R. G. Morales, L. G. Straus, C. Neugebauer-Maresch, M. Teschler-Nicola, S. Constantin, O. T. Moldovan, S. Benazzi, M. Peresani, D. Coppola, M. Lari, S. Ricci, A. Ronchitelli, F. Valentin, C. Thevenet, K. Wehrberger, D. Grigorescu, H. Rougier, I. Crevecoeur, D. Flas, P. Semal, M. A. Mannino, C. Cupillard, H. Bocherens, N. J. Conard, K. Harvati, V. Moiseyev, D. G. Drucker, J. Svoboda, M. P. Richards, D. Caramelli, R. Pinhasi, J. Kelso, N. Patterson, J. Krause, S. Pääbo, D. Reich, The genetic history of Ice Age Europe. Nature 534, 200–205 (2016).
10
J. K. Kozłowski, B. Ginter, Excavation in the Bacho Kiro Cave (Bulgaria). Final report (Państwowe Wydawnictwo Naukowe, 1982).
11
V. V. Popov, in The Pleistocene: Geography, GEOLOGY, and Fauna, G. Huard, J. Gareau, Eds. (2018), pp. 109–236.
12
Z. Chen, A. S. Auler, M. Bakalowicz, D. Drew, F. Griger, J. Hartmann, G. Jiang, N. Moosdorf, A. Richts, Z. Stevanovic, G. Veni, N. Goldscheider, The World Karst Aquifer Mapping project: Concept, mapping procedure and map of Europe. Hydrogeol. J. 25, 771–785 (2017).
13
L. Scherler, T. Tütken, D. Becker, Carbon and oxygen stable isotope compositions of late Pleistocene mammal teeth from dolines of Ajoie (Northwestern Switzerland). Quat. Res. 82, 378–387 (2014).
14
D. Pushkina, H. Bocherens, R. Ziegler, Unexpected palaeoecological features of the Middle and Late Pleistocene large herbivores in southwestern Germany revealed by stable isotopic abundances in tooth enamel. Quat. Int. 339-340, 164–178 (2014).
15
L. M. Arppe, J. A. Karhu, Oxygen isotope values of precipitation and the thermal climate in Europe during the middle to late Weichselian ice age. Quat. Sci. Rev. 29, 1263–1275 (2010).
16
G. J. Bowen, J. Revenaugh, Interpolating the isotopic composition of modern meteoric precipitation. Water Resour. Res. 39, 1299 (2003).
17
G. J. Bowen, L. I. Wassenaar, K. A. Hobson, Global application of stable hydrogen and oxygen isotopes to wildlife forensics. Oecologia 143, 337–348 (2005).
18
IAEA/WMO, Global Network of Isotopes in Precipitation. The GNIP Database (2020); https://nucleus.iaea.org/wiser.
19
I. Obreht, C. Zeeden, U. Hambach, D. Veres, S. B. Marković, J. Bösken, Z. Svirčev, N. Bačević, M. B. Gavrilov, F. Lehmkuhl, Tracing the influence of Mediterranean climate on Southeastern Europe during the past 350,000 years. Sci. Rep. 6, 36334 (2016).
20
A. Wegwerth, A. Ganopolski, G. Ménot, J. Kaiser, O. Dellwig, E. Bard, F. Lamy, H. W. Arz, Black Sea temperature response to glacial millennial-scale climate variability. Geophys. Res. Lett. 42, 8147–8154 (2015).
21
A. Feurdean, A. Perşoiu, I. Tanţău, T. Stevens, E. K. Magyari, B. P. Onac, S. Marković, M. Andrič, S. Connor, S. Fărcaş, M. Gałka, T. Gaudeny, W. Hoek, P. Kolaczek, P. Kuneš, M. Lamentowicz, E. Marinova, D. J. Michczyńska, I. Perşoiu, M. Płóciennik, M. Słowiński, M. Stancikaite, P. Sumegi, A. Svensson, T. Tămaş, A. Timar, S. Tonkov, M. Toth, S. Veski, K. J. Willis, V. Zernitskaya, Climate variability and associated vegetation response throughout Central and Eastern Europe (CEE) between 60 and 8ka. Quat. Sci. Rev. 106, 206–224 (2014).
22
E. Kjellström, J. Brandefelt, J. O. Näslund, B. Smith, G. Strandberg, A. H. L. Voelker, B. Wohlfarth, Simulated climate conditions in Europe during the Marine Isotope Stage 3 stadial. Boreas 39, 436–456 (2010).
23
N. Zwyns, C. H. Paine, B. Tsedendorj, S. Talamo, K. E. Fitzsimmons, A. Gantumur, L. Guunii, O. Davakhuu, D. Flas, T. Dogandžić, N. Doerschner, F. Welker, J. C. Gillam, J. B. Noyer, R. S. Bakhtiary, A. F. Allshouse, K. N. Smith, A. M. Khatsenovich, E. P. Rybin, G. Byambaa, J. J. Hublin, The Northern route for human dispersal in central and northeast asia: New evidence from the site of Tolbor-16, Mongolia. Sci. Rep. 9, 11759 (2019).
24
D. Mihailović, Push-and-pull factors of the Middle to Upper Paleolithic transition in the Balkans. Quat. Int. 551, 47–62 (2020).
25
P. R. Nigst, P. Haesaerts, F. Damblon, C. Frank-Fellner, C. Mallol, B. Viola, M. Götzinger, L. Niven, G. Trnka, J.-J. Hublin, Early modern human settlement of Europe north of the Alps occurred 43,500 years ago in a cold steppe-type environment. Proc. Natl. Acad. Sci. U.S.A. 111, 14394–14399 (2014).
26
D. Richter, G. Tostevin, P. Škrdla, W. Davies, New radiometric ages for the Early Upper Palaeolithic type locality of Brno-Bohunice (Czech Republic): Comparison of OSL, IRSL, TL and 14C dating results. J. Archaeol. Sci. 36, 708–720 (2009).
27
N. Zwyns, thesis, Leiden University (2012).
28
W. Davies, P. Gollop, The human presence in Europe during the Last Glacial Period II. Climate tolerance and climate preferences of Mid-and Late Glacial hominids, in Neanderthals and Modern Humans in the European Landscape of the Last Glaciation (The McDonald Institute for Archaeological Research, 2003), pp. 131–146.
29
S. Constantin, A. V. Bojar, S. E. Lauritzen, J. Lundberg, Holocene and Late Pleistocene climate in the sub-Mediterranean continental environment: A speleothem record from Poleva Cave (Southern Carpathians, Romania). Palaeogeogr. Palaeoclimatol. Palaeoecol. 243, 322–338 (2007).
30
D. L. Dettman, M. J. Kohn, J. Quade, F. J. Ryerson, T. P. Ojha, S. Hamidullah, Seasonal stable isotope evidence for a strong Asian monsoon throughout the past 10.7 m.y. Geology 29, 31–34 (2001).
31
T. Tütken, T. W. Vennemann, H. Janz, E. P. J. Heizmann, Palaeoenvironment and palaeoclimate of the Middle Miocene lake in the Steinheim basin, SW Germany: A reconstruction from C, O, and Sr isotopes of fossil remains. Palaeogeogr. Palaeoclimatol. Palaeoecol. 241, 457–491 (2006).
32
K. Britton, B. T. Fuller, T. Tütken, S. Mays, M. P. Richards, Oxygen isotope analysis of human bone phosphate evidences weaning age in archaeological populations. Am. J. Phys. Anthropol. 157, 226–241 (2015).
33
W. A. Brand, T. B. Coplen, A. T. Aerts-Bijma, J. K. Böhlke, M. Gehre, H. Geilmann, M. Gröning, H. G. Jansen, H. A. J. Meijer, S. J. Mroczkowski, H. Qi, K. Soergel, H. Stuart-Williams, S. M. Weise, R. A. Werner, Comprehensive inter-laboratory calibration of reference materials for δ18O versus VSMOW using various on-line high-temperature conversion techniques. Rapid Commun. Mass Spectrom. 23, 999–1019 (2009).
34
E. Pucéat, M. M. Joachimski, A. Bouilloux, F. Monna, A. Bonin, S. Motreuil, P. Morinière, S. Hénard, J. Mourin, G. Dera, D. Quesne, Revised phosphate-water fractionation equation reassessing paleotemperatures derived from biogenic apatite. Earth Planet. Sci. Lett. 298, 135–142 (2010).
35
S. R. Copeland, M. Sponheimer, P. J. le Roux, V. Grimes, J. A. Lee-Thorp, D. J. de Ruiter, M. P. Richards, Strontium isotope ratios (87Sr/86Sr) of tooth enamel: A comparison of solution and laser ablation multicollector inductively coupled plasma mass spectrometry methods. Rapid Commun. Mass Spectrom. 22, 3187–3194 (2008).
36
C. M. Johnson, C. J. Fridrich, Non-monotonic chemical and O, Sr, Nd, and Pb isotope zonations and heterogeneity in the mafic- to silicic-composition magma chamber of the Grizzly Peak Tuff, Colorado. Contrib. Mineral. Petrol. 105, 677–690 (1990).
37
A. J. E. Pryor, R. E. Stevens, T. C. O’Connell, J. R. Lister, Quantification and propagation of errors when converting vertebrate biomineral oxygen isotope data to temperature for palaeoclimate reconstruction. Palaeogeogr. Palaeoclimatol. Palaeoecol. 412, 99–107 (2014).
38
B. H. Passey, T. E. Cerling, G. T. Schuster, T. F. Robinson, B. L. Roeder, S. K. Krueger, Inverse methods for estimating primary input signals from time-averaged isotope profiles. Geochim. Cosmochim. Acta 69, 4101–4116 (2005).
39
R Core Team, R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, Austria, 2020; https://r-project.org/).
40
J. Allaire, J. Cheng, Y. Xie, J. McPherson, W. Chang, J. Allen, R. Hyndman, rmarkdown: Dynamic documents for R (R package version 1.0) (2018).
41
N. Center Aviation Weather, Gorna Oryahovitsa airport (LBGO) METAR (2012–2020) (2020); https://aviationweather.gov/metar.
42
K. A. Hoppe, S. M. Stover, J. R. Pascoe, R. Amundson, Tooth enamel biomineralization in extant horses: Implications for isotopic microsampling. Palaeogeogr. Palaeoclimatol. Palaeoecol. 206, 355–365 (2004).
43
C. Gadbury, L. Todd, A. H. Jahren, R. Amundson, Spatial and temporal variations in the isotopic composition of bison tooth enamel from the Early Holocene Hudson-Meng Bone Bed, Nebraska. Palaeogeogr. Palaeoclimatol. Palaeoecol. 157, 79–93 (2000).
44
D. A. E. Garrod, B. Howe, J. H. Gaul, Excavations in the cave of Bacho Kiro, north-east Bulgaria. Bulletin of the American School of Prehistoric Research. 15, 46–126 (1939).
45
C. Bronk Ramsey, Bayesian analysis of radiocarbon dates. Radiocarbon 51, 337–360 (2009).
46
M. J. Aitken, Introduction to Optical Dating: The Dating of Quaternary Sediments by the Use of Photon-stimulated Luminescence (Clarendon Press, 1998).
47
J. Buylaert, M. Jain, A. S. Murray, K. J. Thomsen, C. Thiel, R. Sohbati, A robust feldspar luminescence dating method for Middle and Late Pleistocene sediments. Boreas 41, 435–451 (2012).
48
J.-P. Buylaert, A. S. Murray, K. J. Thomsen, M. Jain, Testing the potential of an elevated temperature IRSL signal from K-feldspar. Radiat. Meas. 44, 560–565 (2009).
49
C. Thiel, J.-P. Buylaert, A. Murray, B. Terhorst, I. Hofer, S. Tsukamoto, M. Frechen, Luminescence dating of the Stratzing loess profile (Austria)—Testing the potential of an elevated temperature post-IR IRSL protocol. Quat. Int. 234, 23–31 (2011).
50
K. J. Thomsen, A. S. Murray, M. Jain, L. Bøtter-Jensen, Laboratory fading rates of various luminescence signals from feldspar-rich sediment extracts. Radiat. Meas. 43, 1474–1486 (2008).
51
D. J. Huntley, M. Lamothe, Ubiquity of anomalous fading in K-feldspars and the measurement and correction for it in optical dating. Can. J. Earth Sci. 38, 1093–1106 (2001).
52
S. Kreutzer, C. Schmidt, R. DeWitt, M. Fuchs, The a-value of polymineral fine grain samples measured with the post-IR IRSL protocol. Radiat. Meas. 69, 18–29 (2014).
53
I. D. Clark, P. Fritz, Environmental Isotopes in Hydrogeology (Lewis Publishers, 1997).
54
K. Rozanski, L. Araguás-Araguás, R. Gonfiantini, Isotopic patterns in modern global precipitation. Clim. Chang. Cont. 78, 1–36 (1993).
55
S. Müller, C. Stumpp, J. H. Sørensen, S. Jessen, Spatiotemporal variation of stable isotopic composition in precipitation: Post-condensational effects in a humid area. Hydrol. Process. 31, 3146–3159 (2017).
56
C. Stumpp, J. Klaus, W. Stichler, Analysis of long-term stable isotopic composition in German precipitation. J. Hydrol. 517, 351–361 (2014).
57
W. Dansgaard, Stable isotopes in precipitation. Tellus 16, 436–468 (1964).
58
M. J. Kohn, J. M. Welker, On the temperature correlation of δ18O in modern precipitation. Earth Planet. Sci. Lett. 231, 87–96 (2005).
59
M. J. Kohn, M. J. Schoeninger, J. W. Valley, Herbivore tooth oxygen isotope compositions: Effects of diet and physiology. Geochim. Cosmochim. Acta 60, 3889–3896 (1996).
60
B. Luz, Y. Kolodny, M. Horowitz, Fractionation of oxygen isotopes between mammalian bone-phosphate and environmental drinking water. Geochim. Cosmochim. Acta 48, 1689–1693 (1984).
61
A. Longinelli, Oxygen isotopes in mammal bone phosphate: A new tool for paleohydrological and paleoclimatological research? Geochim. Cosmochim. Acta 48, 385–390 (1984).
62
P. Iacumin, H. Bocherens, Oxygen isotope analyses of co-existing carbonate and phosphate in biogenic apatite: A way to monitor diagenetic alteration of bone phosphate? Earth Planet. Sci. Lett. 142, 1–6 (1996).
63
T. Tütken, H. Furrer, T. Walter Vennemann, Stable isotope compositions of mammoth teeth from Niederweningen, Switzerland: Implications for the Late Pleistocene climate, environment, and diet. Quat. Int. 164-165, 139–150 (2007).
64
G. Skrzypek, R. Sadler, A. Wiśniewski, Reassessment of recommendations for processing mammal phosphate δ18O data for paleotemperature reconstruction. Palaeogeogr. Palaeoclimatol. Palaeoecol. 446, 162–167 (2016).
65
K. A. Hoppe, Correlation between the oxygen isotope ratio of North American bison teeth and local waters: Implication for paleoclimatic reconstructions. Earth Planet. Sci. Lett. 244, 408–417 (2006).
66
K. Rozanski, L. Araguás-Araguás, R. Gonfiantini, Relation between long-term trends of oxygen-18 isotope composition of precipitation and climate. Science 258, 981–985 (1992).
67
S. Pederzani, K. Britton, Oxygen isotopes in bioarchaeology: Principles and applications, challenges and opportunities. Earth Sci. Rev. 188, 77–107 (2019).
68
G. M. O. Maloiy, Water metabolism of East African ruminants in arid and semi-arid regions. Z. Tierzuecht. Zuechtungsbiol. 90, 219–228 (1973).
69
R. A. Arias, T. L. Mader, Environmental factors affecting daily water intake on cattle finished in feedlots. J. Anim. Sci. 89, 245–251 (2011).
70
T. McHugh, Social behavior of the American buffalo (Bison bison bison). Fortschr. Zool. 43, 1–40 (1958).
71
C. F. Winchester, M. J. Morris, Water intake rates of cattle. J. Anim. Sci. 15, 722–740 (1956).
72
S. Groenendyk, P. B. English, I. Abetz, External balance of water and electrolytes in the horse. Equine Vet. J. 20, 189–193 (1988).
73
K. M. Scheibe, K. Eichhorn, B. Kalz, W. J. Streich, A. Scheibe, Water consumption and watering behavior of Przewalski horses (Equus ferus przewalskii) in a semireserve. Zoo Biol. 17, 181–192 (1998).
74
K. A. Houpt, A. Eggleston, K. Kunkle, T. R. Houpt, Effect of water restriction on equine behaviour and physiology. Equine Vet. J. 32, 341–344 (2000).
75
M. Hinton, On the watering of horses: A review. Equine Vet. J. 10, 27–31 (1978).
76
D. D’Angela, A. Longinelli, Oxygen isotopes in living mammal’s bone phosphate: Further results. Chem. Geol. 86, 75–82 (1990).
77
J. D. Bryant, B. Luz, P. N. Froelich, Oxygen isotopic composition of fossil horse tooth phosphate as a record of continental paleoclimate. Palaeogeogr. Palaeoclimatol. Palaeoecol. 107, 303–316 (1994).
78
J. D. Bryant, P. N. Froelich, W. J. Showers, B. J. Genna, A tale of two quarries: Biologic and taphonomic signatures in the oxygen isotope composition of tooth enamel phosphate from modern and Miocene equids. Palaios 11, 397–408 (1996).
79
K. A. Hoppe, R. Amundson, M. Vavra, M. P. McClaran, D. L. Anderson, Isotopic analysis of tooth enamel carbonate from modern North American feral horses: Implications for paleoenvironmental reconstructions. Palaeogeogr. Palaeoclimatol. Palaeoecol. 203, 299–311 (2004).
80
B. Sánchez Chillón, M. T. Alberdi, G. Leone, F. P. Bonadonna, B. Stenni, A. Longinelli, Oxygen isotopic composition of fossil equid tooth and bone phosphate: An archive of difficult interpretation. Palaeogeogr. Palaeoclimatol. Palaeoecol. 107, 317–328 (1994).
81
J. R. Gat, Isotope Hydrology: A Study of the Water Cycle (Imperila College Press, 2010), vol. 6 of Series on Environmental Science and Management.
82
W. G. Darling, A. H. Bath, J. C. Talbot, The O & H stable isotopic composition of fresh waters in the British Isles. 2. Surface waters and groundwater. Hydrol. Earth Syst. Sci. 7, 183–195 (2003).
83
J. Halder, S. Terzer, L. I. Wassenaar, L. Araguás-Araguás, P. K. Aggarwal, The Global Network of Isotopes in Rivers (GNIR): Integration of water isotopes in watershed observation and riverine research. Hydrol. Earth Syst. Sci. 19, 3419–3431 (2015).
84
R. Gonfiantini, Environmental isotopes in lake studies, in Handbook of Environmental Isotope Geochemistry; The Terrestrial Environment (1986), pp. 113–168.
85
Bulgarian Academy of Sciences, Dryanovska river, in Golyama entsiklopediya—Bulgaria (2012), pp. 1870–1871.
86
E. Bojilova, Inter-annual distribution for Yantra river basin, North Bulgaria. 18 Conference of the Danubian Countries on Hydrological Forcasting and Hydrological Bases of Water Management (2017).
87
G. J. Bowen, OIPC: The online isotopes in precipitation calculator, version 3.1 (2020); http://waterisotopes.org.
88
N. M. Kortelainen, J. A. Karhu, Regional and seasonal trends in the oxygen and hydrogen isotope ratios of Finnish groundwaters: A key for mean annual precipitation. J. Hydrol. 285, 143–157 (2004).
89
K. Rozanski, Deuterium and oxygen-18 in European groundwaters—Links to atmospheric circulation in the past. Chem. Geol. 52, 349–363 (1985).
90
R. van Geldern, A. Baier, H. L. Subert, S. Kowol, L. Balk, J. A. C. Barth, Pleistocene paleo-groundwater as a pristine fresh water resource in southern Germany—Evidence from stable and radiogenic isotopes. Sci. Total Environ. 496, 107–115 (2014).
91
P. K. Aggarwal, O. A. Alduchov, K. O. Froehlich, L. Araguás-Araguás, N. C. Sturchio, N. Kurita, Stable isotopes in global precipitation: A unified interpretation based on atmospheric moisture residence time. Geophys. Res. Lett. 39, L11705 (2012).
92
G. J. Bowen, J. B. West, in Tracking Animal Migration with Stable Isotopes, K. A. Hobson, L. I. Wassenaar, Eds. (Elsevier, 2008), vol. 2 of Terrestrial Ecology, pp. 79–105.
93
P. D. Akers, J. M. Welker, G. A. Brook, Reassessing the role of temperature in precipitation oxygen isotopes across the eastern and central United States through weekly precipitation-day data. Water Resour. Res. 53, 7644–7661 (2017).
94
G. A. Schmidt, G. Hoffmann, D. T. Shindell, Y. Hu, Modeling atmospheric stable water isotopes and the potential for constraining cloud processes and stratosphere-troposphere water exchange. J. Geophys. Res.-Atmos. 110, D21314 (2005).
95
A. Zuber, S. M. Weise, J. Motyka, K. Osenbrück, K. Rózański, Age and flow pattern of groundwater in a Jurassic limestone aquifer and related Tertiary sands derived from combined isotope, noble gas and chemical data. J. Hydrol. 286, 87–112 (2004).
96
F. Kaspar, N. Kühl, U. Cubasch, T. Litt, A model-data comparison of European temperatures in the Eemian interglacial. Geophys. Res. Lett. 32, L11703 (2005).
97
G. Skrzypek, A. Winiewski, P. F. Grierson, How cold was it for Neanderthals moving to Central Europe during warm phases of the last glaciation? Quat. Sci. Rev. 30, 481–487 (2011).
98
K. Britton, V. Grimes, L. Niven, T. E. Steele, S. J. P. McPherron, M. Soressi, T. E. Kelly, J. Jaubert, J.-J. Hublin, M. P. Richards, Strontium isotope evidence for migration in late Pleistocene Rangifer: Implications for Neanderthal hunting strategies at the Middle Palaeolithic site of Jonzac, France. J. Hum. Evol. 61, 176–185 (2011).
99
M. C. Dean, Growth layers and incremental markings in hard tissues; a review of the literature and some preliminary observations about enamel structure in Paranthropus boisei. J. Hum. Evol. 16, 157–172 (1987).
100
S. Hillson, Teeth (Cambridge, 2012).
101
H. C. Fricke, J. R. O’Neil, Inter- and intra-tooth variation in the oxygen isotope composition of mammalian tooth enamel phosphate: Implications for palaeoclimatological and palaeobiological research. Palaeogeogr. Palaeoclimatol. Palaeoecol. 126, 91–99 (1996).
102
H. C. Fricke, W. C. Clyde, J. R. O’Neil, Intra-tooth variations in δ18O(PO4) of mammalian tooth enamel as a record of seasonal variations in continental climate variables. Geochim. Cosmochim. Acta 62, 1839–1850 (1998).
103
S. A. Blumenthal, T. E. Cerling, K. L. Chritz, T. G. Bromage, R. Kozdon, J. W. Valley, Stable isotope time-series in mammalian teeth: In situ δ18O from the innermost enamel layer. Geochim. Cosmochim. Acta 124, 223–236 (2014).
104
M. Balasse, Reconstructing dietary and environmental history from enamel isotopic analysis: Time resolution of intra-tooth sequential sampling. Int. J. Osteoarchaeol. 12, 155–165 (2002).
105
D. R. Green, G. M. Green, A. S. Colman, F. B. Bidlack, P. Tafforeau, T. M. Smith, Synchrotron imaging and Markov Chain Monte Carlo reveal tooth mineralization patterns. PLOS ONE 12, e0186391 (2017).
106
M. J. Kohn, Comment: Tooth enamel mineralization in ungulates: Implications for recovering a primary isotopic time-series, by B. H. Passey and T. E. Cerling (2002). Geochim. Cosmochim. Acta 68, 403–405 (2004).
107
B. H. Passey, T. E. Cerling, Tooth enamel mineralization in ungulates: Implications for recovering a primary isotopic time-series. Geochim. Cosmochim. Acta 66, 3225–3234 (2002).
108
R. B. Trayler, M. J. Kohn, Tooth enamel maturation reequilibrates oxygen isotope compositions and supports simple sampling methods. Geochim. Cosmochim. Acta 198, 32–47 (2016).
109
R. Bendrey, D. Vella, A. Zazzo, M. Balasse, S. Lepetz, Exponentially decreasing tooth growth rate in horse teeth: Implications for isotopic analyses. Archaeometry 57, 1104–1124 (2015).
110
D. R. Green, T. M. Smith, G. M. Green, F. B. Bidlack, P. Tafforeau, A. S. Colman, Quantitative reconstruction of seasonality from stable isotopes in teeth. Geochim. Cosmochim. Acta 235, 483–504 (2018).
111
S. A. Blumenthal, T. E. Cerling, T. M. Smiley, C. E. Badgley, T. W. Plummer, Isotopic records of climate seasonality in equid teeth. Geochim. Cosmochim. Acta 260, 329–348 (2019).
112
A. Zazzo, R. Bendrey, D. Vella, A. P. Moloney, F. J. Monahan, O. Schmidt, A refined sampling strategy for intra-tooth stable isotope analysis of mammalian enamel. Geochim. Cosmochim. Acta 84, 1–13 (2012).
113
A. Delgado Huertas, P. Iacumin, B. Stenni, B. Sánchez Chillón, A. Longinelli, Oxygen isotope variations of phosphate in mammalian bone and tooth enamel. Geochim. Cosmochim. Acta 59, 4299–4305 (1995).
114
R. D. Field, Observed and modeled controls on precipitation δ18O over Europe: From local temperature to the Northern Annular Mode. J. Geophys. Res. 115, D12101 (2010).
115
J. M. McArthur, R. J. Howarth, T. R. Bailey, Strontium isotope stratigraphy: LOWESS version 3: Best fit to the marine Sr-isotope curve for 0-509 Ma and accompanying look-up table for deriving numerical age. J. Geol. 109, 155–170 (2001).
116
A. Kassambara, rstatix: Pipe-Friendly Framework for Basic Statistical Tests (2020); https://cran.r-project.org/package=rstatix.
117
H. Wickham, L. Henry, tidyr: Tidy Messy Data (2019); https://cran.r-project.org/package=tidyr.
118
L. Alathea, captioner: Numbers Figures and Creates Simple Captions (2015); https://cran.r-project.org/package=captioner.
119
D. Gohel, flextable: Functions for Tabular Reporting (2020); https://cran.r-project.org/package=flextable.
120
H. Wickham, ggplot2: Elegant Graphics for Data Analysis (Springer New York, 2016).
121
H. Wickham, stringr: Simple, Consistent Wrappers for Common String Operations (2019); https://cran.r-project.org/package=stringr.
122
H. Wickham, R. François, L. Henry, K. Müller, dplyr: A Grammar of Data Manipulation (2019); https://cran.r-project.org/package=dplyr.
123
Y. Xie, knitr: A comprehensive tool for reproducible research in R. Implement Reprod Res. 1, 20 (2014).
124
D. Gohel, officer: Manipulation of Microsoft Word and PowerPoint Documents (2020); https://cran.r-project.org/package=officer.
125
D. Kahle, H. Wickham, ggmap: Spatial Visualization with ggplot2. R J. 5, 144–161 (2013).
126
G. Daróczi, R. Tsegelskyi, pander: An R ‘Pandoc’ Writer (2018); https://cran.r-project.org/package=pander.
127
B. Auguie, gridExtra: Miscellaneous Functions for "Grid" Graphics (2017); https://cran.r-project.org/package=gridExtra.
128
J. Ooms, magick: Advanced Graphics and Image-Processing in R (2020); https://cran.r-project.org/package=magick.
129
T. L. Pedersen, patchwork: The Composer of Plots (2020); https://cran.r-project.org/package=patchwork.
130
A. Kassambara, ggpubr: ‘ggplot2’ Based Publication Ready Plots (2020); https://cran.r-project.org/package=ggpubr.
131
C. O. Wilke, cowplot: Streamlined Plot Theme and Plot Annotations for ‘ggplot2’ (2019); https://cran.r-project.org/package=cowplot.
132
N. Ross, redoc: Reversible Reproducible Documents (2019); https://github.com/noamross/redoc.
133
OpenStreetMap Contributors, Planet dump retrieved from https://planet.osm.org (2020).
Information & Authors
Information
Published In
Science Advances
Volume 7 • Issue 39 • September 2021
History
Received: 19 March 2021
Accepted: 3 August 2021
Published online: 22 September 2021
Copyright
Copyright © 2021 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).
This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.
Authors
Funding Information
http://dx.doi.org/10.13039/501100001659Deutsche Forschungsgemeinschaft: 378496604
http://dx.doi.org/10.13039/501100001871Fundação para a Ciência e a Tecnologia: IF/01157/2015/CP1308/CT0002
http://dx.doi.org/10.13039/501100004189Max-Planck-Gesellschaft: