Late Pleistocene Neanderthal exploitation of stable and mosaic ecosystems in northern Iberia shown by multi-isotope evidence

Abstract

During the last glacial period, rapidly changing environments posed substantial challenges to Neanderthal populations in Europe. Southern continental regions, such as Iberia, have been proposed as important climatic “buffer” zones during glacial phases. Contextualising the climatic and ecological conditions Neanderthals faced is relevant to interpreting their resilience. However, records of the environments and ecosystems they exploited across Iberia exhibit temporal and spatial gaps in coverage. Here we provide new evidence for palaeotemperatures, vegetation structure, and prey herbivore ecology during the late Pleistocene (MIS 5–3) in northern Spain, by applying multiple stable isotope tracers (δ¹⁸O, δ¹³C, δ¹⁵N, δ³⁴S) to herbivore skeletal remains associated with Neanderthal occupations at Axlor Cave, Bizkaia. The results show little change over time and indicate stable climatic conditions and ecosystems across different occupations. Large within-layer isotopic variability in nitrogen and sulphur suggests the presence of a mosaic environment and a variety of isotopic ecotones that were exploited by Neanderthals and their prey. We implement a combination of carbonate and phosphate δ¹⁸O measurements to estimate palaeotemperatures using a cost-effective workflow. We show that the targeted use of phosphate δ¹⁸O measurements to anchor summer peak and winter trough areas enables high-precision seasonal palaeoclimatic reconstructions.

INTRODUCTION

How Neanderthals exploited and adapted to rapidly changing environments during the late Pleistocene is relevant for understanding Neanderthal ecology and of particular importance during the late Marine Isotope Stage (MIS) 3, when they disappeared from the continent. Establishing the ecological niche and adaptability of Neanderthal groups during this time is especially relevant for exploring ecological and evolutionary differences with Homo sapiens to elucidate the contrasting fate of these two species. It has frequently been suggested that Atlantic and Mediterranean coastal areas provided a relatively mild climate and a stable species-rich ecosystem that buffered a number of late Pleistocene plant and mammal species, including Neanderthals, against harsher climates (O'Regan, 2008; González-Sampériz et al., 2010; Jennings et al., 2011; Bradtmöller et al., 2012; López-García et al., 2012; Sánchez Yustos and Diez Martín, 2015; Bicho and Carvalho, 2022). There is abundant evidence of hominin life in this period and region, particularly in the northern Atlantic zone of Iberia (in the areas of modern-day Asturias, Cantabria, and the Basque Country), including several MIS 3 Middle Palaeolithic deposits such as La Güelga, El Castillo, Esquilleu, Morín, El Cuco, Amalda I, and Axlor (Maroto et al., 2012; Higham et al., 2014; Straus and González Morales, 2016; Marín-Arroyo et al., 2018) as well as a small number of assemblages that are consistent with the Châtelperronian technocomplex (Ríos-Garaizar et al., 2022). Extensive archaeological and archaeozoological assemblages have been excavated from many of these sites (Ríos-Garaizar et al., 2012; Yravedra-Sainz de los Terreros et al., 2016; Rasilla et al., 2020; Álvarez-Vena et al., 2021), and key lines of research into late Neanderthal subsistence, behaviour, and resilience are ongoing in the region (Gómez-Olivencia et al., 2018; Marín-Arroyo et al., 2020; Marín-Arroyo and Sanz-Royo, 2022; Vidal-Cordasco et al., 2022).

Despite the abundant record of hominin activities, palaeoclimatic records to support the notion of contemporary climatic stability and mildness are predominantly far removed from archaeological sites (e.g., marine cores in the Bay of Biscay and at the Iberian margin; Roucoux et al., 2005; Martrat et al., 2007; Sánchez Goñi et al., 2008; Moreno et al., 2014). Furthermore, terrestrial records and local reconstructions from archaeological deposits are not abundant and provide inconsistent coverage of the peninsula and are particularly sparse in the region surrounding modern-day Cantabria and the Basque Country (Moreno et al., 2014; Fernández-García et al. 2023).

In recent years, substantial efforts have been made to reconstruct the chronology of archaeological sites and to better understand the local palaeoenvironmental conditions Neanderthals faced (Wood et al., 2014, 2018; Jones et al., 2018, 2020; Marín-Arroyo et al., 2018). One established way of generating information on landscape-scale changes in past climates, environments, and ecosystems is the analysis of stable isotope ratios of macromammal remains from archaeozoological and palaeontological assemblages. Particularly in the case of anthropogenically accumulated faunal assemblages, a direct association with the archaeological record, and therefore, traces of hominin behaviour, can be established. A variety of isotopic tracers can be applied to skeletal faunal remains to yield a wealth of palaeoclimatic, palaeoenvironmental, and palaeoecological information (Drucker et al., 2003, 2011; Stevens and Hedges, 2004; Stevens et al., 2008, 2014; Bocherens et al., 2014; Nehlich, 2015; Jones et al., 2018, 2019, 2020; Pederzani and Britton, 2019; Reade et al., 2020). Stable isotope studies of faunal remains, specifically using carbon and nitrogen stable isotope analyses, have already started to provide valuable information concerning the local environmental conditions experienced by late Neanderthals and Upper Palaeolithic Homo sapiens in the Cantabrian region (Jones et al., 2018, 2019), indicating, for example, a tentative shift to a more open steppe environment with a niche separation between horses (Equus ferus) and red deer (Cervus elaphus) after the end of Middle Palaeolithic and initial Aurignacian (Jones et al., 2018, 2020; Rodríguez-Almagro et al., 2021). Such palaeoecological information has recently increasingly been supplemented by information on spatial ecology from sulphur stable isotope analysis (Jones et al., 2019; Reade et al., 2020). In these examples, addition of sulphur stable isotope analysis as a spatially informative tracer has been very effective in teasing apart dietary and spatial niche separation between different species. Together, these environmental and spatial isotopic tracers effectively track landscape and palaeoecological changes and the presence of isotopic microenvironments, whereas identifying clear climatic trends has been more challenging (Jones et al., 2019).

Larger data sets from more strongly climatically driven stable isotope systems, such as δ ¹⁸O, are needed to better understand the palaeoclimatic context for the late Pleistocene on the northern Spanish coast and to reconstruct palaeotemperatures. Temperature conditions, particularly seasonal temperature variations, are important drivers of plant and animal ecosystems and shape landscapes that hominins inhabit, the subsistence resources that are available for their consumption, and the conditions they need to shelter from. The existing stable isotopic proxies indicate a complex ecosystem with an extensive range of different microenvironments coexisting within this particular orographic region (Jones et al., 2018, 2019; Marín-Arroyo et al., 2020). Such strong spatial gradients in the Cantabrian environment necessitate more local climatic information to disentangle how the northern Iberian coast may have acted as an area that sheltered populations of late Pleistocene mammal species, including Neanderthals, from the rapid and dramatic climatic changes of this time period.

Axlor is a particularly suitable site to reveal the climatic and environmental conditions Neanderthals faced from MIS 5 to early MIS 3 in northern Iberia, as it preserves a long Middle Palaeolithic sequence containing lithics, fauna, and hominin skeletal remains. Using ungulate remains from the anthropogenically accumulated faunal assemblage, we apply an extensive and complementary suite of stable isotopic analyses (δ ¹⁸O, δ ¹³C, δ ¹⁵N, δ ³⁴S) to provide new insights into local palaeoclimate, including seasonal palaeotemperatures—a first for the Middle Palaeolithic of this region. We employ an innovative combination of carbonate and phosphate oxygen stable isotope analyses of sequentially sampled tooth enamel, in which carbonate measurements are used to first identify intra-tooth δ ¹⁸O peaks and troughs, followed by phosphate oxygen isotope analysis of a limited sample selection for palaeotemperature reconstructions. Not only does our approach allow the reconstruction of seasonal palaeotemperatures in a streamlined and cost-effective way, but it also allows us to explore the amount of uncertainty introduced from converting δ ¹⁸O values between measurements of the different fractions, permitting us to provide recommendations for optimal extraction of palaeotemperature information from combined carbonate/phosphate oxygen isotope time series. In parallel, δ ¹³C, δ ¹⁵N, and δ ³⁴S are measured on red deer, horse, and large bovine (aurochs/bison [Bos primigenius/Bison priscus]) bone elements to extract information on landscape and vegetation structure and herbivore feeding ecology (δ ¹³C, δ ¹⁵N), as well as herbivore spatial ecology (δ ³⁴S) associated with the late Pleistocene Neanderthal occupations at the site.

THE SITE OF AXLOR

Axlor Cave is located in Dima municipality, Bizkaia Province, in the Basque Autonomous Community located in northern Iberia (Fig. 1). The region is characterised by a narrow coastal band, bounded by the Cantabrian Mountains to the south, with mountain ridges often extending up to the shore. The valleys in this coastal corridor are typically deep and narrow, running perpendicular to the coastline. Bizkaia is ecologically and topographically distinct from the western parts of the Cantabrian region and is situated between the Atlantic coast and a relatively low section of the Cantabrian Cordillera, with winding, narrow, steep-sided valleys, creating a structurally complex, “banded” topography. Axlor is located close to one of the lowest mountain passes linking the Cantabrian basins and the Alavese Plateau (43°7.27′N, 2°43.68′W; see Fig. 1), at an elevation of approximately 320 m above sea level (m asl), on the southwestern slope of the Dima valley. The area is characterised by rugged, mountainous terrain (Ríos-Garaizar and Moreno, 2015). Several streams and rivers exist in the vicinity of the site, which is relevant for the interpretation of oxygen stable isotope dynamics. The valley bottom is formed by the stream of the Indusi River, which originates in the nearby Altungane (765 m asl) and Saibigain (800 m asl) mountains (Basaguren and Orive, 1991). An additional small stream runs close to the site, which later merges into the Indusi River (Fig. 1). To the southwest, the Dima valley is bordered by the Sierra de Ugatxa (Urragiko Atxa, 588 m asl), which forms the separation from the parallel valley containing the Arratia River (Basaguren and Orive, 1991). The two valleys, as well as the Indusi and Arratia Rivers, merge at a point ~9 km northwest of Axlor. These rivers form part of the watershed of the Ibaizabal River (Ocio et al., 2016).

Figure 1. (A) Location of Axlor in the Cantabrian region, northern Spain. (B) Hydrotopographic setting of the area surrounding Axlor. Elevation contour lines are spaced 200 m apart. Elevation data from the European Digital Elevation Model v. 1.1 (European Environment Agency, 2016). Site location and waterways imported from OpenStreetMap (OpenStreetMap Contributors, 2020).

Axlor preserves an important Middle Palaeolithic sequence originally excavated by Barandiarán from 1967 to 1974; more recently from 2000 to 2008 by González-Urquijo, Ríos-Garaizar, and Ibáñez-Estévez; and since 2019, by González-Urquijo and Lazuén. These new excavations have yielded a stratigraphic sequence overlapping with, but not identical to, the one described by Barandiarán (1980; for further discussion, see Gómez-Olivencia et al. 2018; Table 1). The material used in this study stems from the faunal collection of the Barandiarán excavation curated at the Biscay Archaeological Museum in Bilbao. In this early excavation, nine layers (I–IX) were identified. Layers III to VIII were attributed to the Middle Palaeolithic based on the lithic record, but with some technological differences between the upper and lower parts of this sequence. In the upper part of the Middle Palaeolithic sequence (Layers III–V), the lithic record has been described as consistent with the Quina Mousterian technocomplex, while the lithic record of the lower sequence (Layers VI–VIII) is dominated by Levallois blank production (Ríos-Garaizar, 2017). The upper part of the Middle Palaeolithic sequence also contains a substantial number of bone tools (Mozota Holgueras, 2012).

Table 1. Overview of the Middle Palaeolithic layers identified at Axlor, including published information on their ages.

^a Only samples marked with UF-AMS (ultrafiltration and accelerator mass spectrometry) have been produced on ultrafiltered collagen. Sample Beta-144,262 was initially attributed to Layer D (González-Urquijo et al., 2002), but has recently been included in Layer B (González-Urquijo et al., 2021).

The faunal assemblage appears to be predominantly of anthropogenic origin, with very little evidence of carnivore activity. In both the Barandiarán faunal collection (Altuna, 1989) and the collection from the recent re-excavations (Castaños, 2005), a change in the relative abundance of prey species is evident. In Layers III–V (~MIS 4 to early MIS 3) the most abundant taxa, to an approximately equal extent, are red deer (Cervus elaphus), aurochs/bison (Bos primigenius/Bison priscus), and Iberian ibex (Capra pyrenaica) (Altuna, 1989; Castaños, 2005; Gómez-Olivencia et al., 2018). In Layer VII (possible MIS 5c) the faunal spectrum is dominated by Capra pyrenaica and in Layer VIII (possible MIS 5c) a substantial proportion of the faunal remains stem from Cervus elaphus. A dominance of Capra pyrenaica could indicate cold-arid conditions, as it has been shown in northern Iberian examples that they may have been more resistant to arid conditions or phases of low primary productivity than other taxa such as large bovines (Jones et al., 2020, 2021; but see Yravedra and Cobo-Sánchez, 2015). It should be noted that the Barandiarán faunal collection exhibits significant collection bias in favour of morphologically identifiable fragments and likely also mixes to some extent the material of different deposits due to excavation in artificial spits identified lying horizontally (Gómez-Olivencia et al., 2018), but species abundances in the material here are very similar to those in the more recent excavations (Castaños, 2005) and appear not significantly affected by this collection bias.

Radiocarbon dates obtained from ultrafiltered collagen extracts from Layer IV yielded infinite ages, making this layer and the sequence below (Layers V–IX) older than 51 ka BP (before 1950), but otherwise of unknown age today (Marín-Arroyo et al., 2018; Table 1). It should be noted that many of the radiocarbon dates that have been obtained from the site, including the finite dates from Layers B and D, represent minimum ages. Three of four dates obtained from Layer IV (corresponding Layers B and D) from the 2000–2008 excavations returned infinite radiocarbon ages (Table 1). Two of the infinite ages have minimum age cutoffs that are 3000 yr older than the date recently obtained from Layer III. The only finite date from Layer D (Beta-203,107: 44,920 ± 1950) calibrates to 54,805–44,624 cal yr BP, which does overlap to some extent with the new date for Layer III. However, the confidence in the Beta dates is undermined by the fact that most of the series is performed by accelerator mass spectrometry (AMS) without ultrafiltration (Gómez-Olivencia et al., 2020). In addition, those dates lack information on provenance within the site. Therefore, those AMS Beta dates performed without the ultrafiltration protocol must be considered with caution. When dating materials near the limit of radiocarbon dating, if there is a discrepancy between the results, the older dates are generally considered more accurate, assuming that contaminants were not properly removed from the younger ones (Higham, 2011). Accordingly, our examination of chronological hygiene favours the Oxford Radiocarbon Accelerator Unit (ORAU) ultrafiltration bone ¹⁴C ages for final dating assessments owing to their exact provenance, completeness of the bone collagen suitability indicators, and the more conservative results and uncertainties obtained for replicate samples in Axlor. Based on the abundance of red deer in the lower layers of the Middle Palaeolithic sequence and the archaeological similarities in stone tool technology and fire structures between Axlor Layers VI and VII and Arrillor's Medium Complex (Amk and Smkl), it has been suggested that Layers VI–VIII can be attributed to early MIS 3/MIS 4 and MIS 4 (Ríos-Garaizar, 2017), but new data suggest an older attribution for these levels to ~60 ka BP and an attribution of the lower levels VII and VIII to MIS 5c (~100 ka BP) (Sánchez Hernández, 2021).

A small number of hominin fossil remains were recovered during excavations, some of which were published soon after the completion of the Barandiarán excavations (Basabe, 1973), while others were only identified as hominin fossil remains during later reanalyses of the collections (Gómez-Olivencia et al., 2020). These recent analyses identified two teeth stemming from Layers IV and V as belonging to Neanderthals and also confirmed the taxonomic and stratigraphic identification of a Neanderthal cranial fragment from Layer VIII identified during the Barandiarán excavations (Gómez-Olivencia et al., 2020). Other hominin remains, described initially by Basabe (1973), have recently come under debate, with some researchers suggesting a Homo sapiens taxonomic identification and an origin from the upper layers of the stratigraphy that may have been redeposited due to earlier disturbances of the site during the nineteenth and early twentieth centuries (Gómez-Olivencia et al., 2020, 2023), while others support the original assignment to Neanderthals and argue against any influence of disturbance (González-Urquijo et al., 2021). A recent reply to the latter marks the end of this debate (Gómez-Olivencia et al., 2023), which will likely only be definitively solved using biomolecular evidence and direct dating.

ISOTOPE ARCHAEOZOOLOGY APPROACHES TO PAST PALAEOENVIRONMENTS

Oxygen isotope analyses as a proxy of seasonal climate

The use of the oxygen stable isotope composition (δ ¹⁸O) of faunal tooth enamel as a proxy of past climates rests on two underlying relationships between δ ¹⁸O of tooth enamel (δ ¹⁸O_enamel) and δ ¹⁸O of animal drinking water (δ ¹⁸O_dw) and between δ ¹⁸O of environmental water (δ ¹⁸O_ew) and air temperature, which is a dominant effect in mid- to high latitudes (see below in this paragraph) (Longinelli, 1984; Luz et al., 1984; Rozanski et al., 1993; Iacumin et al., 1996; Kohn et al., 1996; Clark and Fritz, 1997; Müller et al., 2017). Tooth enamel is formed in isotopic equilibrium with animal body water (Longinelli, 1984; Luz et al., 1984; Iacumin et al., 1996; Daux et al., 2008), which in turn mirrors variability in δ ¹⁸O of drinking water for large-bodied, obligate drinking terrestrial mammals such as horses, large bovines, rhinocerotids, or probiscideans, due to the large amounts of drinking water consumed by these animals (D'Angela and Longinelli, 1990; Bryant et al., 1994, 1996a; Sánchez Chillón et al., 1994; Hoppe et al., 2004; Hoppe, 2006). Consequently, δ ¹⁸O_enamel of these animals empirically shows a strong linear correlation with δ ¹⁸O of drinking water (Longinelli, 1984; Luz et al., 1984; Rozanski et al., 1992; Hoppe, 2006; Tütken et al., 2007; Arppe and Karhu, 2010; Pryor et al., 2014; Skrzypek et al., 2016). The isotopic composition of surface waters that are available for drinking is determined by several complex effects, but the oxygen isotopic composition of precipitation (δ ¹⁸O_precip) at mid- to high-latitude locations in temperate climate zones generally shows a strong dependence on air temperature within a single location, and this translates into precipitation-fed surface waters (Dansgaard, 1964; Rozanski et al., 1993; Clark and Fritz, 1997; Kohn and Welker, 2005; Gat, 2010; Stumpp et al., 2014; Müller et al., 2017). Through this effect, the role of δ ¹⁸O_enamel as a proxy for palaeotemperatures at mid- to high-latitude sites is derived, where higher δ ¹⁸O in tooth enamel indicates higher air temperatures, while low δ ¹⁸O values indicate lower temperatures. This is true on longer timescales of climatic change as well as on subannual timescales (Stumpp et al., 2014; Müller et al., 2017), and such subannual δ ¹⁸O change is recorded as a sinusoidal seasonal pattern of δ ¹⁸O in sequential tooth enamel samples of animal teeth, with peaks indicating the summer season and troughs the winter season (Fricke and O'Neil, 1996; Fricke et al., 1998). Such an extraction of a time-dependent measurement series is possible due to the incremental formation of tooth enamel from the crown to the enamel–root junction (ERJ) with no remodelling of the tissue after mineralisation has been completed (Dean, 1987; Hillson, 2005). The stable isotope measurement series obtained from one tooth, therefore, covers the time of tooth formation, which is usually several months to over a year in Bos/Bison sp. (Brown et al., 1960). The use of tooth enamel δ ¹⁸O as a palaeotemperature proxy can be formalised for palaeotemperature estimation by establishing species- and site-specific relationships between δ ¹⁸O_enamel and δ ¹⁸O_dw and between δ ¹⁸O_precip and air temperature via linear regression of modern calibration data (Tütken et al., 2007; Arppe and Karhu, 2010; Skrzypek et al., 2011, 2016; Pryor et al., 2014). It is assumed that the linear relationships established in this way are applicable to the late Pleistocene, as other factors such as animal metabolism, physiology, and drinking requirements are thought to be stable through time (see Supplementary Text 1). For the relationship between δ ¹⁸O_precip and air temperature, a certain degree of stability in atmospheric circulation is assumed, which, based on circulation models and measurements of old groundwater, seems sufficiently well supported for the timescales covered here (see Supplementary Text 1). Other assumptions, such as a rough isotopic equivalency between animal drinking water and precipitation should be validated on a case-by-case basis. A more detailed breakdown of these assumptions can be found in Supplementary Text 1.

Bioapatite phosphate and carbonate groups as isotope ratio mass spectrometry (IRMS) analytes

Many palaeotemperature studies focus on the oxygen-bearing phosphate (PO₄³⁻) component of tooth enamel (a biogenic carbonate–substituted hydroxyapatite, or “bioapatite”), but the carbonate (CO₃²⁻) moiety presents another oxygen-bearing group that is commonly analysed for δ ¹⁸O. The oxygen isotopic composition of these different moieties reflects the same environmental proxy information, but each has specific advantages and disadvantages in regard to stable isotope analysis. While analysis of carbonates is generally simpler and comparatively inexpensive, the analysis of the more laborious to analyse phosphate component offers important advantages in preservation and resistance to diagenesis (Iacumin et al., 1996; Zazzo et al., 2004; Lee-Thorp, 2008). Analytical precision can be significantly better for carbonate analyses (between ~0.05‰ and 0.1‰ for dual-inlet IRMS and an automated carbonate preparation device setup) but is practically often similar to the precision achieved for phosphate analyses (0.2–0.3‰) due to widespread use of continuous-flow (CF) IRMS with GasBench peripherals. Analyses of bioapatite phosphate, furthermore, may have an advantage in the context of palaeotemperature reconstruction specifically, as the vast majority of modern calibration data has so far been generated from this analytical fraction (Pederzani and Britton, 2019). The use of bioapatite carbonate δ ¹⁸O (δ ¹⁸O_carb) for palaeotemperature estimation therefore usually requires an additional conversion step to predicted δ ¹⁸O of bioapatite phosphate (δ ¹⁸O_phos), which incurs an increased estimation uncertainty that may offset the lower analytical precision of some δ ¹⁸O_carb analyses (Skrzypek et al., 2016).

Due to the complementary advantages of each fraction, a combination of both types of analysis is desirable. A conversion between δ ¹⁸O_carb and δ ¹⁸O_phos is possible, as these measures show a strong linear relationship with a slope close to 1 and an intercept offset (δ ¹⁸O_carb-phos) of approximately 9‰ that shows little species-specific variation (Longinelli and Nuti, 1973; Bryant et al., 1996b; Iacumin et al., 1996; Martin et al., 2008; Pellegrini et al., 2011; Chenery et al., 2012; Trayler and Kohn, 2016). This linear relationship can be used to convert between measurements of the different moieties. While δ ¹⁸O_carb-phos is relatively consistent in empirical studies and matches well with modelled values (Aufort et al., 2017), a certain degree of interindividual variability causes a prediction uncertainty for converting measurements that is not insignificant. Based on a large collection of empirical data, Pellegrini et al. (2011) showed that δ ¹⁸O_carb-phos can vary by a magnitude of approximately 1.5‰. For this reason, the use of the spacing between the fractions as a diagenetic indicator on the level of individual samples has been discouraged by some (Pellegrini et al., 2011; France and Owsley, 2015). In contrast to this, Pryor et al. (2014) argued that in the context of palaeotemperature estimation, the error introduced from this conversion is insignificant compared with the calibration errors of the temperature conversion. However, the propagation of this error through the additional temperature calibration steps was not quantitatively considered in Pryor et al. (2014), and the particular concerns related to sequentially sampled δ ¹⁸O time-series data were not discussed. To explore this further, we specifically investigate how the uncertainty introduced from converting δ ¹⁸O_carb and δ ¹⁸O_phos affects seasonal palaeotemperature estimates in a workflow that includes inverse modelling of oxygen isotopic time series (e.g., Passey et al., 2005) before temperature estimation.

Carbon, nitrogen, and sulphur isotopes

The carbon (δ ¹³C) and nitrogen (δ ¹⁵N) isotopic composition of animal bone collagen reflects the isotopic composition of the animal's diet with approximately fixed offsets of ~1‰ for carbon and 3–5‰ for nitrogen on a plant diet (DeNiro and Epstein, 1978, 1981). For herbivores, this means that these systems reflect isotopic dynamics in the plant biome, where δ ¹³C varies with plant photosynthetic pathway (C₃, C₄, or CAM); plant functional type (graze vs. browse, woody plants vs. herbaceous plants); vegetation density; and climatic and atmospheric factors such as temperature, aridity, and atmospheric CO₂ concentration (Park and Epstein, 1961; Tieszen and Boutton, 1989; Merwe and Medina, 1991; Tieszen, 1991). Plant δ ¹⁵N is driven by climatic factors such as aridity/rainfall and temperature, as well as nitrogen isotopic effects in soils related to soil moisture, bacterial activity, and nutrient availability, which in the late Pleistocene appear to be partially driven by the presence or absence of permafrost (Amundson, 2003; Stevens and Hedges, 2004; Stevens et al., 2008, 2014; Drucker et al., 2011; Craine et al., 2015). These factors create isotopic differences between different plant communities, and through this herbivore δ ¹³C and δ ¹⁵N act as proxies for dietary specialisation. Moreover, the climatic factors influencing δ ¹³C and δ ¹⁵N, as well as the climatic tolerances of different plants and the resulting relationship between landscape change and climate change, mean that plant and herbivore δ ¹³C and δ ¹⁵N also track underlying changes in climate and landscape structure.

The sulphur isotopic composition (δ ³⁴S) of herbivore collagen reflects plant δ ³⁴S with little biological fractionation and, through this, underlying isotopic variability in the geosphere and hydrosphere. Differences in δ ³⁴S exist between different bedrock types and soils, which are often strongly driven by the isotopic differentiation between sulphur sources of marine and nonmarine origin (Thode, 1991; Nehlich, 2015). Furthermore, wetland plants can show very distinct δ ³⁴S signatures due to their uptake of sulphides that accumulate in anoxic growing conditions (Guiry et al., 2021). These dynamics create spatial differences in δ ³⁴S of bioavailable sulphur, enabling the use of δ ³⁴S as a movement or location tracer, particularly in coastal regions (McArdle et al., 1998; Nehlich, 2015; Bataille et al., 2021). However, spatial patterns of bioavailable sulphur, and therefore spatial variability of δ ³⁴S of plants, are complex and do not exactly mirror those in bedrock. At the same time, modern anthropogenic impacts (e.g., from fertilizers) make it additionally difficult to establish spatial variability in δ ³⁴S in a way that is useful for archaeological applications (Mizota and Sasaki, 1996; Szpak et al., 2019; Bataille et al., 2021). At the same time, an impact of permafrost presence and climatically induced changes in baseline δ ³⁴S have been proposed (Reade et al., 2020). Specific assignment of spatial ecologies based on δ ³⁴S therefore remains experimental but can be useful in indicating changes in spatial ecology or an overlap (or lack thereof) in the habitats used by different species within a single food web (e.g., Britton et al., 2023). This is particularly the case in coastal environments, where larger δ ³⁴S gradients are to be expected.

MATERIALS AND METHODS

Materials and sampling methodology

Multiple isotopic approaches were applied to herbivore skeletal remains, including bones and teeth, to obtain climatic, environmental, and ecological information for the Middle Palaeolithic landscapes around Axlor. Faunal specimens were sampled from the collections of the Arkeologi Museoa Bilbao, where the material is curated. All specimens were photographed before and after sampling. A total of seven suitable aurochs/bison (Bos primigenius/Bison priscus) teeth from Layers III, IV, and VI were chosen for oxygen stable isotope analysis of sequentially sampled tooth enamel from the entire sequence (Table 2). The sample obtained here represents a mixture of lower second (M2) and third molars (M3) (Table 2). Sequential samples were drilled from the occlusal surface to the ERJ along the full crown height. Each sample was removed as a small strip of ~1 mm × 5 mm dimension perpendicular to the tooth growth using a diamond tipped drill bit. Drill bits were cleaned with ~1 M nitric acid and two cycles of ultrasonication in Milli-Q ultra-pure water (5 minutes each) to prevent cross-contamination between samples. A total of 106 sequential samples were obtained and analysed for δ ¹⁸O_carb. Summer peak and winter trough areas of the δ ¹⁸O time series were identified by visual examination, and three samples each per summer peak and winter trough (n = 47) were chosen to be analysed for δ ¹⁸O_phos as well.

Table 2. Contextual information of teeth sampled for oxygen isotope analysis of tooth enamel bioapatite carbonate and bioapatite phosphate.

^a MIS, Marine Isotope Stage.

Additionally, bone samples were taken from aurochs/bison, horses, and red deer for Layers III, IV, VI, and VIII (Table 3) to generate new carbon and nitrogen (n = 24) and sulphur (n = 65) stable isotope data from bone collagen. An existing data set of carbon and nitrogen stable isotope data from 42 bone samples from Axlor published in Jones et al. (2018) was expanded to a total sample size of 66 (Table 3). Sulphur stable isotope data were generated from the same collagen samples extracted for this study, and bone originally used by Jones et al. (2018) was resampled to generate sulphur stable isotope data for the existing carbon and nitrogen stable isotope data set. A tibia fragment from a large bovine (AXL55) from Layer III was additionally sent for radiocarbon dating at the ORAU, University of Oxford, UK.

Table 3. Number of data points used in this study for different isotopic analyses of bone collagen.

^a Carbon and nitrogen stable isotope data points included from Jones et al. (2018) are marked with an asterisk (*).

Stable isotope analysis of bioapatites (δ ¹⁸O_phos, δ ¹⁸O_carb, δ ¹³C_carb)

For carbon and oxygen stable isotope analysis of tooth enamel carbonates, approximately 5 mg of untreated tooth enamel powder was sent for analysis at IsoAnalytical (Crewe, UK), where samples were weighed into clean exetainer tubes, which were then flushed with 99.995% helium. Phosphoric acid was prepared following guidelines in Coplen et al. (1983) and was then injected through vial septa, and samples were left to react overnight to allow for complete conversion to CO₂. The resulting CO₂ was then analysed for δ ¹³C and δ ¹⁸O by CF-IRMS using a Europa Scientific 20–20 isotope ratio mass spectrometer. Stable isotope measurements were calibrated to the VPDB (Vienna Pee Dee Belemnite) scale and quality controlled using an in-house calcium carbonate standard (IA-R022, accepted values: δ ¹³C = −28.63‰ VPDB and δ ¹⁸O = −22.69‰ VPDB), the international carbonatite standard NBS-18 (δ ¹³C = −5.01‰ VPDB and δ ¹⁸O = −23.2‰ VPDB), and an in-house chalk standard (IA-R066, accepted values: δ ¹³C = 2.33‰ VPDB and δ ¹⁸O = −1.52‰ VPDB). The accepted values of the in-house standards IA-R022 and IA-R066 were obtained by calibrating against international reference materials. Values of IA-R022 were obtained by calibrating against NBS-18 and NBS-19 (TS-limestone, δ ¹³C = 1.95 ‰ VPDB and δ ¹⁸O = −2.20 ‰ VPDB), while values of IA-R066 were obtained using NBS-18 and IAEA-CO-1 (Carrara marble, δ ¹³C = 2.5‰ VPDB and δ ¹⁸O = −2.4‰ VPDB), respectively. A limestone in-house standard, ILC1, was used as an additional quality-control standard. Values obtained for this standard during the analyses of Axlor samples of δ ¹³C = 2.21 ± 0.09‰ VPDB (1 SD, n = 2) and δ ¹⁸O = −3.83‰ VPDB (1 SD, n = 2) compare well with the long-term mean values of this standard, which were δ ¹³C = 2.13‰ VPDB (1 SD, n = 2) and δ ¹⁸O = −3.99‰ VPDB. Analytical precision of quality-control standard replicates was better than 0.09‰ for δ ¹³C and better than 0.14‰ for δ ¹⁸O.

Tooth enamel powder samples were converted to silver phosphate for oxygen stable isotope analysis of bioapatite phosphate using digestion with hydrofluoric acid (HF), followed by crash precipitation of silver phosphates following an adapted version of the protocol developed by Dettman et al. (2001) and modified by Tütken et al. (2006), as described in Britton et al. (2015). A full description of this protocol can be found in Supplementary Text 2.

Oxygen isotope delta measurements of Ag₃PO₄ were conducted using a high-temperature elemental analyser (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). Approximately 0.5 mg of each silver phosphate sample was weighed into cleaned silver capsules (3 mm × 4 mm, IVA Analysentechnik, Meerbusch, Germany) and introduced to the TC/EA using a Costech Zero Blank Autosampler (Costech International, Cernusco sul Naviglio, Italy). Conversion to CO was achieved using a reactor temperature of 1450°C, and gases were separated using an Eurovector E11521 1.4 m × 4 mm × 6 mm stainless steel GC column with 80/100 mesh 5 Å molecular sieve packing (Eurovector Instruments & Software, Pavia, Italy) maintained at 120°C with a carrier gas pressure of 1.3 bar. Samples were measured in triplicate, but in rare cases, individual measurements were rejected if they did not meet quality-control criteria such as appropriate peak area to sample amount relationship. In such cases, δ ¹⁸O values are therefore only based on two replicates.

The δ ¹⁸O values were two-point scale normalised to the VSMOW (Vienna Standard Mean Ocean Water) scale using matrix-matched standards calibrated to international reference materials. Scale normalisation was checked using three separate quality-control standards. Scale normalisation was conducted using the B2207 silver phosphate standard (δ ¹⁸O = 21.7 ± 0.3‰, 1 SD; Elemental Microanalysis, Okehampton, UK) and an in-house silver phosphate standard (KDHP.N, δ ¹⁸O = 4.2 ± 0.3‰, 1 SD). This in-house standard was obtained by equilibrating a KH₂PO₄ solution with Leipzig winter precipitation at ca. 140°C for several days, after which the solution was neutralised using a small amount of NH₄OH and the phosphate precipitated as silver phosphate by addition of AgNO₃ solution. The accepted value of this in-house standard was determined by two-point calibration using B2207 and IAEA-SO-6 (barium sulphate, δ ¹⁸O = −11.35 ± 0.3‰, 1 SD; Brand et al., 2009).

Aliquots of an in-house modern cow enamel standard (BRWE, later replaced by BRWE.2 due to exhaustion of this material) and the standard material NIST SRM 120c (formerly NBS 120c) were precipitated and measured alongside each batch of samples to ensure equal treatment. Additionally, 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 normalisation independent of silver phosphate precipitation. Measurements of these standards gave δ ¹⁸O values of 15.2 ± 0.1‰ for BRWE (1 SD, n = 12), 14.7 ± 0.3 for BRWE.2 (1 SD, n = 10), 21.6 ± 0.7‰ for NIST SRM 120c (1 SD, n = 22), and 13.8 ± 0.3‰ for AS337382 (1 SD, n = 66). This compares well with the consensus value for NIST SRM 120c of 21.7‰ (Pucéat et al., 2010), as well as the long-term averages for BRWE of 15.2 ± 0.3‰, 14.5 ± 0.4 for BRWE.2, and 14.0 ± 0.3‰ for AS337382. Samples were usually measured in triplicate (as 0.5 mg aliquots), and average reproducibility of sample replicate measurements was 0.3‰. Consecutive analysis of sets of standards with widely spaced isotopic value showed no detectable memory effect, and consequently no memory effect correction was used. No effect of the blank, the sample amount, or peak height on the results was observed, and consequently no blank correction or linearity correction was used. If drift across the course of a run was detected in the quality-control standard AS337382, a linear drift correction based on the drift of both normalisation standards was used and checked with the quality-control standard.

Stable isotope analysis of bone collagen (δ ¹³C_coll, δ ¹⁵N_coll, δ ³⁴S_coll)

Typically, <1 g of bone was extracted from selected specimens using a low-vibration micromotor with a diamond-edge cutting wheel. Bones with evidence of anthropogenic modification (e.g., with evidence of fresh fractures or cut marks) were targeted to link specimens to periods of hominin activity. All sampled bones (but not teeth) exhibited such signs of anthropogenic modification. Samples for bone collagen δ ¹³C, δ ¹⁵N, and δ ³⁴S analysis were prepared using facilities in the Institute of Biomedicine and Biotechnology (IBBTEC) at the University of Cantabria in Spain. Collagen was extracted following the protocol of Richards and Hedges (1999) and is described in detail in Supplementary Text 3. Carbon, nitrogen, and sulphur stable isotope analyses of extracted bone collagen were undertaken at IsoAnalytical, in duplicate. Details of the IRMS instrumentation and operating parameters are described in Supplementary Text 3.

Carbon and nitrogen stable isotope measurements were calibrated to the VPDB and AIR (atmospheric N₂) scales and quality checked using in-house standards IA-R068 (soy protein, accepted values: δ ¹³C = −25.22‰ VPDB and δ ¹⁵N = 0.99‰ AIR), IA-R038 (l-alanine, accepted values: δ ¹³C = −24.99‰ VPDB and δ ¹⁵N = −0.65‰ AIR), IA-R069 (tuna protein, accepted values: δ ¹³C = −18.88‰ VPDB and δ ¹⁵N = 11.60‰ AIR), and a mixture of IAEA-C7 (oxalic acid, δ ¹³C = −14.48‰ VPDB) and in-house standard IA-R046 (ammonium sulphate, δ ¹⁵N = 22.04‰ AIR). The accepted values of in-house standards were obtained via calibration against international reference materials: IA-R068, IA-R038, and IA-R069 against IAEA-CH6 (sucrose, δ ¹³C = −10.449‰ VPDB) and IAEA-N-1 (ammonium sulphate, δ ¹⁵N = 0.40‰ AIR). During the analysis of Axlor samples, IA-R069 gave values of δ ¹³C = −18.90 ± 0.04‰ VPDB (1 SD, n = 6) and δ ¹⁵N = −0.65 ± 0.03‰ AIR (1 SD, n = 6). Average reproducibility of samples measured in duplicate was 0.04‰ for δ ¹³C and 0.01‰ for δ ¹⁵N.

For sulphur stable isotope analysis, calibration to the VCDT (Vienna Cañon Diablo Troilite) scale and correction for ¹⁸O contribution were conducted using in-house standards IA-R061 (barium sulphate, δ ³⁴S = 20.33‰ VCDT), IA-R025 (barium sulphate, δ ³⁴S = 8.53‰ VCDT), and IA-R026 (silver sulphide, δ ³⁴S = 3.96‰ VCDT). Accepted values for these in-house standards were obtained by calibrating against international reference materials NBS-17 (barium sulphide, δ ³⁴S = 5.25‰ VCDT) and IAEA-S-1 (silver sulphide, δ ³⁴S = −0.30‰ VCDT). IA-R061, IAEA-SO-5 (barium sulphate, δ ³⁴S = 20.3‰ VCDT), IA-R068 (soy protein, δ ³⁴S = 5.25‰ VCDT), and IA-R069 (tuna protein, δ ³⁴S = 18.91‰ VCDT) were used as quality-control standards. The accepted values of IA-R068 and IA-R069 were obtained via calibration using the international reference materials NBS-17 and IAEA-SO-5. During the analysis of Axlor samples, the following values were obtained for the quality-control standards: IA-R061 δ ³⁴S = 20.32 ± 0.08‰ VCDT (1 SD, n = 6), IAEA-SO-5 δ ³⁴S = 0.47 ± 0.07‰ VCDT (1 SD, n = 2), IA-068 δ ³⁴S = 5.17 ± 0.15‰ VCDT (1 SD, n = 4), and IA-R069 δ ³⁴S = 18.85 ± 0.10‰ VCDT (1 SD, n = 4). Average reproducibility of samples measured in duplicate was 0.2‰.

Radiocarbon dating

A bone sample from a Bos/Bison sp. tibia (AXL55) from Layer III was radiocarbon dated at the ORAU under the sample code OxA-40151. Both collagen extraction for carbon and nitrogen stable isotope analyses and radiocarbon dating for this sample were conducted at the ORAU (see Supplementary Text 4).

Seasonal palaeotemperature estimation

To determine preservation status and establish the linear relationship between the two analytes, δ ¹⁸O_phos and δ ¹⁸O_carb measurements made on the same sample were first checked against each other and against published data sets of the phosphate to carbonate δ ¹⁸O spacing (δ ¹⁸O_carb-phos). The resulting relationship was then used to predict δ ¹⁸O_phos for the samples with only δ ¹⁸O_carb available, and then predicted δ ¹⁸O_phos values were combined with direct δ ¹⁸O_phos measurements in the peak and trough areas to form an aggregated δ ¹⁸O_phos seasonal time series as the basis for palaeotemperature estimation (Fig. 2). For direct comparisons of seasonal δ ¹⁸O_phos values across archaeological layers and with other sites, direct δ ¹⁸O_phos measurements were used (Fig. 3).

Figure 2. (A) δ ¹⁸O_carb and δ ¹⁸O_phos show a strong linear correlation for Bos/Bison sp. samples from this study (blue) with a consistent isotopic spacing of around 9‰ that matches closely with previously reported data from modern enamel samples of large and medium non-human terrestrial mammals (orange; primarily bovines, equids and cervids; data from Bryant et al. [1996b]; Iacumin et al. [1996]; Shahack-Gross et al. [1999]; Zazzo et al. [2004]; Miller et al. [2019]; for discussion of lack of species differences, see Pellegrini et al. [2011]). Grey dashed lines represent the predictive interval of the linear model of the literature data. (B) The relationship between δ ¹⁸O_carb and δ ¹⁸O_phos was used to predict δ ¹⁸O_phos values from δ ¹⁸O_carb measurements and construct a combined δ ¹⁸O time series of δ ¹⁸O_phos measurements (opaque points) and predicted values (transparent points). ERJ, enamel–root junction.

Figure 3. Summer peak, mean annual, and winter trough δ ¹⁸O_phos data points extracted from Bos/Bison sp. seasonal δ ¹⁸O sinusoidal curves are similar across the Middle Palaeolithic sequence of Axlor. This indicates that climatic conditions were similar in these different site occupations. Summer peaks and winter trough points were identified by visual inspection of sinusoidal curves, shown in Fig. 2. Shaded areas depict the maximal spread of the data in each layer.

Before summer and winter palaeotemperature estimation, aggregated δ ¹⁸O_phos time series were put into an inverse model following Passey et al. (2005) to remove the damping of the seasonal δ ¹⁸O amplitude that is caused by time averaging through enamel mineralization and homogenization within each incremental sample. This approach estimates the original δ ¹⁸O input into tooth enamel by inverse modelling of the time averaging introduced through enamel mineralization and the sampling procedure by taking into account the sampling geometry and species-specific parameters of enamel formation. A detailed description of the modelling procedure with associated code can be found in Passey et al. (2005). It should be noted that this inverse model was originally developed for ever-growing teeth and does not take into account the slowing of enamel growth towards the ERJ seen in some molars of large ungulates, particularly horses (Zazzo et al., 2012; Bendrey et al., 2015). However, while such considerations have been used to develop updated inverse models for sheep (Green et al., 2018), no such model is currently available for cattle/bison, and several of the fundamental enamel mineralisation parameters necessary for such a model are unknown for these species. In cattle and bison, this effect is visible but not as pronounced as in other species, and mostly affects the late-growing portions of the third molar. Due to this effect, the amplitudes reconstructed here may represent minimum amplitudes; however, a systematic difference between the amplitudes reconstructed for second molars compared with third molars is not observed. Therefore, we estimate that this effect is relatively small in this study.

Enamel formation input parameters were chosen to reflect Bos/Bison sp. amelogenesis, following values given in Passey and Cerling (2002) and Kohn (2004). The initial mineral content of enamel was set at 25%, enamel appositional length as 1.5 mm, and maturation length as 25 mm. Additionally, sample input variables were given for distance between samples and sample depth. During the modelling procedure, a damping factor describing the damping of the isotopic profile amplitude needs to be chosen using an adjustment of a measured error term (E_meas) to the prediction error (E_pred). The damping factor in this study ranged between 0.01 and 0.06.

The most likely model solutions were used to extract corrected summer peak and winter trough δ ¹⁸O values. Modelled summer and winter values were extracted in the δ ¹⁸O profile locations that correspond to peaks and troughs in the original unmodelled δ ¹⁸O profile. The extracted modelled summer peak and winter trough values were used as the input to generate summer and winter δ ¹⁸O_dw values as well as summer and winter absolute temperatures. Mean annual δ ¹⁸O values were computed directly from raw δ ¹⁸O_phos values, rather than from model corrected curves. If a systematic sampling geometry is applied, the inverse model essentially symmetrically extends the amplitude of the sinusoidal δ ¹⁸O curve, making little changes to the average δ ¹⁸O of an annual cycle (for a detailed discussion, see Pederzani et al., 2021). Therefore, we used annual averages directly from δ ¹⁸O_phos values to avoid the impact of uncertainty on the model output inflating the error ranges on δ ¹⁸O_dw and palaeotemperature estimates.

Palaeotemperatures were then computed following methods outlined in Pryor et al. (2014) using the Excel files published therein for computations of temperature estimates and associated uncertainties. Enamel δ ¹⁸O to δ ¹⁸O_dw regressions used in this study were established using ordinary least square regression (OLS) applied to published data from bison (Hoppe, 2006) and cattle (D'Angela and Longinelli, 1990). Modern calibration data to establish the relationship between precipitation δ ¹⁸O and air temperature were obtained from Global Network of Isotopes in Precipitation (GNIP) stations (IAEA/WMO, 2020) in western, southern, and central Europe, as well as from circum-Mediterranean countries outside Europe that experience precipitation year-round. Usually, only stations with 5 or more years of measurements were used. As the slope of the δ ¹⁸O_precip–air temperature relationship is known to vary seasonally, separate calibration data sets were generated for reconstructing summer, winter, and mean annual temperatures, using GNIP data for the warmest month (commonly July or August; n_stations = 75), coldest month (commonly December or January; n_stations = 83), and long-term mean annual averages (n_stations = 89), respectively. Regression lines for the δ ¹⁸O_precip–air temperature relationships were also obtained using OLS. In this workflow, palaeotemperature estimates are made as an average for each stratigraphic unit, with data from all specimens from the same layer being combined into one temperature estimate. This enables the use of a sample size–dependent calculation of estimation uncertainty but introduces a degree of averaging of between-year climatic variability and biological or behavioural variability between individuals. For this reason, temperature seasonality derived from subtracting such pooled summer and winter temperature estimates may deviate somewhat from the full extent of annual temperature change that occurred in some years.

Predicting δ ¹⁸O_phos from δ ¹⁸O_carb

As we use the δ ¹⁸O_carb to δ ¹⁸O_phos relationship in this study only to convert δ ¹⁸O_carb measurements into δ ¹⁸O_phos values specifically for a group of study specimens with some paired data available, we chose to use the δ ¹⁸O_carb to δ ¹⁸O_phos relationship that is specifically derived from our study specimens to predict δ ¹⁸O_phos from samples that only have δ ¹⁸O_carb measurements, rather than using the relationship established from modern comparative data.

To test prediction performance of this approach, we conducted a number of performance tests. The R code for these tests can be found in the RMarkdown file of the Supplementary Material. The data set of Axlor δ ¹⁸O_phos/δ ¹⁸O_carb pairs (n = 47), was split into a training data set (two-thirds, N = 28) and a validation data set (one-third, N = 19), with random allocation of data points into each data set. A linear model (established using OLS) of the relationship between δ ¹⁸O_carb to δ ¹⁸O_phos in the training data set was used to predict δ ¹⁸O_phos for the validation data set, and the results were compared with measured δ ¹⁸O_phos.

Beyond pure δ ¹⁸O_phos prediction uncertainty, we also evaluate the effects of this uncertainty on seasonal palaeotemperature estimates, the outcome of interest in our study. For this purpose, we tested the impact of prediction uncertainty of δ ¹⁸O_phos from δ ¹⁸O_carb on the outcome of the inverse model, which serves as the input for seasonal palaeotemperature estimation. To simulate the noise potentially generated by the uncertainty in predicting δ ¹⁸O_phos, we chose a single δ ¹⁸O time series (individual AXL66) and randomly generated 10 alternative combined carbonate/phosphate time-series data sets by applying a random error onto the predicted δ ¹⁸O_phos values using random draws from a normal distribution with mean 0 and standard deviation of the root-mean-square error of the model cross-validation test, mimicking the distribution of residuals in δ ¹⁸O_phos prediction (Supplementary Fig. S3). We then processed each of these alternative time-series data sets through the inverse model (Supplementary Fig. S4) and compared the extracted summer peak and winter trough data points (Supplementary Fig. S5).

Software, code, and data

All stable isotope (δ ¹⁸O, δ ¹³C, δ ¹⁵N, δ ³⁴S) and radiocarbon data, as well as palaeotemperature estimates and any other data generated during this study are openly accessible as .csv files at https://github.com/ERC-Subsilience/Axlor_paleoclimatic_data. Also included in the repository are all files needed to run the Passey inverse model and to reproduce the palaeotemperature estimates. This article, including code for all data analyses, was written in R (v. 4.2.0; R Core Team, 2020) on a Windows 10 operating system, and the manuscript was rendered to a .docx file using RMarkdown (Allaire et al., 2018). The RMarkdown files used to produce both the main text and Supplementary Material can also be found in the associated online repository and can be used to reproduce the plots, tables, and statistical analyses of this article. The code for data analyses and production of the article manuscript and SI rendering make use of the following packages: matrixcalc_1.0-5 (Novomestky, 2021), Metrics_0.1.4 (Hamner and Frasco, 2018), ggpubr_0.4.0 (Kassambara, 2020a), patchwork_1.1.1 (Pedersen, 2020), magick_2.7.3 (Ooms, 2020), scales_1.2.0 (Wickham and Seidel, 2020), rstatix_0.7.0 (Kassambara, 2020b), officer_0.4.2 (Gohel, 2020a), stringr_1.4.0 (Wickham, 2019), flextable_0.7.0 (Gohel, 2020b), knitr_1.39 (Xie, 2014), wesanderson_0.3.6 (Ram and Wickham, 2018), captioner_2.2.3 (Alathea, 2015), cowplot_1.1.1 (Wilke, 2019), tidyr_1.2.0 (Wickham and Henry, 2019), dplyr_1.0.9 (Wickham et al., 2019), and ggplot2_3.3.6 (Wickham, 2016).

RESULTS

Radiocarbon dating

A new ¹⁴C date obtained for Layer III (OxA-40151) provided a radiocarbon age of 46,200 ± 3,000 yr BP, which calibrates between 45,510 cal yr BP and the end of the calibration curve at >55,000 cal yr BP (calibrated using IntCal20 [Reimer et al., 2020], OxCal version 4.4 [Bronk Ramsey, 2009]; see Supplementary Fig. S1). The date gives a median calibrated age of 49,690 cal yr BP, but it should be noted that due to the truncation of the posterior probability curve at the end of the calibration curve, this may represent a slight age underestimation. The collagen extract falls within acceptable levels of quality-control criteria, with a C/N ratio of 3.3, 32.4% C and 11.5% N, and stable carbon and nitrogen isotope delta values of −20.0‰ δ ¹³C and 6.4‰ δ ¹⁵N are consistent with expected values for a large bovine in this time period, indicating that the collagen is well preserved in this sample (Supplementary Table S1).

Carbonate and phosphate oxygen

As expected from empirical and theoretical data on the relationship between δ ¹⁸O_carb and δ ¹⁸O_phos, we observe a clear linear relationship between the two fractions in our data with a slope close to 1, matching previously published data on modern—and therefore pristinely preserved—bioapatites from a variety of mammal species (Fig. 2). With the exception of a single sample (AXL65.4), all Axlor data points fall within the predictive interval of the modern comparative data (Bryant et al., 1996b; Iacumin et al., 1996; Shahack-Gross et al., 1999; Zazzo et al., 2004; Miller et al., 2019), and even this outlier sample still falls within the variation observed in these modern data (Fig. 2). The regression of our data results in an equation with a slope slightly higher than 1:

$$\delta ^{18}O_{\,phos} = 1.1\cdot \delta ^{18}O_{carb}-11.8$$

while the regression equation of the modern comparative data shows a slope slightly below 1, with:

$$\delta ^{18}O_{\,phos} = 0.9\cdot \delta ^{18}O_{carb}-7.2$$

In a test of δ ¹⁸O_phos prediction performance using a linear model of the relationship between δ ¹⁸O_carb and δ ¹⁸O_phos from a randomly selected subset training data set of δ ¹⁸O_carb and δ ¹⁸O_phos measurement pairs obtained in this study (see section “Predicting δ ¹⁸O_phos from δ ¹⁸O_carb”), measured and predicted δ ¹⁸O_phos values fall close to a 1:1 relationship with minimal scatter (R² = 0.87; Supplementary Fig. S2). The prediction shows a root-mean-square error of 0.56‰, approximately twice the measurement uncertainty of δ ¹⁸O_phos (Supplementary Fig. S3). While this error is substantially larger than measurement uncertainty, it has negligible influence in the context of obtaining seasonal palaeotemperature information. Using the previously determined δ ¹⁸O_~phos prediction error (Supplementary Fig. S3), we generated random variants of one existing δ ¹⁸O time series (AXL66; Supplementary Fig. S4) and evaluated the variability between inverse model outcomes for these within-error prediction alternatives (see “Predicting δ ¹⁸O_phos from δ ¹⁸O_carb”); Supplementary Fig. S5). The root-mean-square error of extracted seasonal δ ¹⁸O data points from the alternative generated data sets compared with the original AXL66 time series is 1.21. This is much smaller than the mean inverse model uncertainty of ~4.5 (95% confidence interval).

Oxygen isotope palaeoclimatology

Summer peak, winter trough, and mean δ ¹⁸O_phos were extracted from sinusoidal curves to yield information about changes in summer, winter, and mean annual temperatures during the use of Axlor by late Neanderthals. Across Layers III, IV, and VI, δ ¹⁸O values are similar for summer, winter, and mean values (Fig. 3) and can be seen in detail in Table 4.

Table 4. Overview of summer peak, winter trough, and annual midpoint δ ¹⁸O_phos values for each specimen, as well as the means and standard deviations summarised for each layer.^a

^a SDs of individual specimens represent the measurement uncertainty, while SDs of the layer summaries represent the SD of the mean.

The single individual analysed from Layer IV falls on the higher end of the range of summer δ ¹⁸O values observed in the other layers, but within the variability observed in Layer III. A lack of detectable isotopic difference (given the available sample size) between analysed layers is supported by a lack of statistically significant (0.05 alpha level) differences between layers in summer, winter, or mean δ ¹⁸O (P _summer = 0.44, P _winter = 0.98, P _{mean annual} = 0.75).

Given the limited amount of published δ ¹⁸O data from Bos/Bison sp. that can serve as comparisons for our data, a comparison with other species is useful. To facilitate this, a conversion of enamel δ ¹⁸O values to drinking water δ ¹⁸O values is necessary, using a species-specific relationship. At the same time, this facilitates a comparison with modern measurements of δ ¹⁸O in precipitation and environmental waters. The δ ¹⁸O_dw for the Middle Palaeolithic Bos/Bison sp. from Axlor falls between ~−4‰ and −2‰ in summer and at ~−14‰ in winter, with mean δ ¹⁸O_dw between ~ −8‰ and −7‰ (Supplementary Fig. S6, Supplementary Table S2). To capture the pronounced orographic and isotopic variability in the vicinity of Axlor, modern comparison δ ¹⁸O_precip data were estimated for the site location as well the as the locations with lowest and highest elevations in a 50 km radius (Supplementary Fig. S6). For summer and annual means, the δ ¹⁸O_dw estimates from the Axlor data fall within the local variability of modern-day δ ¹⁸O_precip, while winter δ ¹⁸O_dw estimates are lower than modern δ ¹⁸O_precip values (Supplementary Fig. S6).

Results of seasonal and annual temperature estimations can be found in Supplementary Table S2. Mean annual palaeotemperature estimations for Layers III, IV, and VI fall close to modern-day conditions in the area of Axlor, with mean annual temperature estimates of ~12°C and estimation uncertainties between 2.4°C and 4.6°C depending on the number of teeth analysed per layer (Fig. 4). Modern comparisons made here were founded on a spatial grid of temperature data based on observations between 1981 and 2010 provided by the Spanish State Meteorological Agency (AEMET, 2021). It should be noted that spatial variability in microclimate is very pronounced in the Cantabrian and Bizkaia region due to pronounced altitudinal gradients in close proximity to the coast. For this reason, we include as comparisons both temperature data for the site of Axlor and any areas within 50 km of the site. In contrast to mean annual temperatures, summer and winter palaeotemperature estimates fall outside the variability observed within 50 km of the site, with summer temperatures being higher and winter temperatures lower than modern-day conditions (Fig. 4). This results in a stronger temperature seasonality estimated for the Middle Palaeolithic surroundings of Axlor compared with modern-day climates in the vicinity of the site. Temperature estimates are around ~26°C in summer and −5°C in winter, with uncertainties of ~3–5°C depending on the number of teeth analysed per layer (Supplementary Table S2). Temperature seasonality appears slightly more pronounced in Layer IV, but only a single individual was analysed, so a meaningful difference compared with the other layers cannot be established. Overall, temperature results are stable across the three layers that were analysed, mirroring raw δ ¹⁸O results.

Figure 4. Similar to raw oxygen isotope trends, reconstructed palaeotemperatures for summer, winter, and mean annual temperatures show little difference across the Middle Palaeolithic sequence of Axlor. Temperatures are close to and overlap with modern-day temperatures for summer and mean annual, but not winter, in the study region (based on data from the site location [dark bands] and a 50 km radius [light bands] obtained from a spatial grid of temperature data based on observations between 1981 and 2010 provided by the Spanish State Meteorological Agency [AEMET, 2021]).

Carbon, nitrogen, and sulphur stable isotopes

The carbon, nitrogen, and sulphur isotopic composition of bone collagen was analysed for red deer (Cervus elaphus), large bovines (Bos primigenius/Bison priscus), and to a lesser extent, horses (Equus ferus). An overview of results by taxon and layer can be found in Figure 5. Collagen was well preserved across specimens (n = 56 of N_total = 66 with collagen yield >1%). Collagen extractions used for stable isotope analysis had C/N ratios between 3.0 and 3.4, consistent with in vivo collagen (Van Klinken, 1999). C/S and N/S ratios were 640 ± 110 and 199 ± 34, respectively (mean and 1 SD), which is also consistent with good collagen preservation (Nehlich and Richards, 2009).

Figure 5. Carbon, nitrogen, and sulphur stable isotope values show little diachronic change for large bovines (Bos/Bison sp.), red deer (Cervus elaphus), or horses (Equus ferus). Intra-layer variability is generally high for all taxa and layers, but this is consistent across the sequence. Red deer are consistently lower in δ ¹⁵N compared with the other taxa, but otherwise different taxa show consistent isotopic similarity. The lack of diachronic change indicates a similarity of environments and plant and herbivore ecosystems between the different Middle Palaeolithic layers analysed here. Data points represent individual analysed bone fragments, connecting lines depict the group mean, while shaded violin plots show the shape of data distribution. Depictions of group means were removed for horses, due to the very low sample size and consequent lack of robustness for any apparent diachronic trends.

Results show little isotopic change in any of the isotopic systems and taxa across different layers at Axlor. In red deer, a slightly lower δ ¹³C can be seen in Layer VIII compared with the other layers, but this is very strongly driven by a single outlier. If this outlier is excluded, analysis of variance does not detect any statistically significant differences between layers for either C. elaphus or Bos/Bison sp. in any of the isotopic systems ($P_{{\rm B/B, }\delta ^{13}{\rm C}}$ = 0.96, $P_{{\rm B/B, }\delta ^{15}{\rm N}}$ = 1, $P_{{\rm B/B, }\delta ^{34}{\rm S}}$ = 0.7, P $_{{\rm deer, }\delta ^{13}{\rm C}}$ = 0.08, P $_{{\rm deer, }\delta ^{15}{\rm N}}$ = 0.58, P $_{{\rm deer, }\delta ^{34}{\rm S}}$ = 0.32). The mean isotopic delta values for horses appear to change between layers; however, this could be a product of a small sample size, and more data would be needed to draw any conclusions about isotopic change. Due to the low number of individuals analysed for horses, no statistical tests comparing the layers were conducted.

It is interesting to note that Bos/Bison sp. and C. elaphus are comparable to each other in their carbon and sulphur stable isotope values, with means of ~ −20‰ δ ¹³C and ~11‰ δ ³⁴S, respectively, but are noticeably distinct from each other in δ ¹⁵N, with means of 6.2‰ δ ¹⁵N in Bos/Bison sp. and 3.3‰ δ ¹⁵N in C. elaphus, a separation that is consistent across all analysed layers (Figs. 5–7). Intra-individual variability is pronounced for all taxa and isotopic systems, with values ranging over a ~2‰ spread in δ ¹³C, ~5‰ in δ ¹⁵N, and ~6‰ in δ ³⁴S (excluding a C. elaphus outlier). Horses also show such large variability between individuals despite the small number of samples analysed here. Scatter plots of the different isotopic systems (Figs. 6 and 7) do not reveal any correlations between carbon, nitrogen, or sulphur stable isotope delta values of large bovines or red deer.

Figure 6. A plot of carbon and nitrogen stable isotope values of Bos/Bison sp. and Cervus elaphus shows that the two herbivore taxa occupy distinct nitrogen isotopic spaces, while being very similar in δ ¹³C. This suggests a different feeding ecology or habitat use between these species, which remains stable across the layers analysed here. Data points represent measurements from individual bone fragments, while ellipses depict 1 SD (68%) of the data spread.

Figure 7. A plot of carbon and sulphur stable isotope values of Bos/Bison sp. and Cervus elaphus shows that these taxa heavily overlap in these isotopic systems. While the baseline variability in sulphur isotope delta values in Bizkaia is not well understood, a similarity in this isotopic system suggests that Bos/Bison sp. and Cervus elaphus did not systematically occupy spatially distinct habitats, particularly with regards to distance from the coast. Data points represent measurements from individual bone fragments, while ellipses depict 1 SD (68%) of the data spread.

DISCUSSION

Using combined carbonate-phosphate δ¹⁸O time series for palaeclimatology

The isotopic offset in δ ¹⁸O between the carbonate and the phosphate moiety (Δ¹⁸O_carb-phos) observed in Bos/Bison sp. tooth enamel samples in this study of 8.7 ± 0.7‰ (1 SD) closely matches the average Δ¹⁸O_carb-phos offset in published modern comparative data of 8.7 ± 0.9‰ (mean and 1 SD) (Bryant et al., 1996b; Pellegrini et al., 2011; Trayler and Kohn, 2016) as well as the theoretically predicted offset of ~8‰ (Aufort et al., 2017). The regression equation in the Axlor sample between δ ¹⁸O_carb and δ ¹⁸O_phos differs slightly from published equations and the cross-study sample equation derived here from modern comparative data. Such small discrepancies are expected and commonly occur between different studies, due to interindividual or interpopulation differences and limited coverage of isotopic space in different studies (Longinelli and Nuti, 1973; Bryant et al., 1996b; Iacumin et al., 1996; Martin et al., 2008; Pellegrini et al., 2011; Chenery et al., 2012; Trayler and Kohn, 2016). Overall, the good match between the Axlor samples and the theoretically and empirically documented relationship between δ ¹⁸O_carb and δ ¹⁸O_phos indicates that the enamel samples are well preserved and that it may be possible to predict δ ¹⁸O_phos from δ ¹⁸O_carb with reasonable accuracy and precision to reduce the amount of costly δ ¹⁸O_phos measurements needed.

The consistency of δ ¹⁸O_carb-phos between samples, individuals, species, or sites has been discussed at length in the literature, but predominantly in the context of diagenetic quality control (see section “Oxygen Isotope Analyses as a Proxy of Seasonal Climate”). Particularly in this context, several studies caution against the use of δ ¹⁸O_carb-phos as a fixed numeric indicator of preservation given the level of variability even between samples of the same individual (e.g., Pellegrini et al., 2011). However, an acceptable level of noise in this relationship differs depending on how data are used.

In this study, we investigated specifically the use case of predicting δ ¹⁸O_phos from δ ¹⁸O_carb for a subset of samples in a sequentially sampled δ ¹⁸O time series, which serves as the input for seasonal inverse modelling and subsequent palaeotemperature estimation of summer and winter temperatures (see sections “Oxygen Isotope Analyses as a Proxy of Seasonal Climate” and “Seasonal Palaeotemperature Estimation”). Using cross-validation in samples with both δ ¹⁸O_phos and δ ¹⁸O_carb measurements, it has been possible to determine that prediction uncertainty of δ ¹⁸O_phos is nonnegligible for raw (unmodelled) stable isotope data, approximately twice the measurement uncertainty. However, this uncertainty becomes much less significant when looking at the outcome of the inverse model. When reconstructing summer and winter palaeotemperatures, the summer peak and winter trough δ ¹⁸O values are extracted after inverse modelling, and only these peak and trough areas are ultimately used as input for palaeotemperature estimation (see section “Seasonal Palaeotemperature Estimation”). The intermediate areas of the sinusoidal δ ¹⁸O curve do influence the outcome of the inverse modelling procedure but are not themselves immediate areas of interest. In this approach, actual δ ¹⁸O_phos measurements were conducted on samples in the summer peak and winter trough areas. This apparently sufficiently serves to anchor these parts of the δ ¹⁸O curve and protects areas of interest for palaeotemperature estimation from being substantially influenced by the uncertainty of the δ ¹⁸O_phos to δ ¹⁸O_carb prediction. We simulated the effect of prediction uncertainty on δ ¹⁸O curve shape and inverse model outcomes by generating alternative δ ¹⁸O time series that represent possible variations of the δ ¹⁸O curve within the prediction uncertainty of δ ¹⁸O_phos from δ ¹⁸O_carb and compared the extracted seasonal δ ¹⁸O values of the inverse model outcome between each trial. Very little difference is noted between the summer and winter δ ¹⁸O values extracted from each trial, with variability less than half of the inverse model 95% confidence interval. Therefore, it can be concluded that the impact on summer and winter temperatures estimated from δ ¹⁸O curves, where intermediate (non-peak) curve areas have been predicted from δ ¹⁸O_carb, is negligible. This validates the use of a combination of the two measurement types that retains both robustness of palaeotemperature estimates and a enables a comparatively inexpensive workflow.

However, this conclusion applies specifically to δ ¹⁸O time-series data that include actual δ ¹⁸O_phos measurements in the peak and trough areas—the relevant sections for seasonal palaeoclimatic reconstructions. As would be expected, it can be seen in the inverse models of alternative simulated δ ¹⁸O curves (Supplementary Figs. S3 and S4) that predicted (non-measured) curve areas show significantly more variability between trials than the curve areas that are “anchored” by actual measurements. Consequently, one should be cautious when using completely predicted δ ¹⁸O_phos curves, and we recommend that any δ ¹⁸O_phos measurements should be taken from peak and trough areas. In this study, it is demonstrated that three δ ¹⁸O_phos measurements around each peak or trough area are necessary to sufficiently anchor these curve areas against excessive influence of the prediction uncertainty. Finally, systematic lagging sometimes proposed for δ ¹⁸O_carb vs δ ¹⁸O_phos time series has not been consistently observed (Pellegrini et al., 2011; Trayler and Kohn, 2016), but small variations in the location of the most extreme δ ¹⁸O values commonly exist between time series of the two fractions. Thus, using approximately three values rather than only individual samples for each peak or trough is recommended.

It should be noted that these conclusions are only valid in the case of inverse modelled sequential δ ¹⁸O series, but not necessarily for palaeotemperature estimates from unmodelled time series or bulk samples, as our arguments rest on the preservation of the δ ¹⁸O seasonal curve shape through the inverse modelling procedure due to the use of strategic δ ¹⁸O_phos anchor points. In case studies where δ ¹⁸O_phos values are completely predicted and/or unmodelled, it is our opinion that a substantially larger number of individuals would need to be analysed per archaeological layer to offset the level of noise introduced by the δ ¹⁸O_carb to δ ¹⁸O_phos prediction uncertainty (for further discussion of bulk samples, see Pryor et al., 2014; Skrzypek et al., 2016).

Reconstructing the environmental and ecological setting of the Middle Palaeolithic at Axlor

Chronology

Layers III to VIII of Axlor have previously been suggested to have formed during early MIS 3 and earlier, but the Beta dates generated from the site were produced without ultrafiltration removal of possible contaminants—the standard method currently applied in Middle Palaeolithic sites with a chronology close to the limit of the radiocarbon—and lack provenance information within the site (Ríos-Garaizar, 2017; Gómez-Olivencia et al., 2018; Marín-Arroyo et al., 2018; González-Urquijo et al., 2021; Sánchez Hernández, 2021; Table 1). An additional two infinite radiocarbon ages from Layer IV, produced by Marín-Arroyo et al. (2018), already indicated that the sequence may be older than previously suggested. This is supported again in this study by new data produced on ultrafiltered collagen from Layer III, which falls into early MIS 3. Of the dates that have been obtained from Layer D (equivalent to Layer IV) (e.g., Beta-203,107 and Beta-225,486; González-Urquijo et al., 2021) the former overlaps with the new age for Layer III presented here, while the latter provides an infinite age of >43,000 radiocarbon years. However, due to the concerns with the Beta dates mentioned earlier, we consider the ORAU dates to be more robust. With the upper Middle Palaeolithic sequence falling into early MIS 3 according to the newly presented radiocarbon date, an age of at least very early MIS 3 but possibly MIS 4 or older for the underlying layers is more likely. New data might confirm an older attribution for these lower levels VII and VIII corresponding to MIS 5c (Sánchez Hernández, 2021). However, additional studies of the site chronology, for instance using luminescence or electron spin resonance dating, are necessary to clarify the age of the lower part of the Middle Palaeolithic sequence of this important site.

Palaeoclimate

While temperature is the predominant driver of precipitation δ ¹⁸O values in northern Spain, other secondary impacts on δ ¹⁸O_precip, and therefore δ ¹⁸O_enamel, such as aridity, rainfall amount, differences in drinking water sources, changes in atmospheric circulation, or spatial variability from altitudinal gradients, exist as well. As discussed previously (see section “Oxygen Isotope Analyses as a Proxy of Seasonal Climate” and Supplementary Text 1), we argue that circulation changes likely play a minor role in this case study, but other factors should be examined before making climatic inferences.

In C₃ ecosystems, δ ¹⁸O and δ ¹³C can both be impacted by aridity and rainfall amount, a positive correlation between the two tracers is normally observed in settings where aridity and rainfall amount are the main drivers of both isotopic systems (Drucker et al., 2008; Diefendorf et al., 2010; Feranec et al., 2010; Richards et al., 2017). In the Axlor Bos/Bison sp. samples, however, neither annual means nor sequential data of δ ¹⁸O and δ ¹³C systematically correlate with each other (Supplementary Figs. S7 and S8), indicating that aridity is not a major driver of δ ¹⁸O and δ ¹³C in these individuals. Palaeoclimatic reconstruction could also be biased by Bos/Bison sp. accessing water from sources that are not isotopically representative of local precipitation. Pronounced isotopic deviation from seasonally variable local precipitation occurs most commonly in waterbodies such as large lakes, large rivers, snowmelt, or ground water, as these undergo substantial evaporation, are significantly time averaged, introduce a time lag, or transport water over large distances (Gonfiantini, 1986; Darling et al., 2003; Gat, 2010; Halder et al., 2015; Pederzani and Britton, 2019). Rivers in the study area are part of the watershed of the Ibaizabal River and are typically relatively small and short, in the range of tens of kilometres (Ocio et al., 2016). River water in the study region is, therefore, unlikely to be significantly time averaged, but rather is mostly fed from relatively local precipitation. While two lakes exist to the south of the site, the regular use of drinking water from substantially seasonally buffered water sources, such as larger lakes or ground water, is unlikely, given the seasonally variable δ ¹⁸O time-series data in all sampled teeth from Axlor. Due to the presence of higher-altitude areas surrounding Axlor, river water in the watershed likely includes some input of precipitation from a higher elevation, which is characterised by lower δ ¹⁸O values. For this reason, we account for this in our comparisons with modern water isotope data by including δ ¹⁸O_precip estimates from high-elevation locations to represent as much as possible the local variability in δ ¹⁸O of environmental waterbodies (Supplementary Fig. S6). Overall, we conclude that Bos/Bison sp. δ ¹⁸O values at Axlor are a useful proxy of palaeotemperatures, but it is necessary to take into account that local waters imbibed by animals ranging in the site vicinity likely included some higher-elevation precipitation. It should be noted that there was also a presence of small cirque glaciers in the eastern Cantabrian Mountains during the Quaternary (Serrano et al., 2017), which could have supplied low δ ¹⁸O glacial melt into the local watersheds. Potential low δ ¹⁸O inputs, therefore, need to be considered when making comparisons with modern δ ¹⁸O values of precipitation rather than river samples, as is the case in this study.

Oxygen stable isotope values of Bos/Bison sp. tooth enamel from Axlor are similar in all analysed layers (III, IV, and VI), indicating a comparability in temperature conditions between the different deposits (Fig. 3). The number of analysed individuals is relatively small, particularly in Layer IV. However, based on the available data, a consistency of climatic conditions between the different occupations is indicated. At the same time, both conversions to palaeotemperatures and comparisons with other published oxygen isotope data indicate that mean annual climatic conditions at Axlor were characterised by relatively high temperatures, akin to modern-day climates during the Middle Palaeolithic phases that are represented in the archaeological sequence, while temperature seasonality appears to have been more pronounced but stable across the different occupations.

The results from the temperature conversion resemble modern-day temperatures close to the site for mean annual conditions, while summer and winter palaeotemperature estimates fall higher and lower than modern-day, respectively (Fig. 4). This creates a more pronounced temperature seasonality for the Pleistocene sample represented here. It should be noted that summer δ ¹⁸O_dw estimates overlap within error with current summer δ ¹⁸O_precip values (Supplementary Fig. S6). This may indicate that the δ ¹⁸O_precip–T _air relationship in the Cantabrian region today has a lower slope than in the broader European data set used for palaeotemperature estimation in this study. Indeed δ ¹⁸O_precip data from Santander indicate a relatively low δ ¹⁸O_precip–T _air slope of ~0.3 (IAEA/WMO, 2020). Therefore, summer palaeotemperatures in our study may perhaps best be described as similar to modern day or higher, while winter values are lower than modern day for both δ ¹⁸O_dw and temperature estimates, indicating a stronger temperature seasonality in the Pleistocene sample. A possible contribution of glacial melt from Quaternary glaciers in the Cantabrian Mountains (Serrano et al., 2017) means that some palaeotemperatures may be underestimates of actual local temperatures. However, glacial melt and snowmelt would be expected to mostly impact spring and summer surface waters, when temperatures are higher and increase melting of snow and ice. An underestimation of winter temperatures from this effect is therefore not expected. We therefore conclude that most likely a higher temperature seasonality, mostly driven by lower winter temperatures, characterised the climate of the Middle Palaeolithic occupations of Axlor.

A higher temperature seasonality may have been driven by different orbital configuration, resulting in higher insolation during summer and lower insolation during winter. Several phases of orbital configurations producing stronger seasonality of insolation, and thus temperature, fall into the late Pleistocene, for instance, at ~60 ka BP and ~85 ka BP (Laskar et al., 2004). At the same time, the level of the Cantabrian Sea was lower than today during MIS 3 and MIS 4 as well as most of MIS 5 (Ballesteros et al., 2017), meaning that Axlor Cave was situated further inland at the time, which could have contributed to a more continental climate. In current climate conditions, the area around Axlor exhibits pronounced differences in microclimate over short distances, with the coastal areas being substantially wetter and with less pronounced seasonality, while southern valleys experience less precipitation and a temperature seasonality that is up to 6°C higher than at the coast (Pellitero et al., 2019). At the same time, these microclimatic differences over relatively short distances may mean that palaeoclimatic reconstructions could also relate to prey procurement in slightly different areas and microclimates than were present directly at the site. While we made efforts to capture the local isotopic variability in our comparisons of archaeological and modern data, a use of more southern valleys or lower-elevation areas by hunted Bos/Bison sp. cannot be fully excluded. Indeed, it has been previously proposed that Bos/Bison sp. game may have been procured at some distance from the site, due to the steep terrain around the cave (Ríos-Garaizar and Moreno, 2015).

For mean annual conditions, a similarity to Holocene climatic conditions is also supported by more direct comparisons of oxygen stable isotope data. Mean annual δ ¹⁸O_dw values also match those reconstructed for residents of a twelfth-century CE site located ~15 km from Axlor in a neighbouring valley (Guede et al., 2020). Estimates of δ ¹⁸O_dw for human tooth enamel from this site range from ~ −5‰ to −10‰. Strontium stable isotope analyses of the medieval human remains showed use of various areas in the vicinity of the site by the site inhabitants, making their oxygen stable isotope values more broadly representative of the isotopic variability of drinking water sources in the site area. A match of the Axlor data with these δ ¹⁸O_dw values further underscores a climatic similarity with the relatively mild Holocene climatic conditions of the twelfth-century CE.

Comparisons with other herbivore oxygen stable isotope data also show that Axlor Bos/Bison sp. oxygen isotope values match well with values for warmer phases of the last glacial period. The δ ¹⁸O values from Axlor are similar to those generated from Bos/Bison sp. remains from the MIS 3 and MIS 4 layers of La Ferrassie in southwest France (mean δ ¹⁸O_phos = 16.3 ± 0.7‰, mean summer δ ¹⁸O_phos = 17.9 ± 0.8‰, mean winter δ ¹⁸O_phos = 14.2 ± 1.0‰; Pederzani et al., 2021). Our reconstructed δ ¹⁸O_dw values are also slightly higher, but within error of δ ¹⁸O_dw values calculated from δ ¹⁸O_carb values of horse tooth enamel from the Mousterian layer 20E (undated MIS 3) of El Castillo (calculated in this study as −9.8 ± 2.8‰), a Cantabrian site located ~150 km west of Axlor with a similar distance to the coast (Jones et al., 2018). We used only the individual CAS61 as a comparison, as the other individuals from El Castillo do not show clearly sinusoidal δ ¹⁸O curves. Axlor δ ¹⁸O mean annual values are ~3‰ lower than a Bos/Bison sp. sample from the approximately contemporary site of Valdegoba Cave near Burgos (δ ¹⁸O = −8.6 ± 0.9‰) (Feranec et al., 2010). This difference closely matches the difference between the two localities in modern-day precipitation δ ¹⁸O (~−6‰ in Axlor [Bowen and Revenaugh, 2003] and ~ −9‰ in Burgos [IAEA/WMO, 2020]), showing again that our results conform with the expectations for a warm-phase Middle Pleniglacial climate. The fact that this mainly hydrotopographically driven spatial difference in δ ¹⁸O is well reflected in enamel samples also further supports the effectiveness of using δ ¹⁸O_enamel as a proxy for water δ ¹⁸O. The observed lack of correlation between δ ¹⁸O and δ ¹³C in the Bos/Bison sp. enamel carbonates (Supplementary Fig. S8) or δ ¹³C and δ ¹⁵N in bone collagen in Bos/Bison sp. or Cervus elaphus (Fig. 6) additionally suggests that relatively high δ ¹⁸O are indeed more likely related to warm mean annual temperature conditions rather than an aridity effect.

The similarity in δ ¹⁸O data between Axlor and La Ferrassie indicates that comparable conditions with warm mean annual temperatures characterised both the southwest of France and the northern Spanish coast during the Middle Pleniglacial. Under current climatic conditions in the southwest of France, mean annual temperatures are also similar to those in the Basque Country. However, the temperature seasonality is currently higher in France compared with the weather station closest to Axlor. The similarities between the Axlor and La Ferrassie data indicate that the archaeological deposits at Axlor were formed during more seasonal climatic conditions than today, creating a reduced climate gradient with the southwest of France. Whether Neanderthals occupied Axlor during the winter season is unknown, but seasonality studies in the Cantabrian region suggest that a number of sites were occupied in various seasons or year-round, while climatic and herbivore ecology were stable throughout the Middle Palaeolithic (e.g., Sánchez-Hernández et al., 2019). Lower winter temperatures likely influenced the growing season and availability of particular plant resources especially from forested habitats, but studies of dental calculus have shown that Neanderthals living in more open, colder environments incorporated a variety of plant taxa into their diet, including grass seeds, with similar dietary breadth as individuals living in warmer, more closed environments (Power et al., 2018). At the same time, climatic conditions at Axlor appear to have supported a rich and stable herbivore resource base, which likely formed an important part of Neanderthal subsistence (see section “Palaeoenvironment and Palaeoecology”).

It should be noted that the faunal assemblage from the Barandiaran excavation, which was used here, was recovered using an excavation methodology following artificial spits, likely leading to some imprecision in assigning faunal material to individual layers (Ríos-Garaizar, 2012; Gómez-Olivencia et al., 2018). While not ideal, we argue that the pattern of relatively homogeneous climatic conditions observed in the isotopic results is unlikely to exclusively be artefacts of the excavation methodology. If strong climatic differences had originally existed, but had been obscured by a mixing of the material between layers, a diachronically more homogeneous pattern but with a large intra-layer δ ¹⁸O variability would be expected. Because the intra-layer variability in δ ¹⁸O is not unusually high (e.g., see the comparison with data from La Ferrassie), a strong temporal mixing of climatically different phases seems unlikely. While the other isotopic tracers show larger intra-layer variability, this pattern is typical for other sites in the region (Jones et al., 2018, 2019; Marín-Arroyo et al., 2020).

Overall, our results indicate that the late Neanderthal occupations analysed here took place during mild mean annual climatic conditions that are similar between the different occupations analysed. Due a lack of clear chronological information for some archaeological layers at Axlor, we cannot precisely pinpoint the mechanism behind this climatic similarity, but several scenarios would be consistent with the isotopic data:

1. 1. The Bizkaia region was characterised by relatively mild climates throughout the time when Neanderthals occupied Axlor Cave, and climatic oscillations were less pronounced than in other parts of Europe.

2. 2. Late Pleistocene time periods of colder climate are not represented in the material accumulation of Layers III, IV, and VI, because the site was not used by Neanderthals when climatic conditions were harsher. Instead Layers III, IV, and VI represent repeated occupations during warmer phases/interstadials.

A buffering of some of the more extreme environmental fluctuations of the late Pleistocene has been suggested previously for the Cantabrian region (e.g., Jones et al., 2018). Long-term climatic records from the region do show evidence of millennial-scale climatic change, including episodes of colder, drier climates in MIS 5–3 (e.g., Sánchez Goñi et al., 2013). For example, a speleothem from Cueva de las Perlas, ~100 km west of Axlor, shows a more arid episode around 51–46.5 ka, expressed as growth rate changes (Deeprose, 2018). On the other hand, there are also examples of other terrestrial palaeoclimatic records that are not specifically tied to hominin activity, such as the micromammals from Cueva del Conde in Asturias showing a relative stability in climate in northern Spain during MIS 3 (e.g., López-García et al., 2011).

Due to the presence of sterile layers in the sequence of Layer III and below, an intermittent occupation of the site seems likely (Gómez-Olivencia et al., 2018), but additional chronological clarification and site-specific climatic information is needed to test whether this pattern is related to disuse of the site during unfavourable climatic conditions or not. Similarly, due to the lack of absolute chronological information for much of the Middle Palaeolithic sequence of Axlor, it is currently not possible to determine how much time is represented in these deposits. Nonetheless, some observations about the climatic framework for Neanderthal behaviour can be made within the context of the site itself. Layers III–V of Axlor have previously been proposed to have formed during less favourable climatic conditions, based on the presence of some reindeer (Rangifer tarandus) specimens, mainly in Layer IV, and the presence of Quina Mousterian lithic technology (Gómez-Olivencia et al., 2014; Ríos-Garaizar, 2017). However, it should be kept in mind that the faunal spectrum of these layers is generally dominated by red deer and aurochs/bison, which are erythermic species well adapated to both warm and cold conditions (Ríos-Garaizar and Moreno, 2015). In contrast, the number of reindeer specimens found at the site is very low (n = 21 for all deposits), and their generally sparse presence in the Cantabrian region suggests that the area fell towards the edge of the reindeer biogeographic range (Gómez-Olivencia et al., 2014). Additionally, given the climatic and ecological flexibility of larger mammals such as reindeer, we believe that the presence of low numbers of reindeer in Layers III–V are still compatible with the stable isotopic data from Bos/Bison sp., which suggest generally mild mean annual temperatures in all analysed deposits. It is also possible that the presence of reindeer was supported by lower winter temperatures, as reconstructed in our data set. However, whether reindeer migrated to the Cantabrian region in winter or were present year-round is currently unknown (Gómez-Olivencia et al., 2014), so this point warrants further investigation. On the other hand, it is also possible that accumulations of reindeer and larger bovines originate from climatically different short occupations that have been combined through the palimpsest character of the site, which is thought to represent repeated short and intensive occupations (Ríos-Garaizar, 2017). A more detailed investigation of the spatial distribution of reindeer bones within the site would help to make a stronger case one way or the other, but given that δ ¹⁸O variability between different Bos/Bison sp. individuals assigned to one stratigraphic unit at Axlor is comparable to that of other Middle Palaeolithic sites, we suggest that there is at present no indication of an unusually pronounced palimpsest with accumulation of Bos/Bison sp. individuals across multiple climatic phases into single stratigraphic units.

Palaeoenvironment and palaeoecology

The lack of diachronic changes seen in the carbon and nitrogen isotopic composition of Cervus elaphus and Bos/Bison sp. bone collagen (Fig. 5) indicates a stability of feeding ecology and of the plant ecosystem throughout the analysed layers. This matches with the climatic homogeneity indicated by stable δ ¹⁸O of Bos/Bison sp. tooth enamel discussed in the previous section, drawing overall a picture of climatic, environmental, and ecological homogeneity between the different Middle Palaeolithic occupations of the site. A stability in the geospatially variable sulphur stable isotope system (Fig. 5) might also indicate that animals continued to range over—and be hunted in—similar sulphur isozones through time. The range of sulphur isotope values represented within taxonomic and layer groups is comparatively large (Fig. 5) compared with the common 2–4‰ intra-individual variability observed in modern studies of herbivores from a single population or at other Palaeolithic archaeological sites (e.g., Nehlich, 2015; Jones et al., 2018; Reade et al., 2020). Thus, an exploitation of a relatively broad variety of isotopically distinct habitats by both animals and perhaps hominins throughout the Middle Palaeolithic occupations of Axlor seems likely. Baseline geospatial variability in δ ³⁴S is relatively poorly understood in the northern Spanish coastal region, but δ ³⁴S values are likely to at least vary with distance from the coast, where higher δ ³⁴S values prevail closer to the coast (McArdle et al., 1998; Nehlich, 2015; Bataille et al., 2021). Based on comparison with modelled δ ³⁴S values, coastal δ ³⁴S could be ~4–6‰ higher than in valleys further inland, and δ ³⁴S values of fauna at Axlor seem to not capture isozones that are directly adjacent to the coastline (Bataille et al., 2021). However, it should be noted that the isoscape by Bataille et al. (2021) is based on a relatively small number of samples from Holocene archaeological contexts, and baseline values are likely to have been different in the Pleistocene due to climatic and environmental differences. Therefore, more data are necessary to confirm this and to clarify the geographic scale that is associated with the relatively large range of different δ ³⁴S values seen in the fauna at Axlor. As the local and regional δ ³⁴S variability is largely unknown at the moment, this higher variability may also be related to a relatively high degree of baseline δ ³⁴S heterogeneity in a mosaic landscape characterised by deep narrow valleys that may shield valley floors from sea spray effects despite relative proximity to the coast.

At the same time, the range of isotopic values represented in the nitrogen isotopic system is also quite large, while carbon isotope values vary less strongly within each taxonomic group (Fig. 5). The pronounced within-taxon and within-layer variability in δ ¹⁵N may indicate that a range of different environments or portions of the plant ecosystem were exploited by the herbivore taxa analysed here. Such pronounced variability could be consistent with two scenarios:

1. 1. Each layer represents a relatively pronounced temporal palimpsest, and animals in each stratigraphic unit lived under a variety of environmental conditions, or taxa fulfilled a different ecological role in different time periods.

2. 2. Isotopic variability in herbivores reflects a pronounced spatial variability of stable isotope values in the plant ecosystem due to the presence of a mosaic landscape with heterogeneous vegetation structure.

Large intra-layer isotopic variability has been previously noted at other Cantabrian archaeological sites, but due to a lack of data from more directly climatically driven tracers such as oxygen stable isotope values, the cause for this could previously not be clearly determined (Jones et al., 2018, 2019; Marín-Arroyo et al., 2020). Here, the intra-layer variability in the oxygen isotopic composition of Bos/Bison sp. tooth enamel is not uncommonly large in any of the analysed layers, but is instead comparable with sites in other regions, such as La Ferrassie (Pederzani et al., 2021), that also do not exhibit such large intra-layer variability in other isotopic systems. Consequently, it appears unlikely that the temporal palimpsest of the Axlor archaeological deposits is solely responsible for the large range of values observed in nitrogen and sulphur stable isotope values. Animals should originate from a large variety of climatic contexts and show high δ ¹⁸O variability under such a scenario. Therefore, the combination of the multiple types of stable isotope data here lends more support to the presence of a relatively heterogeneous, isotopically variable, mosaic landscape that offered a breadth of vegetation types and habitats during the late Pleistocene.

At the same time, our data also shed light on the feeding ecology of different herbivore taxa. C. elaphus show lower δ ¹⁵N than Bos/Bison sp. in all analysed layers (Fig. 6), indicating that they consistently occupied a different isotopic niche. This could be caused by these taxa either fulfilling a different ecological role with different feeding behaviours but within the same or similar landscapes and habitats, or exploiting spatially distinct landscapes. Higher δ ¹⁵N values of Bos/Bison sp. could be linked to a higher proportion of grass in the diet, compared with the more browsing C. elaphus, as higher δ ¹⁵N values have been observed in grasses and forbs compared with shrubs or trees in arctic and peri-arctic environments (Drucker, 2022). A selective feeding in drier microenvironments may also be possible, but the differences between large bovines and red deer do not map onto any differences in sulphur stable isotopes (Fig. 7). Therefore, a spatial separation of the two taxa seems less likely. Thus, an explanation in relation to feeding ecology seems more likely, with the two species using (within isotopically detectable bounds) the same ranges but occupying separate isotopic feeding niches, likely along a grazer–browser spectrum. The consistency of how these herbivores relate to each other in their isotope ecology also supports the interpretation that Neanderthals occupying Axlor did so within a relatively consistent climate, environment, and ecosystem.

CONCLUSIONS

Oxygen, carbon, nitrogen, and sulphur stable isotope analysis of macromammal remains, with evidence of hominin manipulation, provides a snapshot of climates, environments, and ecosystems in the surroundings of Axlor Cave. Climatic and environmental conditions appear to have been remarkably stable across different punctuated episodes when Neanderthals inhabited the site during the late Pleistocene. Oxygen stable isotope values from large bovine tooth enamel indicate that mean annual temperatures were mild and similar to modern-day climates in the region, while winter temperatures were lower and temperature seasonality more pronounced (but stable) across occupations. A similarity in the environment across layers is also indicated by carbon and nitrogen isotope ratios in herbivore bone collagen. Herbivore taxa such as red deer and aurochs/bison appear to have occupied slightly different feeding niches, with large bovines likely consuming more grass, but their relationship to each other is consistent in all layers, indicating niche partitioning but also niche conservation through time. Likewise, noticeable changes in spatial ecology are not observed, as indicated by homogeneous sulphur stable isotope values in all analysed layers. Intra-layer ranges in sulphur and nitrogen stable isotope values are high, a phenomenon also observed at other sites in the region. Due to a lack of such large variability in the oxygen isotope sample, it is most plausible that large ranges in sulphur and nitrogen stable isotopes are driven by the presence of isotopically diverse habitats in mosaic environments within the complex topography of the site surroundings, rather than by the effects of the temporal palimpsest character of the archaeological deposits.

The somewhat unclear chronology of the site and a lack of naturally accumulated fauna that could serve as a comparison mean that it is challenging to clarify whether the pattern of climatic similarity between layers results from a buffered climate in the region during this time or from a punctuated settlement pattern in which sampled deposits formed during similar mild climatic phases. Further study of the site's chronology and the expression of climatic oscillations in terrestrial palaeoclimatic archives in the region are needed to clarify this. Such work would substantially increase the specificity of the palaeoclimatic data generated here. Nevertheless, the isotopic data allow some site-specific inferences about the archaeological record. For example, the data presented here may tentatively support a hypothesis that differences in the technological record of the site, between the lower and upper part of the Middle Palaeolithic sequence, are not related to strong climatic or environmental differences between these different occupation phases. Other factors, such as changes in mobility patterns, technological tradition, or prey species and hunting technique could all be considered as alternative explanations to be explored in future research. This case study adds to expanding research indicating that late Pleistocene Neanderthal groups in southwest Europe may have preferred exploiting stable and productive ecosystems with mild mean annual temperatures (e.g., Pederzani et al., 2021; Vidal-Cordasco et al., 2022). At the same time, the repeated occupations of Axlor and climatically similar Neanderthal occupations at sites in both northern Iberia and southwest France (discussed earlier) suggest that Neanderthal groups were well adapted to the conditions represented here, including the more pronounced temperature seasonality and colder winter temperatures, the effects of which may have been buffered by a stable ecosystem of herbivore prey. More broadly, our results illustrate the utility of approaches such as multi-isotope archaeozoology to generate the site-specific and hominin-associated data that are needed to characterise the ecological adaptations of late Pleistocene hominins.

In addition to making palaeoclimatic and palaeoenvironmental inferences, this study also explores methodological aspects of palaeotemperature estimation related to the use of oxygen stable isotope measurements of bioapatite carbonate and bioapatite phosphate. An approach that combines measurements of the two analytical fractions is presented. Statistical simulations show that uncertainty introduced by predicting δ ¹⁸O_phos from δ ¹⁸O_carb has a negligible effect on summer and winter palaeotemperature estimations, if primary δ ¹⁸O_phos measurements are used to anchor summer peak and winter trough areas of the δ ¹⁸O curve. Therefore, this study suggests such an approach is a viable option for case studies that are looking to achieve robust seasonal palaeotemperature estimation using a comparatively low-cost workflow.