Reconstructing Late Pleistocene paleoclimate at the scale of human behavior: an example from the Neandertal occupation of La Ferrassie (France)

Exploring the role of changing climates in human evolution is currently impeded by a scarcity of climatic information at the same temporal scale as the human behaviors documented in archaeological sites. This is mainly caused by high uncertainties in the chronometric dates used to correlate long-term climatic records with archaeological deposits. One solution is to generate climatic data directly from archaeological materials representing human behavior. Here we use oxygen isotope measurements of Bos/Bison tooth enamel to reconstruct summer and winter temperatures in the Late Pleistocene when Neandertals were using the site of La Ferrassie. Our results indicate that, despite the generally cold conditions of the broader period and despite direct evidence for cold features in certain sediments at the site, Neandertals used the site predominantly when climatic conditions were mild, similar to conditions in modern day France. We suggest that due to millennial scale climate variability, the periods of human activity and their climatic characteristics may not be representative of average conditions inferred from chronological correlations with long-term climatic records. These results highlight the importance of using direct routes, such as the high-resolution archives in tooth enamel from anthropogenically accumulated faunal assemblages, to establish climatic conditions at a human scale.


Site background
The site of La Ferrassie is situated in the Dordogne region of southwest France, close to the commune of Le Bugue. La Ferrassie comprises a complex of several sites: the 'Grand Abri' (large rock shelter), a smaller rock shelter and a small cave, (Upper Cave). In this study we focus on the Grand Abri (hereafter La Ferrassie) which was excavated by Capitan and Peyrony and later by Delporte. [1][2][3][4] Our study is based on samples from the recent excavations by Turq and colleagues. 5  However, in the unmodeled ages there was an age inversion in the radiocarbon dates of Layers 4 and 5 with no obvious resolution. 8 While the quality indicators of the 14 C dates and extracted collagen do not suggest any issues with the reliability of the dates, the geological features of the sediments do not indicate a reworking of Layers 4 and 5. 8 At the same time, the 14 C ages for each layer are internally very consistent, indicating a lack of mixing with other deposits and making a reworking hypothesis less likely. Unfortunately, this contradiction cannot be solved at this time, meaning that the ages of Layer 4 should be treated with caution at this point. Based on the geological analyses of the layer sediments we retain the stratigraphic order of these layers when presenting stable isotope results. It should also be noted that for the purposes of the interpretations in our paper, smaller shifts of the age of  this discrepancy is likely due to selection bias present in these older collections. 9 Supplementary Figure S1:Excavations of the Western Sector uncovered a stratigraphic section (A) of 9 Middle and Upper Paleolithic layers (Layer 9 not shown due to spatially restricted extent of this layer) with ages covering late MIS 5 to early MIS 2 (A). Solid black lines in the stratigraphy schematic denote layer boundaries, while dashed white lines demarcate sub-unit divisions. The Middle Paleolithic layers examined in this study are dominated by Levallois blank production with moderate to high proportions of scrapers (B; the orange portion indicates percent scrapers out of the total number of tools per layer). Radiocarbon dates are presented following Model 2 in 8. Note that there was an age inversion in the unmodelled 14 C dates between Layer 4 and 5 that remains unresolved. Radiocarbon ages for Layer 4 therefore should be treated with some caution.
The density of anthropogenic artifacts as exemplified by piece plotted (larger than 2.5 cm) bone fragments and stone tools shows a noticeable change across the Middle Paleolithic sequence, with a pronounced increase in artifact density from Layer 4 onward, with the highest artifact density in Layer 5B (Supplementary Figure S2). Indeed, an increase in the amount of anthropogenic inputs was used to differentiate between Layer 5A and 5B, layers that in terms of their sedimentary matrix are otherwise considered to be the same stratigraphic unit.
Supplementary Figure S2:The density (count per liter of sediment) of piece plotted bone fragments and lithics steadily increases throughout the Middle Paleolithic sequence of La Ferrassie, with highest densities in Layer 5B.

Geology
During the latest excavations, we recognized nine stratigraphic units based on standard lithological criteria, such as color, composition, texture, and the geometry of the sedimentary body. We provide in Supplementary Table S1 descriptions and working interpretations of   Figure S3 and Supplementary Figure S4). Stratigraphic contact and orientations analyses, alongside the microscopic evidence of silt cappings and banded fabric and field observations, clearly associates the deposition of Layer 2 with pronounced cold conditions by solifluction -a downslope slow creep process typical of periglacial landscapes. It can be excluded that cryoturbation features in this Layer represent a reworking by later freeze-thaw processes. The last 'stratigraphic phase' is the formation of ice lenses and silt cappings, which show no signs of disturbance since their formation.
Supplementary Figure S3:a) Field view of Layers 1 (reddish brown) and 2 (yellowish white) in Squares H5-G5 showing cryoturbated aspect of Layer 2 and its sharp, undulating contact with underlying Layer 1 indicative of solifluction. Note the finely comminuted nature of the limestone particles in Layer 2. Scale has 10 cm major increments, 1 cm minor increments; b) Photomicrograph of Layer 2 deposits showing ice lensing (red arrows) that produces banded fabric (lenticular aggregates separated by segregation voids) in silty granular deposits rich in limestone sand; a bone (B) is at upper right. Plane polarized light (PPL), scale is 1 mm; c) Photomicrograph of Layer 2 deposits in cross-polarized light (XPL). Note the abundance of angular calcareous sand (red arrows) indicative of mechanical disintegration of limestone produced by freeze-thaw. Scale is 1 mm.
Supplementary Figure S4:a) Thin section scan (dark field illumination) of sample LF10-1 from Layer 2 illustrating silt cappings on fresh limestone clasts (white arrows). Size of thin section is 75 x 50 mm. b) Photomicrograph of Layer 2 deposits showing thick calcareous silt capping, indicatirng cold, freeze-thaw conditions, on an oyster fragment derived from the local limestone (yellow arrow). Plane-polarized light (PPL); scale is 1 mm.
Layer 5 is subdivided into two subunits, 5A and 5B, that are similar in their sedimentary matrix but differ in terms of fine fraction and density of anthropogenic inputs, the latter are more abundant in 5A. The deposits vary in thickness from 50 cm to ~30 cm near the wall and are comprised of silty sands with dm-sized roof fall. Deposition relates to slope dry-fall processes with a depositional cone emanating from the platform existing somewhere in front of the Upper Cave. Both field and micromorphological observations do not reveal any indications of cold-associated features.

Faunal spectrum
The relative proportion of the main prey taxa in the faunal assemblage (

Bone surface modifications
Surface modifications on bone fragments (> 2.5 cm size cut-off) are dominated by human modifications (cut-marks, scraping, percussion notches, anvil-marks, use as a retoucher) This suggests a relative stability in the intensity of marrow extraction and carcass processing intensity to extract nutrients. This may also indicate that nutritional stress of Neandertals at the site varied little over time.
Supplementary Figure S6:Relative frequencies of human (blue) and carnivore (red) surface modifications on ungulate bones relative to total NISP (printed below bars) grouped by layer (panels) underline the anthropogenic nature of the La Ferrassie faunal assemblage. Bar heights indicate the proportion of fragments with a particular modification and total NISP counts exclude isolated teeth and antler fragments. Small ungulates include roe deer, medium ungulates include red deer and reindeer, large ungulates include Bos/Bison and horses. Prior to sampling, teeth dental calculus was removed in the sampling area using gentle abrasion with a diamond tipped drill bit. Teeth were then repeatedly sonicated in Milli-Q ultrapure water to remove sediment residue. To obtain tooth enamel samples for oxygen and strontium isotope analysis sequential samples were drilled in small strips (ca. 8 x 1.5 x 0.7 mm) perpendicular to the tooth growth axis using a diamond tipped drill bit. Series of sequential samples were drilled covering the complete crown length of one tooth loph. An example of a tooth before and after sampling can be seen in Supplementary Figure  solution. The resulting precipitate was pelleted by centrifugation (12000 rpm for 5 Min) and the silver nitrate supernatant removed. The silver phosphate was then washed three times with MilliQ ultrapure water using centrifugation and vortex mixing steps between rinses to eliminate any silver nitrate from the sample. Silver phosphate samples were then dried over night at 50 °C and stored over desiccant until further analysis.

Oxygen stable isotope methodology
Oxygen isotope ratios of Ag3PO4 were analyzed using a high temperature elemental analyzer

Oxygen isotope intratooth profiles
Using a serial sampling approach, series of time-dependent 18 Figure S8). As not all teeth preserve full annual cycles, not every tooth could yield both summer and winter seasonal 18 O information, which leads to differences in the number of total data points extracted for summer and winter 18

Inverse modeling of intratooth oxygen stable isotope profiles
Sequential sampling of herbivore tooth enamel necessarily introduces a certain amount of time averaging of the 18 O signal. This is due to the nature of enamel formation (amelogenesis) and the sampling strategy. Due to the complex and extended multi-phase process of enamel formation and maturation, where most of the mineral mass -and therefore the dominant isotopic inputs into tooth enamel -is added in a diffuse process that does not follow the visible incremental geometry of enamel apposition, but is continuously averaged to a varying extent over the enamel maturation process.

Oxygen isotope analyses as a temperature proxy
While raw 18  The secondary assumption is then, that surface waters consumed by Bos/Bison are broadly reflective of 18 Oprecip, which is the case for the majority but not all surface waters. 55 Notable exceptions are deep groundwaters, large rivers and large lakes, which can be isotopically decoupled from local precipitation due to effects of water transport, water residence time and evaporation. [56][57][58] As there are no large lakes present in the area of the site this leaves three potential scenarios of water source choice could lead to divergence of 18  where strontium isotope analysis of river samples and the regional isoscape also suggest slightly higher 87  We can therefore not confidently exclude that Bos/Bison in this study were consuming water from the Vézère River. However, the impacts of river water consumption on oxygen isotope patterns would be quite limited in this case. As no direct oxygen isotope measurements are available for the Vézére river, we make an estimation of the maximum isotopic deviation from precipitation in the low elevation plateaus around the La Ferrassie locality. The waters with lowest 18 O that contribute to the Vézère River, will originate from the areas with the highest elevation in the source area of the river. Water with higher 18 O will flow into the river at later points, making it progressively higher in 18 18 O of precipitation in the plateau regions around the site of La Ferrassie. 61 Considering that this isotopic difference is of comparable magnitude only as large as typical intra-individual variability within one archaeological layer and is at least 3 times smaller than the temperature induced isotopic difference expected between a glacial phase and an interglacial phase in the Late Pleistocene, we regard this effect to be negligible.
The reconstruction of air temperature from precipitation 18 O in turn relies on a strong linear relationship between air temperature and precipitation 18 O that occurs mid to high latitude environments. [16][17][18][62][63][64] The exact relationship between 18 Oprecip and temperature can vary slightly depending on geographic location due to secondary influences of atmospheric circulation on 18 Oprecip, but the slopes of temperature to 18 Oprecip regressions are spatially remarkably stable with only small differences in slope between different mid-to high latitude regions. 55,64-68 At the same time, due to the influence of atmospheric circulation on the 18 Oprecip-temperature relationship, air temperature reconstruction based on modern calibration data sets therefore also assumes a certain degree of equivalency in atmospheric circulation between modern day calibration data sets and past environments. The robustness of the 18 Oprecip-temperature relationship between different regions with different atmospheric circulation regimes in modern day data sets suggests however, that the 18 Oprecip-temperature relationship is stable to moderate circulation changes. In this study we employ a modern calibration data set that encompasses data from a range of locations across Europe, and the variability in the 18 Oprecip-temperature relationship contributes to the uncertainty estimates given for the final paleotemperature reconstructions. 46 Our model therefore already incorporates a numerical representation of a certain degree of variability related to impacts of atmospheric circulation on the 18 Opreciptemperature relationship to the extent that it is present across modern day localities.
Although Pleistocene data on the 18 45 Given this data we conclude that the assumption that modern day 18 Oprecip-temperature relationships can be applied to the Late Pleistocene is sufficiently met.
In addition to potential effects on 18 Oprecip from differences in atmospheric circulation, there can also be some isotopic impact from changes in sea water 18 O. Sea water 18 O changes between glacial and interglacial phases due to the differences in global ice volume. 72 While the isotopic composition of precipitation is predominantly driven by circulation dynamics and precipitation conditions, the isotopic composition of sea water is does also affect 18 Oprecip, as sea water forms the original source of water vapor that produces clouds and eventually precipitation. 55

Control for migratory behavior using 87 Sr/ 86 Sr
To confirm that samples Bos/Bison did not exhibit long distance migratory behavior, seven individuals were sampled for Sr isotope analysis. For each of the seven individuals two data points from different seasons (as determined from 18 O) were obtained to assess whether seasonal movements were undertaken by these animals. All data points can be found in Supplementary

Correlations between different isotopic systems
To investigate how different environmental and ecological changes may been connected throughout the sequence of La Ferrassie, and to explore the most likely major drivers of different isotopic systems in this case study, we conduct correlation tests of the different isotopic systems studied here (carbon, nitrogen, oxygen). A comparatively good correlation can be observed between diachronic changes in 13

Sample provenance and spatial patterns in the stratigraphy
To exclude that the mismatch between climatic conditions indicated by Bos/Bison oxygen isotope values and the cold climate features in the sediments of Layer 2 is the result of a mixing of several occupations with differing climatic conditions within the same layer, we explore the spatial relationship between our sampled teeth, their oxygen stable isotope value, the location of reindeer bones in the sequence and the location of cold climate indicator in the layer sediments. A spatial plot of reindeer bones within the sequence reveals that they are spread evenly across Layer 2 (as well as the other studied layers) and can be found also in close proximity to Bos/Bison teeth that were sampled for this study (Supplementary Figure S15). However, the number of bone fragments that could be specifically identified to Rangifer tarandus is low in all layers studied here, so we additionally also plot bone fragments identified as Cervid/Rangifer, which contain a mixture of large cervids and rangifer fragments, which show the same trend as the bone fragments specifically identified to Rangifer tarandus (Supplementary Figure S15).
Supplementary Figure Figure S17). However, this does not appear to be the case as a samples with lower (I4-792) and higher (I4-785) 18  As we sample both teeth and bones from Bos/Bison that are not anatomically connected to each other, we explore the spatial relationship between tooth and bone samples to establish that both sample types represent comparable groups of animals and can be interpreted as evidence of climatic conditions for the same time period. For all three layer where teeth were analyzed, at least one tooth was recovered in close spatial proximity to a bone sample (Supplementary Figure S18). Additionally, bone samples originate from several different parts of each layer and should therefore not be substantially biased towards any particular group of faunal remains within any layer. The bone samples that are in close proximity to tooth samples are likely to represent the same populations as these teeth. The lack of isotopic variability within both bones and teeth then further indicates that this connection also extends to the tooth samples that are not themselves in proximity to a bone sample.
Supplementary Figure S18:Spatial projection of all teeth (circles) and bones (diamonds) sampled for stable isotope analysis show that for each layer at least one tooth sample was obtained from close proximity to a bone sample. Tooth samples are labelled with Sample IDs and all sample points are colored by layer. Shown is the west section of the Main area of the western sector, with the North wall face of the Abri on the right of the image. Note that spatial projections of samples can plot slightly past layer boundaries that can be seen in the section orthophotograph as some samples were recovered at some distance from the section. One bone sample (F7-107; Layer 3) was excluded from this projection as it originates from a greater distance from the section and plotted far outside the layer boundaries of the section photograph.
Additionally, we plot the location of tooth samples in relation to all recorded Bos/Bison bones larger than 2.5 cm in Supplementary Figure S19. This plot shows that Bos/Bison bones are evenly spread throughout Layers 2, 5A and 5B, including all areas where teeth were obtained.
We therefore overall believe that -to the extent that can be determined by spatial association -the Bos/Bison bones and teeth from La Ferrassie are sufficiently well connected to each other to assume they represent a comparable collection of animals.
Supplementary Figure S19:Spatial plot of the find location of Bos/Bison teeth sampled for oxygen isotope analysis (large orange points, labelled by Sample ID, and colored by Layer) in relation to all recorded Bos/Bison bones > 2.5 cm (small red circles). All bone fragments > 2.5 cm (small black circles) are plotted for context in the background. Shown is the west section of the Main area of the western sector, with the North wall face of the Abri on the right of the image.

SI Tables
Supplementary