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, Reference O'Regan2008; González-Sampériz et al., Reference González-Sampériz, Leroy, Carrión, Fernández, García-Antón, Gil-García, Uzquiano, Valero-Garcés and Figueiral2010; Jennings et al., Reference Jennings, Finlayson, Fa and Finlayson2011; Bradtmöller et al., Reference Bradtmöller, Pastoors, Weninger and Weniger2012; López-García et al., Reference López-García, Blain, Sanz and Daura2012; Sánchez Yustos and Diez Martín, Reference Sánchez Yustos and Diez Martín2015; Bicho and Carvalho, Reference Bicho and Carvalho2022). 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., Reference Maroto, Vaquero, Arrizabalaga, Baena, Baquedano, Jordá and Julià2012; Higham et al., Reference Higham, Douka, Wood, Ramsey, Brock, Basell and Camps2014; Straus and González Morales, Reference Straus and González Morales2016; Marín-Arroyo et al., Reference Marín-Arroyo, Rios-Garaizar, Straus, Jones, De la Rasilla, González Morales, Richards, Altuna, Mariezkurrena and Ocio2018) as well as a small number of assemblages that are consistent with the Châtelperronian technocomplex (Ríos-Garaizar et al., Reference Ríos-Garaizar, Iriarte, Arnold, Sánchez-Romero, Marín-Arroyo, Emeterio and Gómez-Olivencia2022). Extensive archaeological and archaeozoological assemblages have been excavated from many of these sites (Ríos-Garaizar et al., Reference Ríos-Garaizar, Arrizabalaga and Villaluenga2012; Yravedra-Sainz de los Terreros et al., Reference Yravedra-Sainz de los Terreros, Gómez-Castanedo, Aramendi-Picado, Montes-Barquín and Sanguino-González2016; Rasilla et al., Reference Rasilla, Duarte, Sanchis, Carrión, Cañaveras, Marín-Arroyo and Real2020; Álvarez-Vena et al., Reference Álvarez-Vena, Álvarez-Lao, Laplana, Quesada, Rojo, García-Sánchez and Menéndez2021), and key lines of research into late Neanderthal subsistence, behaviour, and resilience are ongoing in the region (Gómez-Olivencia et al., Reference Gómez-Olivencia, Sala, Núñez-Lahuerta, Sanchis, Arlegi and Rios-Garaizar2018; Marín-Arroyo et al., Reference Marín-Arroyo, Geiling, Jones, González Morales, Straus and Richards2020; Marín-Arroyo and Sanz-Royo, Reference Marín-Arroyo and Sanz-Royo2022; Vidal-Cordasco et al., Reference Vidal-Cordasco, Ocio, Hickler and Marín-Arroyo2022).
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., Reference Roucoux, Abreu, Shackleton and Tzedakis2005; Martrat et al., Reference Martrat, Grimalt, Shackleton, Abreu, Hutterli and Stocker2007; Sánchez Goñi et al., Reference Sánchez Goñi, Landais, Fletcher, Naughton, Desprat and Duprat2008; Moreno et al., Reference Moreno, Svensson, Brooks, Connor, Engels, Fletcher and Genty2014). 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., Reference Moreno, Svensson, Brooks, Connor, Engels, Fletcher and Genty2014; Fernández-García et al. Reference Fernández-García, Vidal-Cordasco, Jones and Marín-Arroyo2023).
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., Reference Wood, Arrizabalaga, Camps, Fallon, Iriarte-Chiapusso, Jones and Maroto2014, Reference Wood, Bernaldo de Quirós, Maíllo-Fernández, Tejero, Neira and Higham2018; Jones et al., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018, Reference Jones, Marín-Arroyo, Straus and Richards2020; Marín-Arroyo et al., Reference Marín-Arroyo, Rios-Garaizar, Straus, Jones, De la Rasilla, González Morales, Richards, Altuna, Mariezkurrena and Ocio2018). 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., Reference Drucker, Bocherens and Billiou2003, Reference Drucker, Kind and Stephan2011; Stevens and Hedges, Reference Stevens and Hedges2004; Stevens et al., Reference Stevens, Jacobi, Street, Germonpré, Conard, Münzel and Hedges2008, Reference Stevens, Hermoso-Buxán, Marín-Arroyo, González-Morales and Straus2014; Bocherens et al., Reference Bocherens, Drucker and Madelaine2014; Nehlich, Reference Nehlich2015; Jones et al., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018, Reference Jones, Richards, Reade, Bernaldo de Quirós and Marín-Arroyo2019, Reference Jones, Marín-Arroyo, Straus and Richards2020; Pederzani and Britton, Reference Pederzani and Britton2019; Reade et al., Reference Reade, Tripp, Charlton, Grimm, Leesch, Müller and Sayle2020). 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., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018, Reference Jones, Richards, Reade, Bernaldo de Quirós and Marín-Arroyo2019), 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., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018, Reference Jones, Marín-Arroyo, Straus and Richards2020; Rodríguez-Almagro et al., Reference Rodríguez-Almagro, Sala, Wißing, Arriolabengoa, Etxeberria, Rios-Garaizar and Gómez-Olivencia2021). Such palaeoecological information has recently increasingly been supplemented by information on spatial ecology from sulphur stable isotope analysis (Jones et al., Reference Jones, Richards, Reade, Bernaldo de Quirós and Marín-Arroyo2019; Reade et al., Reference Reade, Tripp, Charlton, Grimm, Leesch, Müller and Sayle2020). 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., Reference Jones, Richards, Reade, Bernaldo de Quirós and Marín-Arroyo2019).
Larger data sets from more strongly climatically driven stable isotope systems, such as δ 18O, 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., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018, Reference Jones, Richards, Reade, Bernaldo de Quirós and Marín-Arroyo2019; Marín-Arroyo et al., Reference Marín-Arroyo, Geiling, Jones, González Morales, Straus and Richards2020). 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 (δ 18O, δ 13C, δ 15N, δ 34S) 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 δ 18O 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 δ 18O 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, δ 13C, δ 15N, and δ 34S 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 (δ 13C, δ 15N), as well as herbivore spatial ecology (δ 34S) 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, Reference Ríos-Garaizar, Moreno, Conard and Delagnes2015). 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, Reference Basaguren and Orive1991). 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, Reference Basaguren and Orive1991). 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., Reference Ocio, Stocker, Eraso and Cowpertwait2016).
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 (Reference Barandiarán and Barandiarán1980; for further discussion, see Gómez-Olivencia et al. Reference Gómez-Olivencia, Sala, Núñez-Lahuerta, Sanchis, Arlegi and Rios-Garaizar2018; 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, Reference Ríos-Garaizar2017). The upper part of the Middle Palaeolithic sequence also contains a substantial number of bone tools (Mozota Holgueras, Reference Mozota Holgueras2012).
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., Reference González-Urquijo, Ríos Garaizar and Ibáñez Estévez2002), but has recently been included in Layer B (González-Urquijo et al., Reference González-Urquijo, Bailey and Lazuen2021).
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, Reference Altuna, Pathou and Freeman1989) and the collection from the recent re-excavations (Castaños, Reference Castaños, Montes Barquín and Lasheras Corruchaga2005), 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, Reference Altuna, Pathou and Freeman1989; Castaños, Reference Castaños, Montes Barquín and Lasheras Corruchaga2005; Gómez-Olivencia et al., Reference Gómez-Olivencia, Sala, Núñez-Lahuerta, Sanchis, Arlegi and Rios-Garaizar2018). 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., Reference Jones, Marín-Arroyo, Straus and Richards2020, Reference Jones, Marín-Arroyo, Corchón Rodríguez and Richards2021; but see Yravedra and Cobo-Sánchez, Reference Yravedra and Cobo-Sánchez2015). 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., Reference Gómez-Olivencia, Sala, Núñez-Lahuerta, Sanchis, Arlegi and Rios-Garaizar2018), but species abundances in the material here are very similar to those in the more recent excavations (Castaños, Reference Castaños, Montes Barquín and Lasheras Corruchaga2005) 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., Reference Marín-Arroyo, Rios-Garaizar, Straus, Jones, De la Rasilla, González Morales, Richards, Altuna, Mariezkurrena and Ocio2018; 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., Reference Gómez-Olivencia, López-Onaindia, Sala, Balzeau, Pantoja-Pérez, Arganda-Carreras, Arlegi, Rios-Garaizar and Gómez-Robles2020). 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, Reference Higham2011). Accordingly, our examination of chronological hygiene favours the Oxford Radiocarbon Accelerator Unit (ORAU) ultrafiltration bone 14C 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, Reference Ríos-Garaizar2017), 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, Reference Sánchez Hernández2021).
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, Reference Basabe1973), while others were only identified as hominin fossil remains during later reanalyses of the collections (Gómez-Olivencia et al., Reference Gómez-Olivencia, López-Onaindia, Sala, Balzeau, Pantoja-Pérez, Arganda-Carreras, Arlegi, Rios-Garaizar and Gómez-Robles2020). 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., Reference Gómez-Olivencia, López-Onaindia, Sala, Balzeau, Pantoja-Pérez, Arganda-Carreras, Arlegi, Rios-Garaizar and Gómez-Robles2020). Other hominin remains, described initially by Basabe (Reference Basabe1973), 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., Reference Gómez-Olivencia, López-Onaindia, Sala, Balzeau, Pantoja-Pérez, Arganda-Carreras, Arlegi, Rios-Garaizar and Gómez-Robles2020, Reference Gómez-Olivencia, López-Onaindia, Sala, Balzeau, Pantoja-Pérez, Arganda-Carreras, Arlegi, Rios-Garaizar and Gómez-Robles2023), while others support the original assignment to Neanderthals and argue against any influence of disturbance (González-Urquijo et al., Reference González-Urquijo, Bailey and Lazuen2021). A recent reply to the latter marks the end of this debate (Gómez-Olivencia et al., Reference Gómez-Olivencia, López-Onaindia, Sala, Balzeau, Pantoja-Pérez, Arganda-Carreras, Arlegi, Rios-Garaizar and Gómez-Robles2023), 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 (δ 18O) of faunal tooth enamel as a proxy of past climates rests on two underlying relationships between δ 18O of tooth enamel (δ 18Oenamel) and δ 18O of animal drinking water (δ 18Odw) and between δ 18O of environmental water (δ 18Oew) and air temperature, which is a dominant effect in mid- to high latitudes (see below in this paragraph) (Longinelli, Reference Longinelli1984; Luz et al., Reference Luz, Kolodny and Horowitz1984; Rozanski et al., Reference Rozanski, Araguás-Araguás and Gonfiantini1993; Iacumin et al., Reference Iacumin, Bocherens, Mariotti and Longinelli1996; Kohn et al., Reference Kohn, Schoeninger and Valley1996; Clark and Fritz, Reference Clark and Fritz1997; Müller et al., Reference Müller, Stumpp, Sørensen and Jessen2017). Tooth enamel is formed in isotopic equilibrium with animal body water (Longinelli, Reference Longinelli1984; Luz et al., Reference Luz, Kolodny and Horowitz1984; Iacumin et al., Reference Iacumin, Bocherens, Mariotti and Longinelli1996; Daux et al., Reference Daux, Lécuyer, Héran, Amiot, Simon, Fourel, Martineau, Lynnerup, Reychler and Escarguel2008), which in turn mirrors variability in δ 18O 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, Reference D'Angela and Longinelli1990; Bryant et al., Reference Bryant, Luz and Froelich1994, Reference Bryant, Froelich, Showers and Genna1996a; Sánchez Chillón et al., Reference Sánchez Chillón, Alberdi, Leone, Bonadonna, Stenni and Longinelli1994; Hoppe et al., Reference Hoppe, Amundson, Vavra, McClaran and Anderson2004; Hoppe, Reference Hoppe2006). Consequently, δ 18Oenamel of these animals empirically shows a strong linear correlation with δ 18O of drinking water (Longinelli, Reference Longinelli1984; Luz et al., Reference Luz, Kolodny and Horowitz1984; Rozanski et al., Reference Rozanski, Araguás-Araguás and Gonfiantini1992; Hoppe, Reference Hoppe2006; Tütken et al., Reference Tütken, Furrer and Walter Vennemann2007; Arppe and Karhu, Reference Arppe and Karhu2010; Pryor et al., Reference Pryor, Stevens, O'Connell and Lister2014; Skrzypek et al., Reference Skrzypek, Sadler and Wiśniewski2016). The isotopic composition of surface waters that are available for drinking is determined by several complex effects, but the oxygen isotopic composition of precipitation (δ 18Oprecip) 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, Reference Dansgaard1964; Rozanski et al., Reference Rozanski, Araguás-Araguás and Gonfiantini1993; Clark and Fritz, Reference Clark and Fritz1997; Kohn and Welker, Reference Kohn and Welker2005; Gat, Reference Gat2010; Stumpp et al., Reference Stumpp, Klaus and Stichler2014; Müller et al., Reference Müller, Stumpp, Sørensen and Jessen2017). Through this effect, the role of δ 18Oenamel as a proxy for palaeotemperatures at mid- to high-latitude sites is derived, where higher δ 18O in tooth enamel indicates higher air temperatures, while low δ 18O values indicate lower temperatures. This is true on longer timescales of climatic change as well as on subannual timescales (Stumpp et al., Reference Stumpp, Klaus and Stichler2014; Müller et al., Reference Müller, Stumpp, Sørensen and Jessen2017), and such subannual δ 18O change is recorded as a sinusoidal seasonal pattern of δ 18O in sequential tooth enamel samples of animal teeth, with peaks indicating the summer season and troughs the winter season (Fricke and O'Neil, Reference Fricke and O'Neil1996; Fricke et al., Reference Fricke, Clyde and O'Neil1998). 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, Reference Dean1987; Hillson, Reference Hillson2005). 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., Reference Brown, Christofferson, Massler and Weiss1960). The use of tooth enamel δ 18O as a palaeotemperature proxy can be formalised for palaeotemperature estimation by establishing species- and site-specific relationships between δ 18Oenamel and δ 18Odw and between δ 18Oprecip and air temperature via linear regression of modern calibration data (Tütken et al., Reference Tütken, Furrer and Walter Vennemann2007; Arppe and Karhu, Reference Arppe and Karhu2010; Skrzypek et al., Reference Skrzypek, Winiewski and Grierson2011, Reference Skrzypek, Sadler and Wiśniewski2016; Pryor et al., Reference Pryor, Stevens, O'Connell and Lister2014). 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 δ 18Oprecip 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 (PO43−) component of tooth enamel (a biogenic carbonate–substituted hydroxyapatite, or “bioapatite”), but the carbonate (CO32−) moiety presents another oxygen-bearing group that is commonly analysed for δ 18O. 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., Reference Iacumin, Bocherens, Mariotti and Longinelli1996; Zazzo et al., Reference Zazzo, Lécuyer, Sheppard, Grandjean and Mariotti2004; Lee-Thorp, Reference Lee-Thorp2008). 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, Reference Pederzani and Britton2019). The use of bioapatite carbonate δ 18O (δ 18Ocarb) for palaeotemperature estimation therefore usually requires an additional conversion step to predicted δ 18O of bioapatite phosphate (δ 18Ophos), which incurs an increased estimation uncertainty that may offset the lower analytical precision of some δ 18Ocarb analyses (Skrzypek et al., Reference Skrzypek, Sadler and Wiśniewski2016).
Due to the complementary advantages of each fraction, a combination of both types of analysis is desirable. A conversion between δ 18Ocarb and δ 18Ophos is possible, as these measures show a strong linear relationship with a slope close to 1 and an intercept offset (δ 18Ocarb-phos) of approximately 9‰ that shows little species-specific variation (Longinelli and Nuti, Reference Longinelli and Nuti1973; Bryant et al., Reference Bryant, Koch, Froelich, Showers and Genna1996b; Iacumin et al., Reference Iacumin, Bocherens, Mariotti and Longinelli1996; Martin et al., Reference Martin, Bentaleb, Kaandorp, Iacumin and Chatri2008; Pellegrini et al., Reference Pellegrini, Lee-Thorp and Donahue2011; Chenery et al., Reference Chenery, Pashley, Lamb, Sloane and Evans2012; Trayler and Kohn, Reference Trayler and Kohn2016). This linear relationship can be used to convert between measurements of the different moieties. While δ 18Ocarb-phos is relatively consistent in empirical studies and matches well with modelled values (Aufort et al., Reference Aufort, Ségalen, Gervais, Paulatto, Blanchard and Balan2017), 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. (Reference Pellegrini, Lee-Thorp and Donahue2011) showed that δ 18Ocarb-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., Reference Pellegrini, Lee-Thorp and Donahue2011; France and Owsley, Reference France and Owsley2015). In contrast to this, Pryor et al. (Reference Pryor, Stevens, O'Connell and Lister2014) 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. (Reference Pryor, Stevens, O'Connell and Lister2014), and the particular concerns related to sequentially sampled δ 18O time-series data were not discussed. To explore this further, we specifically investigate how the uncertainty introduced from converting δ 18Ocarb and δ 18Ophos affects seasonal palaeotemperature estimates in a workflow that includes inverse modelling of oxygen isotopic time series (e.g., Passey et al., Reference Passey, Cerling, Schuster, Robinson, Roeder and Krueger2005) before temperature estimation.
Carbon, nitrogen, and sulphur isotopes
The carbon (δ 13C) and nitrogen (δ 15N) 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, Reference DeNiro and Epstein1978, Reference DeNiro and Epstein1981). For herbivores, this means that these systems reflect isotopic dynamics in the plant biome, where δ 13C varies with plant photosynthetic pathway (C3, C4, 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 CO2 concentration (Park and Epstein, Reference Park and Epstein1961; Tieszen and Boutton, Reference Tieszen, Boutton, Rundel, Ehleringer and Nagy1989; Merwe and Medina, Reference Merwe and Medina1991; Tieszen, Reference Tieszen1991). Plant δ 15N 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, Reference Amundson2003; Stevens and Hedges, Reference Stevens and Hedges2004; Stevens et al., Reference Stevens, Jacobi, Street, Germonpré, Conard, Münzel and Hedges2008, Reference Stevens, Hermoso-Buxán, Marín-Arroyo, González-Morales and Straus2014; Drucker et al., Reference Drucker, Kind and Stephan2011; Craine et al., Reference Craine, Elmore, Wang, Augusto, Baisden, Brookshire, Cramer, Hasselquist, Hobbie and Kahmen2015). These factors create isotopic differences between different plant communities, and through this herbivore δ 13C and δ 15N act as proxies for dietary specialisation. Moreover, the climatic factors influencing δ 13C and δ 15N, as well as the climatic tolerances of different plants and the resulting relationship between landscape change and climate change, mean that plant and herbivore δ 13C and δ 15N also track underlying changes in climate and landscape structure.
The sulphur isotopic composition (δ 34S) of herbivore collagen reflects plant δ 34S with little biological fractionation and, through this, underlying isotopic variability in the geosphere and hydrosphere. Differences in δ 34S 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, Reference Thode, Krouse and Grinenko1991; Nehlich, Reference Nehlich2015). Furthermore, wetland plants can show very distinct δ 34S signatures due to their uptake of sulphides that accumulate in anoxic growing conditions (Guiry et al., Reference Guiry, Noël and Fowler2021). These dynamics create spatial differences in δ 34S of bioavailable sulphur, enabling the use of δ 34S as a movement or location tracer, particularly in coastal regions (McArdle et al., Reference McArdle, Liss and Dennis1998; Nehlich, Reference Nehlich2015; Bataille et al., Reference Bataille, Jaouen, Milano, Trost, Steinbrenner, Crubézy and Colleter2021). However, spatial patterns of bioavailable sulphur, and therefore spatial variability of δ 34S 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 δ 34S in a way that is useful for archaeological applications (Mizota and Sasaki, Reference Mizota and Sasaki1996; Szpak et al., Reference Szpak, Longstaffe, Macdonald, Millaire, White and Richards2019; Bataille et al., Reference Bataille, Jaouen, Milano, Trost, Steinbrenner, Crubézy and Colleter2021). At the same time, an impact of permafrost presence and climatically induced changes in baseline δ 34S have been proposed (Reade et al., Reference Reade, Tripp, Charlton, Grimm, Leesch, Müller and Sayle2020). Specific assignment of spatial ecologies based on δ 34S 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., Reference Britton, Jimenez, Le Corre, Pederzani, Daujeard, Jaouen, Vettese, Tütken, Hublin and Moncel2023). This is particularly the case in coastal environments, where larger δ 34S 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 δ 18Ocarb. Summer peak and winter trough areas of the δ 18O 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 δ 18Ophos 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. (Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018) 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. (Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018) 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. (Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018) are marked with an asterisk (*).
Stable isotope analysis of bioapatites (δ 18Ophos, δ 18Ocarb, δ 13Ccarb)
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. (Reference Coplen, Kendall and Hopple1983) and was then injected through vial septa, and samples were left to react overnight to allow for complete conversion to CO2. The resulting CO2 was then analysed for δ 13C and δ 18O 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: δ 13C = −28.63‰ VPDB and δ 18O = −22.69‰ VPDB), the international carbonatite standard NBS-18 (δ 13C = −5.01‰ VPDB and δ 18O = −23.2‰ VPDB), and an in-house chalk standard (IA-R066, accepted values: δ 13C = 2.33‰ VPDB and δ 18O = −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, δ 13C = 1.95 ‰ VPDB and δ 18O = −2.20 ‰ VPDB), while values of IA-R066 were obtained using NBS-18 and IAEA-CO-1 (Carrara marble, δ 13C = 2.5‰ VPDB and δ 18O = −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 δ 13C = 2.21 ± 0.09‰ VPDB (1 SD, n = 2) and δ 18O = −3.83‰ VPDB (1 SD, n = 2) compare well with the long-term mean values of this standard, which were δ 13C = 2.13‰ VPDB (1 SD, n = 2) and δ 18O = −3.99‰ VPDB. Analytical precision of quality-control standard replicates was better than 0.09‰ for δ 13C and better than 0.14‰ for δ 18O.
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. (Reference Dettman, Kohn, Quade, Ryerson, Ojha and Hamidullah2001) and modified by Tütken et al. (Reference Tütken, Vennemann, Janz and Heizmann2006), as described in Britton et al. (Reference Britton, Fuller, Tütken, Mays and Richards2015). A full description of this protocol can be found in Supplementary Text 2.
Oxygen isotope delta measurements of Ag3PO4 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, δ 18O values are therefore only based on two replicates.
The δ 18O 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 (δ 18O = 21.7 ± 0.3‰, 1 SD; Elemental Microanalysis, Okehampton, UK) and an in-house silver phosphate standard (KDHP.N, δ 18O = 4.2 ± 0.3‰, 1 SD). This in-house standard was obtained by equilibrating a KH2PO4 solution with Leipzig winter precipitation at ca. 140°C for several days, after which the solution was neutralised using a small amount of NH4OH and the phosphate precipitated as silver phosphate by addition of AgNO3 solution. The accepted value of this in-house standard was determined by two-point calibration using B2207 and IAEA-SO-6 (barium sulphate, δ 18O = −11.35 ± 0.3‰, 1 SD; Brand et al., Reference Brand, Coplen, Aerts-Bijma, Böhlke, Gehre, Geilmann and Gröning2009).
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 δ 18O 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., Reference Pucéat, Joachimski, Bouilloux, Monna, Bonin, Motreuil and Morinière2010), 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 (δ 13Ccoll, δ 15Ncoll, δ 34Scoll)
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 δ 13C, δ 15N, and δ 34S 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 (Reference Richards and Hedges1999) 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 N2) scales and quality checked using in-house standards IA-R068 (soy protein, accepted values: δ 13C = −25.22‰ VPDB and δ 15N = 0.99‰ AIR), IA-R038 (l-alanine, accepted values: δ 13C = −24.99‰ VPDB and δ 15N = −0.65‰ AIR), IA-R069 (tuna protein, accepted values: δ 13C = −18.88‰ VPDB and δ 15N = 11.60‰ AIR), and a mixture of IAEA-C7 (oxalic acid, δ 13C = −14.48‰ VPDB) and in-house standard IA-R046 (ammonium sulphate, δ 15N = 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, δ 13C = −10.449‰ VPDB) and IAEA-N-1 (ammonium sulphate, δ 15N = 0.40‰ AIR). During the analysis of Axlor samples, IA-R069 gave values of δ 13C = −18.90 ± 0.04‰ VPDB (1 SD, n = 6) and δ 15N = −0.65 ± 0.03‰ AIR (1 SD, n = 6). Average reproducibility of samples measured in duplicate was 0.04‰ for δ 13C and 0.01‰ for δ 15N.
For sulphur stable isotope analysis, calibration to the VCDT (Vienna Cañon Diablo Troilite) scale and correction for 18O contribution were conducted using in-house standards IA-R061 (barium sulphate, δ 34S = 20.33‰ VCDT), IA-R025 (barium sulphate, δ 34S = 8.53‰ VCDT), and IA-R026 (silver sulphide, δ 34S = 3.96‰ VCDT). Accepted values for these in-house standards were obtained by calibrating against international reference materials NBS-17 (barium sulphide, δ 34S = 5.25‰ VCDT) and IAEA-S-1 (silver sulphide, δ 34S = −0.30‰ VCDT). IA-R061, IAEA-SO-5 (barium sulphate, δ 34S = 20.3‰ VCDT), IA-R068 (soy protein, δ 34S = 5.25‰ VCDT), and IA-R069 (tuna protein, δ 34S = 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 δ 34S = 20.32 ± 0.08‰ VCDT (1 SD, n = 6), IAEA-SO-5 δ 34S = 0.47 ± 0.07‰ VCDT (1 SD, n = 2), IA-068 δ 34S = 5.17 ± 0.15‰ VCDT (1 SD, n = 4), and IA-R069 δ 34S = 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, δ 18Ophos and δ 18Ocarb measurements made on the same sample were first checked against each other and against published data sets of the phosphate to carbonate δ 18O spacing (δ 18Ocarb-phos). The resulting relationship was then used to predict δ 18Ophos for the samples with only δ 18Ocarb available, and then predicted δ 18Ophos values were combined with direct δ 18Ophos measurements in the peak and trough areas to form an aggregated δ 18Ophos seasonal time series as the basis for palaeotemperature estimation (Fig. 2). For direct comparisons of seasonal δ 18Ophos values across archaeological layers and with other sites, direct δ 18Ophos measurements were used (Fig. 3).
Figure 2. (A) δ 18Ocarb and δ 18Ophos 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. [Reference Bryant, Koch, Froelich, Showers and Genna1996b]; Iacumin et al. [Reference Iacumin, Bocherens, Mariotti and Longinelli1996]; Shahack-Gross et al. [Reference Shahack-Gross, Tchernov and Luz1999]; Zazzo et al. [Reference Zazzo, Lécuyer, Sheppard, Grandjean and Mariotti2004]; Miller et al. [Reference Miller, Chenery, Lamb, Sloane, Carden, Atici and Sykes2019]; for discussion of lack of species differences, see Pellegrini et al. [Reference Pellegrini, Lee-Thorp and Donahue2011]). Grey dashed lines represent the predictive interval of the linear model of the literature data. (B) The relationship between δ 18Ocarb and δ 18Ophos was used to predict δ 18Ophos values from δ 18Ocarb measurements and construct a combined δ 18O time series of δ 18Ophos measurements (opaque points) and predicted values (transparent points). ERJ, enamel–root junction.
Figure 3. Summer peak, mean annual, and winter trough δ 18Ophos data points extracted from Bos/Bison sp. seasonal δ 18O 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 δ 18Ophos time series were put into an inverse model following Passey et al. (Reference Passey, Cerling, Schuster, Robinson, Roeder and Krueger2005) to remove the damping of the seasonal δ 18O amplitude that is caused by time averaging through enamel mineralization and homogenization within each incremental sample. This approach estimates the original δ 18O 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. (Reference Passey, Cerling, Schuster, Robinson, Roeder and Krueger2005). 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., Reference Zazzo, Bendrey, Vella, Moloney, Monahan and Schmidt2012; Bendrey et al., Reference Bendrey, Vella, Zazzo, Balasse and Lepetz2015). However, while such considerations have been used to develop updated inverse models for sheep (Green et al., Reference Green, Smith, Green, Bidlack, Tafforeau and Colman2018), 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 (Reference Passey and Cerling2002) and Kohn (Reference Kohn2004). 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 (Emeas) to the prediction error (Epred). 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 δ 18O values. Modelled summer and winter values were extracted in the δ 18O profile locations that correspond to peaks and troughs in the original unmodelled δ 18O profile. The extracted modelled summer peak and winter trough values were used as the input to generate summer and winter δ 18Odw values as well as summer and winter absolute temperatures. Mean annual δ 18O values were computed directly from raw δ 18Ophos 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 δ 18O curve, making little changes to the average δ 18O of an annual cycle (for a detailed discussion, see Pederzani et al., Reference Pederzani, Aldeias, Dibble, Goldberg, Hublin, Madelaine and McPherron2021). Therefore, we used annual averages directly from δ 18Ophos values to avoid the impact of uncertainty on the model output inflating the error ranges on δ 18Odw and palaeotemperature estimates.
Palaeotemperatures were then computed following methods outlined in Pryor et al. (Reference Pryor, Stevens, O'Connell and Lister2014) using the Excel files published therein for computations of temperature estimates and associated uncertainties. Enamel δ 18O to δ 18Odw regressions used in this study were established using ordinary least square regression (OLS) applied to published data from bison (Hoppe, Reference Hoppe2006) and cattle (D'Angela and Longinelli, Reference D'Angela and Longinelli1990). Modern calibration data to establish the relationship between precipitation δ 18O 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 δ 18Oprecip–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; nstations = 75), coldest month (commonly December or January; nstations = 83), and long-term mean annual averages (nstations = 89), respectively. Regression lines for the δ 18Oprecip–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 δ 18Ophos from δ 18Ocarb
As we use the δ 18Ocarb to δ 18Ophos relationship in this study only to convert δ 18Ocarb measurements into δ 18Ophos values specifically for a group of study specimens with some paired data available, we chose to use the δ 18Ocarb to δ 18Ophos relationship that is specifically derived from our study specimens to predict δ 18Ophos from samples that only have δ 18Ocarb 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 δ 18Ophos/δ 18Ocarb 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 δ 18Ocarb to δ 18Ophos in the training data set was used to predict δ 18Ophos for the validation data set, and the results were compared with measured δ 18Ophos.
Beyond pure δ 18Ophos 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 δ 18Ophos from δ 18Ocarb 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 δ 18Ophos, we chose a single δ 18O time series (individual AXL66) and randomly generated 10 alternative combined carbonate/phosphate time-series data sets by applying a random error onto the predicted δ 18Ophos 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 δ 18Ophos 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 (δ 18O, δ 13C, δ 15N, δ 34S) 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., Reference Allaire, Cheng, Xie, McPherson, Chang, Allen and Hyndman2018). 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, Reference Novomestky2021), Metrics_0.1.4 (Hamner and Frasco, Reference Hamner and Frasco2018), ggpubr_0.4.0 (Kassambara, Reference Kassambara2020a), patchwork_1.1.1 (Pedersen, Reference Pedersen2020), magick_2.7.3 (Ooms, Reference Ooms2020), scales_1.2.0 (Wickham and Seidel, Reference Wickham and Seidel2020), rstatix_0.7.0 (Kassambara, Reference Kassambara2020b), officer_0.4.2 (Gohel, Reference Gohel2020a), stringr_1.4.0 (Wickham, Reference Wickham2019), flextable_0.7.0 (Gohel, Reference Gohel2020b), knitr_1.39 (Xie, Reference Xie, Stodden, Leisch and Peng2014), wesanderson_0.3.6 (Ram and Wickham, Reference Ram and Wickham2018), captioner_2.2.3 (Alathea, Reference Alathea2015), cowplot_1.1.1 (Wilke, Reference Wilke2019), tidyr_1.2.0 (Wickham and Henry, Reference Wickham and Henry2019), dplyr_1.0.9 (Wickham et al., Reference Wickham, François, Henry and Müller2019), and ggplot2_3.3.6 (Wickham, Reference Wickham2016).
RESULTS
Radiocarbon dating
A new 14C 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., Reference Reimer, Austin, Bard, Bayliss, Blackwell, Bronk Ramsey, Butzin, Cheng, Edwards and Friedrich2020], OxCal version 4.4 [Bronk Ramsey, Reference Bronk Ramsey2009]; 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‰ δ 13C and 6.4‰ δ 15N 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 δ 18Ocarb and δ 18Ophos, 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., Reference Bryant, Koch, Froelich, Showers and Genna1996b; Iacumin et al., Reference Iacumin, Bocherens, Mariotti and Longinelli1996; Shahack-Gross et al., Reference Shahack-Gross, Tchernov and Luz1999; Zazzo et al., Reference Zazzo, Lécuyer, Sheppard, Grandjean and Mariotti2004; Miller et al., Reference Miller, Chenery, Lamb, Sloane, Carden, Atici and Sykes2019), 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:
while the regression equation of the modern comparative data shows a slope slightly below 1, with:
In a test of δ 18Ophos prediction performance using a linear model of the relationship between δ 18Ocarb and δ 18Ophos from a randomly selected subset training data set of δ 18Ocarb and δ 18Ophos measurement pairs obtained in this study (see section “Predicting δ 18Ophos from δ 18Ocarb”), measured and predicted δ 18Ophos values fall close to a 1:1 relationship with minimal scatter (R2 = 0.87; Supplementary Fig. S2). The prediction shows a root-mean-square error of 0.56‰, approximately twice the measurement uncertainty of δ 18Ophos (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 δ 18O~phos prediction error (Supplementary Fig. S3), we generated random variants of one existing δ 18O time series (AXL66; Supplementary Fig. S4) and evaluated the variability between inverse model outcomes for these within-error prediction alternatives (see “Predicting δ 18Ophos from δ 18Ocarb”); Supplementary Fig. S5). The root-mean-square error of extracted seasonal δ 18O 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 δ 18Ophos 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, δ 18O 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 δ 18Ophos 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 δ 18O 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 δ 18O (P summer = 0.44, P winter = 0.98, P mean annual = 0.75).
Given the limited amount of published δ 18O 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 δ 18O values to drinking water δ 18O values is necessary, using a species-specific relationship. At the same time, this facilitates a comparison with modern measurements of δ 18O in precipitation and environmental waters. The δ 18Odw for the Middle Palaeolithic Bos/Bison sp. from Axlor falls between ~−4‰ and −2‰ in summer and at ~−14‰ in winter, with mean δ 18Odw 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 δ 18Oprecip 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 δ 18Odw estimates from the Axlor data fall within the local variability of modern-day δ 18Oprecip, while winter δ 18Odw estimates are lower than modern δ 18Oprecip 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 δ 18O 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 Ntotal = 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, Reference Van Klinken1999). 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, Reference Nehlich and Richards2009).
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 δ 15N 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 δ 13C 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‰ δ 13C and ~11‰ δ 34S, respectively, but are noticeably distinct from each other in δ 15N, with means of 6.2‰ δ 15N in Bos/Bison sp. and 3.3‰ δ 15N 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 δ 13C, ~5‰ in δ 15N, and ~6‰ in δ 34S (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 δ 13C. 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 δ18O time series for palaeclimatology
The isotopic offset in δ 18O between the carbonate and the phosphate moiety (Δ18Ocarb-phos) observed in Bos/Bison sp. tooth enamel samples in this study of 8.7 ± 0.7‰ (1 SD) closely matches the average Δ18Ocarb-phos offset in published modern comparative data of 8.7 ± 0.9‰ (mean and 1 SD) (Bryant et al., Reference Bryant, Koch, Froelich, Showers and Genna1996b; Pellegrini et al., Reference Pellegrini, Lee-Thorp and Donahue2011; Trayler and Kohn, Reference Trayler and Kohn2016) as well as the theoretically predicted offset of ~8‰ (Aufort et al., Reference Aufort, Ségalen, Gervais, Paulatto, Blanchard and Balan2017). The regression equation in the Axlor sample between δ 18Ocarb and δ 18Ophos 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, Reference Longinelli and Nuti1973; Bryant et al., Reference Bryant, Koch, Froelich, Showers and Genna1996b; Iacumin et al., Reference Iacumin, Bocherens, Mariotti and Longinelli1996; Martin et al., Reference Martin, Bentaleb, Kaandorp, Iacumin and Chatri2008; Pellegrini et al., Reference Pellegrini, Lee-Thorp and Donahue2011; Chenery et al., Reference Chenery, Pashley, Lamb, Sloane and Evans2012; Trayler and Kohn, Reference Trayler and Kohn2016). Overall, the good match between the Axlor samples and the theoretically and empirically documented relationship between δ 18Ocarb and δ 18Ophos indicates that the enamel samples are well preserved and that it may be possible to predict δ 18Ophos from δ 18Ocarb with reasonable accuracy and precision to reduce the amount of costly δ 18Ophos measurements needed.
The consistency of δ 18Ocarb-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 δ 18Ocarb-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., Reference Pellegrini, Lee-Thorp and Donahue2011). 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 δ 18Ophos from δ 18Ocarb for a subset of samples in a sequentially sampled δ 18O 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 δ 18Ophos and δ 18Ocarb measurements, it has been possible to determine that prediction uncertainty of δ 18Ophos 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 δ 18O 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 δ 18O curve do influence the outcome of the inverse modelling procedure but are not themselves immediate areas of interest. In this approach, actual δ 18Ophos measurements were conducted on samples in the summer peak and winter trough areas. This apparently sufficiently serves to anchor these parts of the δ 18O curve and protects areas of interest for palaeotemperature estimation from being substantially influenced by the uncertainty of the δ 18Ophos to δ 18Ocarb prediction. We simulated the effect of prediction uncertainty on δ 18O curve shape and inverse model outcomes by generating alternative δ 18O time series that represent possible variations of the δ 18O curve within the prediction uncertainty of δ 18Ophos from δ 18Ocarb and compared the extracted seasonal δ 18O values of the inverse model outcome between each trial. Very little difference is noted between the summer and winter δ 18O 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 δ 18O curves, where intermediate (non-peak) curve areas have been predicted from δ 18Ocarb, 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 δ 18O time-series data that include actual δ 18Ophos 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 δ 18O 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 δ 18Ophos curves, and we recommend that any δ 18Ophos measurements should be taken from peak and trough areas. In this study, it is demonstrated that three δ 18Ophos 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 δ 18Ocarb vs δ 18Ophos time series has not been consistently observed (Pellegrini et al., Reference Pellegrini, Lee-Thorp and Donahue2011; Trayler and Kohn, Reference Trayler and Kohn2016), but small variations in the location of the most extreme δ 18O 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 δ 18O series, but not necessarily for palaeotemperature estimates from unmodelled time series or bulk samples, as our arguments rest on the preservation of the δ 18O seasonal curve shape through the inverse modelling procedure due to the use of strategic δ 18Ophos anchor points. In case studies where δ 18Ophos 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 δ 18Ocarb to δ 18Ophos prediction uncertainty (for further discussion of bulk samples, see Pryor et al., Reference Pryor, Stevens, O'Connell and Lister2014; Skrzypek et al., Reference Skrzypek, Sadler and Wiśniewski2016).
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, Reference Ríos-Garaizar2017; Gómez-Olivencia et al., Reference Gómez-Olivencia, Sala, Núñez-Lahuerta, Sanchis, Arlegi and Rios-Garaizar2018; Marín-Arroyo et al., Reference Marín-Arroyo, Rios-Garaizar, Straus, Jones, De la Rasilla, González Morales, Richards, Altuna, Mariezkurrena and Ocio2018; González-Urquijo et al., Reference González-Urquijo, Bailey and Lazuen2021; Sánchez Hernández, Reference Sánchez Hernández2021; Table 1). An additional two infinite radiocarbon ages from Layer IV, produced by Marín-Arroyo et al. (Reference Marín-Arroyo, Rios-Garaizar, Straus, Jones, De la Rasilla, González Morales, Richards, Altuna, Mariezkurrena and Ocio2018), 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., Reference González-Urquijo, Bailey and Lazuen2021) 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, Reference Sánchez Hernández2021). 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 δ 18O values in northern Spain, other secondary impacts on δ 18Oprecip, and therefore δ 18Oenamel, 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 C3 ecosystems, δ 18O and δ 13C 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., Reference Drucker, Bridault, Hobson, Szuma and Bocherens2008; Diefendorf et al., Reference Diefendorf, Mueller, Wing, Koch and Freeman2010; Feranec et al., Reference Feranec, García, DÍez and Arsuaga2010; Richards et al., Reference Richards, Pellegrini, Niven, Nehlich, Dibble, Turq and McPherron2017). In the Axlor Bos/Bison sp. samples, however, neither annual means nor sequential data of δ 18O and δ 13C systematically correlate with each other (Supplementary Figs. S7 and S8), indicating that aridity is not a major driver of δ 18O and δ 13C 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, Reference Gonfiantini, Fritz and Fontes1986; Darling et al., Reference Darling, Bath and Talbot2003; Gat, Reference Gat2010; Halder et al., Reference Halder, Terzer, Wassenaar, Araguás-Araguás and Aggarwal2015; Pederzani and Britton, Reference Pederzani and Britton2019). 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., Reference Ocio, Stocker, Eraso and Cowpertwait2016). 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 δ 18O 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 δ 18O values. For this reason, we account for this in our comparisons with modern water isotope data by including δ 18Oprecip estimates from high-elevation locations to represent as much as possible the local variability in δ 18O of environmental waterbodies (Supplementary Fig. S6). Overall, we conclude that Bos/Bison sp. δ 18O 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., Reference Serrano, González-Trueba, Pellitero and Gómez-Lende2017), which could have supplied low δ 18O glacial melt into the local watersheds. Potential low δ 18O inputs, therefore, need to be considered when making comparisons with modern δ 18O 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 δ 18Odw estimates overlap within error with current summer δ 18Oprecip values (Supplementary Fig. S6). This may indicate that the δ 18Oprecip–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 δ 18Oprecip data from Santander indicate a relatively low δ 18Oprecip–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 δ 18Odw 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., Reference Serrano, González-Trueba, Pellitero and Gómez-Lende2017) 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., Reference Laskar, Robutel, Joutel, Gastineau, Correia and Levrard2004). 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., Reference Ballesteros, Rodríguez-Rodríguez, González-Lemos, Giralt, Álvarez-Lao, Adrados and Jiménez-Sánchez2017), 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., Reference Pellitero, Fernández-Fernández, Campos, Serrano and Pisabarro2019). 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, Reference Ríos-Garaizar, Moreno, Conard and Delagnes2015).
For mean annual conditions, a similarity to Holocene climatic conditions is also supported by more direct comparisons of oxygen stable isotope data. Mean annual δ 18Odw 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., Reference Guede, Zuluaga, Ortega, Alonso-Olazabal, Murelaga, Garcia Camino and Iacumin2020). Estimates of δ 18Odw 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 δ 18Odw 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 δ 18O 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 δ 18Ophos = 16.3 ± 0.7‰, mean summer δ 18Ophos = 17.9 ± 0.8‰, mean winter δ 18Ophos = 14.2 ± 1.0‰; Pederzani et al., Reference Pederzani, Aldeias, Dibble, Goldberg, Hublin, Madelaine and McPherron2021). Our reconstructed δ 18Odw values are also slightly higher, but within error of δ 18Odw values calculated from δ 18Ocarb 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., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018). We used only the individual CAS61 as a comparison, as the other individuals from El Castillo do not show clearly sinusoidal δ 18O curves. Axlor δ 18O mean annual values are ~3‰ lower than a Bos/Bison sp. sample from the approximately contemporary site of Valdegoba Cave near Burgos (δ 18O = −8.6 ± 0.9‰) (Feranec et al., Reference Feranec, García, DÍez and Arsuaga2010). This difference closely matches the difference between the two localities in modern-day precipitation δ 18O (~−6‰ in Axlor [Bowen and Revenaugh, Reference Bowen and Revenaugh2003] 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 δ 18O is well reflected in enamel samples also further supports the effectiveness of using δ 18Oenamel as a proxy for water δ 18O. The observed lack of correlation between δ 18O and δ 13C in the Bos/Bison sp. enamel carbonates (Supplementary Fig. S8) or δ 13C and δ 15N in bone collagen in Bos/Bison sp. or Cervus elaphus (Fig. 6) additionally suggests that relatively high δ 18O are indeed more likely related to warm mean annual temperature conditions rather than an aridity effect.
The similarity in δ 18O 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., Reference Sánchez-Hernández, Gourichon, Pubert, Rendu, Montes and Rivals2019). 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., Reference Power, Salazar-García, Rubini, Darlas, Harvati, Walker, Hublin and Henry2018). 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, Reference Ríos-Garaizar2012; Gómez-Olivencia et al., Reference Gómez-Olivencia, Sala, Núñez-Lahuerta, Sanchis, Arlegi and Rios-Garaizar2018). 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 δ 18O variability would be expected. Because the intra-layer variability in δ 18O 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., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018, Reference Jones, Richards, Reade, Bernaldo de Quirós and Marín-Arroyo2019; Marín-Arroyo et al., Reference Marín-Arroyo, Geiling, Jones, González Morales, Straus and Richards2020).
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. 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. 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., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018). 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., Reference Sánchez Goñi, Bard, Landais, Rossignol and d'Errico2013). 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, Reference Deeprose2018). 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., Reference López-García, Cuenca-Bescós, Blain, Álvarez-Lao, Uzquiano, Adán, Arbizu and Arsuaga2011).
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., Reference Gómez-Olivencia, Sala, Núñez-Lahuerta, Sanchis, Arlegi and Rios-Garaizar2018), 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., Reference Gómez-Olivencia, Arceredillo, Álvarez-Lao, Garate, San Pedro, Castaños and Rios-Garaizar2014; Ríos-Garaizar, Reference Ríos-Garaizar2017). 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, Reference Ríos-Garaizar, Moreno, Conard and Delagnes2015). 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., Reference Gómez-Olivencia, Arceredillo, Álvarez-Lao, Garate, San Pedro, Castaños and Rios-Garaizar2014). 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., Reference Gómez-Olivencia, Arceredillo, Álvarez-Lao, Garate, San Pedro, Castaños and Rios-Garaizar2014), 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, Reference Ríos-Garaizar2017). 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 δ 18O 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 δ 18O 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, Reference Nehlich2015; Jones et al., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018; Reade et al., Reference Reade, Tripp, Charlton, Grimm, Leesch, Müller and Sayle2020). 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 δ 34S is relatively poorly understood in the northern Spanish coastal region, but δ 34S values are likely to at least vary with distance from the coast, where higher δ 34S values prevail closer to the coast (McArdle et al., Reference McArdle, Liss and Dennis1998; Nehlich, Reference Nehlich2015; Bataille et al., Reference Bataille, Jaouen, Milano, Trost, Steinbrenner, Crubézy and Colleter2021). Based on comparison with modelled δ 34S values, coastal δ 34S could be ~4–6‰ higher than in valleys further inland, and δ 34S values of fauna at Axlor seem to not capture isozones that are directly adjacent to the coastline (Bataille et al., Reference Bataille, Jaouen, Milano, Trost, Steinbrenner, Crubézy and Colleter2021). However, it should be noted that the isoscape by Bataille et al. (Reference Bataille, Jaouen, Milano, Trost, Steinbrenner, Crubézy and Colleter2021) 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 δ 34S values seen in the fauna at Axlor. As the local and regional δ 34S variability is largely unknown at the moment, this higher variability may also be related to a relatively high degree of baseline δ 34S 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 δ 15N 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. 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. 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., Reference Jones, Richards, Straus, Reade, Altuna, Mariezkurrena and Marín-Arroyo2018, Reference Jones, Richards, Reade, Bernaldo de Quirós and Marín-Arroyo2019; Marín-Arroyo et al., Reference Marín-Arroyo, Geiling, Jones, González Morales, Straus and Richards2020). 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., Reference Pederzani, Aldeias, Dibble, Goldberg, Hublin, Madelaine and McPherron2021), 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 δ 18O 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 δ 15N 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 δ 15N 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 δ 15N values have been observed in grasses and forbs compared with shrubs or trees in arctic and peri-arctic environments (Drucker, Reference Drucker2022). 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., Reference Pederzani, Aldeias, Dibble, Goldberg, Hublin, Madelaine and McPherron2021; Vidal-Cordasco et al., Reference Vidal-Cordasco, Ocio, Hickler and Marín-Arroyo2022). 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 δ 18Ophos from δ 18Ocarb has a negligible effect on summer and winter palaeotemperature estimations, if primary δ 18Ophos measurements are used to anchor summer peak and winter trough areas of the δ 18O 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.
Supplementary material
The Supplementary Material for this article can be found at https://doi.org/10.1017/qua.2023.32.
Data availability statement
Data sets are available on GitHub (https://github.com/ERC-Subsilience/Axlor_paleoclimatic_data).
Acknowledgments
The carbon, nitrogen and sulphur stable isotope analyses and the stable isotope analyses of bioapatite carbonates were funded as part of the ABRUPT project (HAR2017-84997-P) funded by the Ministry of Science, Innovation and Universities and the SUBSILIENCE project funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 818299—ERC-2018-Consolidator), both awarded to ABM-A. SP was supported by the Max Planck Society and the University of Aberdeen during the time of this project, and the oxygen isotope analysis of bioapatite phosphates was funded by the Max Planck Society. Access to the archaeological collections was granted by the Museo de Arqueología de Bizkaia (Basque Government), and initial sampling for carbon and nitrogen isotope analysis was achieved by Hazel Reade funded by the FP7-PEOPLE-2012-CIG- 322112 project) and ABM-A. We appreciate Joseba Rios-Garaizar's advice about Axlor stratigraphy during the sampling process. We thank Ignacio Valera (IBBTEC, University of Cantabria) for kindly allowing the use of his laboratory facilities for collagen extraction. We thank Carlos Revilla Gómez (IBBTEC, University of Cantabria) for laboratory assistance during collagen extraction. Thanks are also due to Manuel Trost (MPI-EVA) for assistance during silver phosphate preparation and to Sven Steinbrenner for assistance with TC/EA-IRMS. KB is supported by a Philip Leverhulme Prize from the Leverhulme Trust (PLP-2019-284).
Author contributions
The project was designed by ABM-A and KB. Specimen identification and documentation were undertaken by ABM-A and LAP. Sampling and sample processing was conducted by JMG, LAP, JRJ, and SP. TC/EA-IRMS analyses were performed by SP. Data analyses and R coding were carried out by SP. SP wrote the paper with input from all authors.
Code availability statement
The R codes used to perform all of the analyses reported in this manuscript are available on GitHub (https://github.com/ERC-Subsilience/Axlor_paleoclimatic_data).
Open access
