Palaeoecology of ungulates in northern Iberia during the Late Pleistocene through isotopic analysis of teeth

During the Late Pleistocene, stadial and interstadial fluctuations affected vegetation, fauna, and human groups that were forced to cope with these pronounced spatial–temporal climatic and environmental changes. These changes were especially abrupt during Marine Isotope Stage (MIS) 3. Here, we reconstruct the climatic trends in northern Iberia considering the stable isotopic composition of ungulate skeletal tissue found in archaeological deposits dated between 80 and 15 ka cal BP. The carbon and oxygen isotopic composition preserved in the carbonate fraction of tooth enamel provides a reliable and high-resolution proxy of the food and water consumed by these animals, which is indirectly related to the local vegetation, environment, and climate, allowing us to estimate palaeotemperatures and rainfall intensity. This study presents new isotope data from 44 bovine, equid, and cervid teeth from five archaeological sites in the Vasco-Cantabrian region (El Castillo, Axlor, Labeko Koba, Aitzbitarte III interior, and El Otero) and one in northeastern Iberia (Canyars), where human evidence is attested from the Mousterian to the Magdalenian. The carbon isotope values reflect animals feeding on diverse C3 plants in open environments and point to differentiated ecological niches for equids and bovines, especially during the Aurignacian in the Vasco-Cantabrian region. Temperature estimations based on oxygen isotopic compositions and rainfall obtained from carbon isotopic compositions indicate colder and more arid conditions than nowadays for the human occupations from the Late Mousterian to the Aurignacian. The contemporary northeastern Iberian site shows slightly lower temperatures related to an arid period when animals mainly grazed in open landscapes. In the Vasco-Cantabrian region, during MIS 2, the Gravettian data reflect a landscape opening, whereas the Magdalenian points to warmer (but still arid) conditions.

1 Introduction

Understanding local and regional climatic variability during the Late Pleistocene in southern Europe is crucial for assessing the potential impact of climate on the adaptation and decline of Neanderthals and the subsequent expansion and resilience of anatomically modern humans during the Upper Palaeolithic (e.g. D'Errico and Sánchez Goñi, 2003; Finlayson and Carrión, 2007; Sepulchre et al., 2007; Staubwasser et al., 2018). During the Late Pleistocene, the climatic records demonstrate stadial and interstadial continuous fluctuations during Marine Isotope Stage 3 (MIS 3; ca. 60–27 ka) and MIS 2 (ca. 27–11 ka). Human groups had to face those episodes, which affected vegetation and fauna to different extents, depending on the region. Northern Iberia is a key study area due to the abundance of well-preserved archaeological caves and rock shelters, where, in the last decade, an updated and multidisciplinary approach has been applied to disentangle how changing environmental conditions affected the subsistence dynamics of Middle and Upper Palaeolithic hominins. Recent chronological, technological, and subsistence studies and ecological reconstructions reveal a more complex regional panorama than previously known (e.g. Sánchez Goñi, 2020; Vidal-Cordasco et al., 2022, 2023; Timmermann, 2020; Klein et al., 2023).

The Vasco-Cantabrian region, located in northwestern Iberia, is subject to the influence of Atlantic climatic conditions, where the impact of the glacial–interglacial oscillations during MIS 3 has recently been evaluated (Vidal-Cordasco et al., 2022). Modelling of traditional environmental proxies (small vertebrates and pollen) associated to archaeo-palaeontological deposits shows a progressive shift in the climatic conditions, with decreasing temperatures and rainfall levels detected during the late Mousterian (Fernández-García et al., 2023). Ecological alterations have been observed in large mammals, such as niche partitioning between horses and cervids (Jones et al., 2018), a decrease in the available biomass for secondary consumers, and, consequently, a reduction in the ungulate carrying capacity (Jones et al., 2018; Vidal-Cordasco et al., 2022). Cold and arid conditions are observed during the Aurignacian and the Gravettian until the onset of MIS 2. Afterwards, during the Last Glacial Maximum (LGM; 23–19 ka), the global climatic deterioration associated with this glacial phase results in colder and more arid conditions in the region, with a predominance of open landscapes. However, this region still provided resources for human survival, acting as a refugium with more humid conditions in comparison to the Mediterranean area (Cascalheira et al., 2021; Fagoaga et al., 2018; Fernández-García et al., 2023; Garcia-Ibaibarriaga et al., 2019; Lécuyer et al., 2021; Posth et al., 2023). By the end of the LGM, a climate amelioration and a moderate expansion of the deciduous forest are documented from the late Solutrean through the Magdalenian (Garcia-Ibaibarriaga et al., 2019; Jones et al., 2021).

In contrast, northeastern Iberia is influenced by the Mediterranean climate. The MIS 3 human settlement in this region has been linked to cooler temperatures and higher rainfall compared to the present but also to climatic fluctuations that are less pronounced compared to the Vasco-Cantabrian region (López-García et al., 2014; Fernández-García et al., 2020; Vidal-Cordasco et al., 2022). Archaeobotanical and small vertebrate evidence indicate relatively stable climatic conditions but also suggest the persistence of open forests during the Middle to Upper Palaeolithic transition, as found in northwestern Iberia (Allué et al., 2018; Ochando et al., 2021). However, certain archaeological records indicate specific climatic episodes, such as increased aridity and landscape opening during Heinrich events 4 and 5 (HE4 and HE5; e.g. Álvarez-Lao et al., 2017; Daura et al., 2013; López-García et al., 2022; Rufí et al., 2018).

These multi-proxy studies have significantly expanded our understanding of the environmental evolution in Iberia, alongside proxies derived from marine core records in Iberian margins (Fourcade et al., 2022; Martrat et al., 2004; Naughton et al., 2007; Roucoux et al., 2001; Sánchez-Goñi et al., 1999, 2009) and other regional palaeoclimatic records sourced from local natural deposits (e.g. Pérez-Mejías et al., 2019; Moreno et al., 2010, 2012; González-Sampériz et al., 2020; Ballesteros et al., 2020). However, the availability of proxies enabling the direct connections between these environmental shifts and human activities remains limited.

In this study, we investigate the palaeoecological and palaeoenvironmental dynamics in northern Iberia during the late Middle and Upper Palaeolithic by measuring the carbon and oxygen isotopic composition of bioapatite carbonates (δ¹³C_carb/δ¹⁸O_carb) preserved in archaeological mammal teeth. These analyses provide high-resolution snapshots of ecological information from animals accumulated during human occupations at the caves. Tooth enamel forms incrementally and does not biologically remodel (Kohn, 2004; Passey and Cerling, 2002), in contrast to other bodily tissue such as bone. The isotope values measured in it reflect the animal diet and water sources consumed during its mineralization, around 1 to 2 years of life for the species included in our study (bovids, equids, and cervids) (e.g. Hoppe et al., 2004; Pederzani and Britton, 2019; Ambrose and Norr, 1993; Luz et al., 1984). The preserved carbon isotope composition relies on animal dietary choices reflecting mainly the type of plant consumed (C3/C4), exposure to light, and humidity levels. Otherwise, the oxygen isotope composition reflects mainly the environmental water consumed by animals, directly by drinking or through diet, which reflects isotopic information derived from water sources and changes in climatic conditions. Both indirectly provide information on the vegetation and climate that allows the estimation of past temperatures, rainfall, and moisture on a sub-annual scale, returning isotopic data of the foraging areas where animals were feeding during tooth formation.

Figure 1 Location of the archaeological sites included in this study. From west to east, in the autonomous community of Cantabria, El Castillo and El Otero; in the Basque Country, Axlor and Aitzbitarte III interior; in Catalonia, Canyars.

Figure 2 Representation of the duration of each archaeological level (dots represent the median values, and bars represent 95 % confidence intervals for ¹⁴C dates and 68 % confidence intervals for ESR and OSL dates) related to technocomplexes in both northwestern (black) and northeastern (green) Iberia and to the δ¹⁸O record from the NGRIP (North Greenland Ice Core Project members, 2004; Rasmussen et al., 2014). Grey bands indicate the Greenland Stadial (GS). Dates from El Castillo (C14 UF, ESR), El Otero (C14 UF), Axlor (C14 UF, OSL), Labeko Koba (C14 UF), Aitzbitarte III interior (C14 AMS), and Canyars (C14 UF, ABA, ABOx-SC) are shown in Supplement Sects. S2 and S3.

By analysing the stable isotopic composition of 44 ungulate teeth obtained from 15 archaeological levels directly associated with human occupation, including El Castillo, Axlor, Labeko Koba, Aitzbitarte III interior, and El Otero in northwestern Iberia and Terrasses de la Riera dels Canyars in northeastern Iberia, this study presents novel insights into local and regional environmental and climatic trends associated to human presence during the Late Pleistocene (Figs. 1 and 2; Sect. S1). Specifically, it focuses on the Middle to Upper Palaeolithic transition in both areas and on the post-LGM period in the Vasco-Cantabrian region.

Table 1 Number of teeth sampled by species, archaeological sites, and cultural periods.

The main objectives of this work are (1) to assess how regional environmental conditions, including changes in moisture and vegetation cover but also temperatures and rainfall, are recorded in the stable isotopic composition of tooth enamel; (2) to characterize animal diets and their ecological niches; (3) to obtain quantitative temperature data to compare with available proxies; and (4) to characterize seasonal patterns of animals found in the archaeological sites by identifying winter and summer fluctuations.

Table 2 Mean, maximum value (Max), minimum value (Min), and standard deviation (SD) of δ¹³C and δ¹⁸O values per archaeological site and level organized by cultural periods. CCE: calcium carbonate equivalent; n: number of intra-tooth subsamples measured. In tooth type: position (U: upper; L: lower), laterality (R: right; L: left), tooth (M: molar; P: premolar).

2 Archaeological sites and sampled material

This study selected a total of 44 ungulate teeth, including 25 bovines (Bos primigenius, Bison priscus, and Bos/Bison sp.), 14 equids (Equus sp. and Equus ferus), and 5 cervids (Cervus elaphus) originating from five archaeological sites in the Vasco-Cantabrian region (El Castillo, El Otero, Axlor, Labeko Koba, and Aitzbitarte III interior) and one in the northeastern area (Terrasses de la Riera dels Canyars, henceforth Canyars). These teeth were recovered from 15 archaeological levels attributed to the following technocomplexes: Mousterian (n = 14), Transitional Aurignacian (n = 10), Châtelperronian (n = 2), Aurignacian (n = 12), Gravettian (n = 1), and Magdalenian (n = 5) (Tables 1 and 2). Archaeozoological studies of the archaeological sites are available (synthesis in Marín-Arroyo and Sanz-Royo, 2022; Daura et al., 2013), and most prove that faunal remains were accumulated by human acquisition during the different cultural phases. Isotopic results for equid teeth and for the bone collagen of other ungulates from El Castillo were previously published by Jones et al. (2019) as well as the combined bioapatite carbonate and phosphate analyses of bovines from Axlor (Pederzani et al., 2023). A comprehensive description of each archaeological site is provided in Sect. S1.

3 Methods

3.1 Dating methods

Individual Bayesian age models were built for Canyars, El Castillo, Labeko Koba, and Aitzbitarte III interior based on radiocarbon dates (AMS UF and non-UF, ABOx-SC, and ABA pretreatments on bones and charcoal remains) using the software OxCal4.4 (Ramsey, 2009) and considering the IntCal20 calibration curve (Reimer et al., 2020) (Sect. S3). The Bayesian model enables the modification of the calibrated probability distribution function (PDF) of individual dates based on the existing relative stratigraphic and other relative age information. A resolution of 20 years was assumed, being a reasonable balance between required accuracy and computational costs. An order function in OxCal was used to calculate the probability that one PDF predated another, providing information to assess synchronicity and temporal overlap of individual archaeological levels and cultural phases in each of the four separate sites modelled. Dates were organized into a “sequence”, and chronological information for each level was grouped into a single “phase” with start and end “boundaries” to bracket each archaeological level. The interval between the start of each level and its end provided the duration of each level. In all cases, convergence was greater than 95 %. CQL codes, individual Bayesian models, and modelled dates per site are reported in Sect. S3.

No chronological models were built for El Otero because only a single date was obtained for level IV and for El Castillo levels 20E and 21 (ESR dated) and Axlor levels III, IV, and VI (OSL dated) because dates go beyond the limit of the radiocarbon. To show the duration of these levels in combination with the other sites and levels, each of these dates was estimated by adding and subtracting the sigma (68 % confidence interval) from the uncalibrated date. In this way, we estimated the duration of levels that go beyond 55 ka cal BP (Sect. 2).

3.2 Tooth sampling

All teeth included were sequentially sampled to reconstruct the complete δ¹⁸O_carb and δ¹³C_carb intra-tooth profiles based on enamel carbonate bioapatite. Intra-tooth sequential sampling was applied to the second and third molars and to the third and fourth premolars. Bovine and horse teeth sampled exceeded 3–4 cm of crown height to ensure that at least a 1-year isotopic record of animal life was obtained (Britton et al., 2019; Hoppe et al., 2004). Samples were taken perpendicular to the growth axis on the tooth where the enamel was best preserved, avoiding taphonomic alterations, such as cracks or post-depositional damages, whenever possible. Samples were performed in the buccal face for the lower teeth and in the lingual part for the upper ones. The outermost enamel surface was abraded to remove the superficial enamel, calculus, cementum, and concretions adhering to the surface and to avoid contaminations. The sequential sampling consisted of straight strips (ca. 8 mm × 1.5 mm × 1 mm) covering the width of the selected lobe, approximately every 2–3 mm, from the crown to the enamel–root junction (ERJ). The sample depth covered around 75 % of the enamel depth, and dentine inclusion was avoided. A low-revolution variable-speed manual drill was used, equipped with 1 mm diamond-coated drill bits of conical and cylindrical shape. About 10–15 mg of enamel powder was collected in each subsample, generating 693 subsamples for isotope-ratio mass spectrometer (IRMS) measurements (see complete intra-tooth profiles in Sect. S4).

3.3 Sample treatment and stable isotope mass spectrometry

Several authors have debated the necessity of chemical pretreatments to remove organic matter and secondary carbonates from bioapatite carbonates before stable isotopic analysis. Some chemical treatments can introduce secondary carbonates, increase carbonate content, and alter the original isotopic signal (Pellegrini and Snoeck, 2016; Snoeck and Pellegrini, 2015). For this reason, in this work, most of the samples were not pretreated, except for the equid and cervid samples from Labeko Koba, El Otero, and El Castillo that were sampled and pretreated in an earlier phase of the project. The absence of pretreatment can elevate the risk of secondary carbonates (Chesson et al., 2021; France et al., 2020). Nonetheless, no pretreatment method can guarantee their complete removal, and the “side effects” may compromise the final isotopic signal to a greater extent. While variations in pretreatment methods exist among samples in this study, the lack of a universally accepted protocol necessitates careful consideration of any potential isotopic effects resulting from these differences.

Pretreatment was followed for abovementioned samples from 14 teeth, where around 7 mg of powdered enamel was prepared and pretreated with 3 % sodium hypochlorite (NaOCl) at room temperature for 24 h (0.1 mL mg⁻¹ sample) and thoroughly rinsed with deionized water before a reaction with 0.1 M acetic acid for 4 h (0.1 mL mg⁻¹ sample) (Balasse et al., 2002; Jones et al., 2019). Samples were then thoroughly rinsed, frozen, and freeze-dried. NaOCl is one of the most common agents used for pretreating carbonates and works as a base that removes organic matter by oxidation. Although it is regarded as one of the most efficient agents for removing organic matter, it can induce the absorption of exogenous carbonates, such as atmospheric CO₂ and secondary carbonates (Pellegrini and Snoeck, 2016; Snoeck and Pellegrini, 2015). It is argued that acetic acid after NaOCl pretreatment can remove exogenous carbonates absorbed during NaOCl application. However, it is unclear if all newly introduced carbonates are finally released and what effect they have on the original isotopic composition. These samples were analysed in the Godwin Laboratory (Department of Earth Sciences, University of Cambridge). Enamel powder samples were reacted with 100 % orthophosphoric acid for 2 h at 70 °C in individual vessels in an automated GasBench interfaced with a Thermo Finnigan MAT 253 IRMS. Results were reported in reference to the international standard VPDB and calibrated using the NBS 19 standard (limestone: δ¹³C = +1.95 ‰, δ¹⁸O = −2.2 ‰; Coplen, 2011), for which the precision is better than 0.08 ‰ for δ¹³C and 0.11 ‰ for δ¹⁸O.

For the non-pretreated samples, carbon and oxygen stable isotopic ratios were measured using continuous-flow isotope-ratio mass spectrometry, specifically a Europa Scientific 20-20 IRMS coupled to a chromatograph, at the Iso-Analytical laboratory in Cheshire, UK. The samples were weighed into clean exetainer tubes after being flushed with 99.995 % helium. Phosphoric acid was then added to the samples, and they were allowed to react overnight to ensure the complete conversion of carbonate to CO₂, following the method outlined by Coplen et al. (1983). The reference materials used for VPDB calibration and quality control of the analysis included IA-R022 (calcium carbonate: δ¹³C = −28.63 ‰, δ¹⁸O = −22.69 ‰), NBS 18 (carbonatite: δ¹³C = −5.01 ‰, δ¹⁸O = −23.2 ‰), and IA-R066 (chalk: δ¹³C = +2.33 ‰, δ¹⁸O = −1.52 ‰). The accepted values of the in-house standards IA-R022 and IA-R066 were obtained by calibrating against IAEA international reference materials (1) NBS 18 and NBS 19 and (2) NBS 18 and IAEA-CO-1 (Carrara marble: δ¹³C = 2.5 ‰, δ¹⁸O = −2.4 ‰), respectively. Additionally, in-house standards of long-term measurements were used: ILC1 (calcite: δ¹³C = 2.13 ‰, δ¹⁸O = −3.99 ‰) and Y-02 (calcite: δ¹³C = 1.48 ‰, δ¹⁸O = −9.59 ‰). The analytical precision of quality control standard replicates was better than 0.09 ‰ for δ¹³C and better than 0.12 ‰ for δ¹⁸O. The calcium carbonate content test of these samples, ranging between 3.9 % and 8.9 %, does not indicate a substantial presence of secondary carbonates, as per Chesson et al. (2021). Additionally, phosphate results in samples from Axlor showed δ¹⁸O_carb-δ¹⁸O_phos offsets within the expected range for well-preserved samples (Pederzani et al., 2023).

3.4 Carbon stable isotopic compositions as environmental and ecological tracers

To unravel animal diet and compare the different species in standardized terms, it is necessary to consider the enrichment factor (ε^∗) between δ¹³C obtained by the animal from its diet (δ¹³C_diet) and δ¹³C recorded in enamel carbonates (δ¹³C_carb) (Bocherens, 2003; Cerling and Harris, 1999). The ε^∗ estimated for large ruminant mammals results in an offset of around 14.1 ‰ between diet and dental enamel, commonly applied to medium-sized herbivores. However, it is well known that this offset varies between species, considering animals' different physiological parameters. Recently, a formal model to predict species-specific diet–consumer isotopic offsets has been proposed, which uses body mass (BM) and digestive physiology as the main factors that regulate the ε^∗ (Tejada-Lara et al., 2018). This model proposes the following prediction equations for ruminant or foregut fermenters (Eq. 1) and hindgut fermenters (Eq. 2):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\epsilon }}^{\ast }=\mathrm{2.34}+\mathrm{0.05}\left(\text{BM}\right)[r^{\mathrm{2}}=\mathrm{0.78};p\text{-value }=\mathrm{0.008}]\text{(2)}& {\displaystyle }{\mathit{\epsilon }}^{\ast }=\mathrm{2.42}+\mathrm{0.032}\left(\text{BM}\right)[r^{\mathrm{2}}=\mathrm{0.74};p\text{-value}=\mathrm{0.003}].\end{array}$$

In this work, we compare species with different digestive physiology: ruminants for bovines and cervids and non-ruminants for equids. The ε^∗ value was adjusted for each animal to avoid bias from digestive physiology when comparing these species. The following enrichment factors have been used: 14.6 ‰ for Bos taurus (Passey et al., 2005a), 13.7 ‰ for Equus caballus (Cerling and Harris, 1999), and 13.2 ‰ for Cervus elaphus (Merceron et al. (2021), following Eq. (1) for ruminants with a mean body mass of 125 kg.

In body tissue, carbon isotopic composition is regarded as a combination of diet (understood as consumed food), environment openness (and associated exposure to light), and the amount of precipitation. Assuming that δ¹³C of past vegetation is close to δ¹³C_diet of ungulates, Lécuyer et al. (2021) proposed to estimate mean annual precipitation (MAP) from δ¹³C_carb, derived from diets based on C3 plants. After transforming δ¹³C_carb to δ¹³C_diet using the enrichment factors established above, this work suggested transforming this value to δ¹³C from vegetation (δ¹³C_leaf). However, the isotopic composition of animals' diets may not directly reflect vegetation cover but rather the food preference of the animal, and this approach should be discussed alongside other environmental data.

The MAP estimation is based on a least-squares regression developed by Rey et al. (2013) and on the dataset of Kohn (2010) (Eq. 4), which first requires estimating the δ¹³C_leaf (Eq. 3). The δ¹³C values of atmospheric CO₂ (δ¹³C_atm) are fixed at −7 ‰ (Lécuyer et al., 2021; Leuenberger et al., 1992; Schmitt et al., 2012). Atmospheric CO₂ levels varied throughout the Late Pleistocene, with δ¹³C_atm ranging between −7 ‰ and −6.4 ‰ (Eggleston et al., 2016), favouring an age-specific correction approach. However, maintaining general corrections is preferred, considering the chronological uncertainty of the studied levels.

$$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{leaf}}\left(\text{VPDB}\right)=& ({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{atm}}-{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{diet}})\\ & /[\mathrm{1}+({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{diet}}/\mathrm{1000}\left)\right]\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}\text{log1}(\text{MAP}+\mathrm{300})=& \mathrm{0.092}(\pm \mathrm{0.004})\times {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{leaf}}\\ & +\mathrm{1.148}(\pm \mathrm{0.074})\end{array}\end{array}$$

Additionally, the equation of Lécuyer et al. (2021) also accounts for the pCO₂ effect on δ¹³C_leaf estimation, which is expected to result in an offset of +1 ‰ from current levels (considering that the pCO₂ was lower than that experienced after the deglaciation period). If this correction was not applied, MAP results could be underestimated by −150 mm. In agreement with the appreciation of Lécuyer et al. (2021), these MAP estimations are a preliminary approximation and should be cross-validated with other environmental proxies. The associated uncertainties range from ± 100 to 200 mm, influencing the interpretation of the final values.

3.5 Oxygen stable isotope compositions as environmental tracers

Stable oxygen isotopes from meteoric water (mainly derived from rainfall) strongly correlate with mean air temperatures in mid- to high latitudes (Dansgaard, 1964; Rozanski et al., 1992) on a regional to local scale. Obligate drinkers, like bovines and horses, acquire this water and record its isotopic composition in their teeth and bones with a fixed but species-specific offset (Pederzani and Britton, 2019). Considering this two-step relationship, past climatic conditions can be estimated. However, most of the temperature reconstructions based on δ¹⁸O have considered the δ¹⁸O from the phosphate fraction of bioapatite enamel (δ¹⁸O_phos) in order to build linear correlations between tooth enamel and drinking water δ¹⁸O and obtain climatic information. For this reason, the δ¹⁸O_carb values obtained in this work were converted into δ¹⁸O_phos. To do so, in order to express in VSMOW notation, the δ¹⁸O_carb was firstly corrected using the following correlation (Brand et al., 2014; Coplen et al., 1983):

$$\begin{array}{}\text{(5)}& {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{carb}}\left(\text{VSMOW}\right)=\mathrm{1.0309}\times {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{carb}}\left(\text{VPDB}\right)+\mathrm{30.91}.\end{array}$$

Secondly, considering the relationship existent in tooth enamel between the carbonate and phosphate fraction (Iacumin et al., 1996; Pellegrini et al., 2011) from a compilation of the existent bibliography of modern animal measurements (Bryant et al., 1996; Pellegrini et al., 2011; Trayler and Kohn, 2017), Pederzani et al. (2023) proposed the following correlation:

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{phos}}\left(\text{VSMOW}\right)=\mathrm{0.941}\times {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{carb}}\left(\text{VSMOW}\right)-\mathrm{7.16}.\end{array}$$

Once the isotopic information is expressed in δ¹⁸O_phos (VSMOW), we can estimate the δ¹⁸O in meteoric waters (δ¹⁸O_mw). It is known that different physiological factors will condition how oxygen isotope composition is fixed in each mammalian group. Thus, the correlations are usually species-specific and developed considering the physiology of each animal group. The obligate drinkers heavily rely on consuming large amounts of liquid drinking water, with the relative contribution of water from plants being negligible, and then minimizing the possible impact of isotopic enrichment through evapotranspiration in plants (Hoppe, 2006; Maloiy, 1973; Pederzani and Britton, 2019). However, certain types of drinking behaviours can impact δ¹⁸O, such as systematic consumption of certain highly buffered water sources (rivers or lakes), and can significantly attenuate the final signal recorded. The correlation employed by this work relies on recent data compilations (Pederzani et al., 2021b, 2023). In the case of horses (Eq. 7), the data combination of Blumenthal et al. (2019), Chillón et al. (1994), Bryant et al. (1994), and Delgado Huertas et al. (1995) is considered; for bovines (Eq. 8), the data from D'Angela and Longinelli (1990) and Hoppe (2006) are put together in Eq. (8). To estimate δ¹⁸O_mw from red deer remains, we selected the correlation of D'Angela and Longinelli (1990) (Eq. 9).

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{mw}}\left(\text{VSMOW}\right)=& \left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{phos}}\right(\text{VSMOW})-\mathrm{22.14})\\ & /\mathrm{0.62}\end{array}\text{(8)}& {\displaystyle }\begin{array}{rl}{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{mw}}\left(\text{VSMOW}\right)=& \left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{phos}}\right(\text{VSMOW})-\mathrm{22.36})\\ & /\mathrm{0.78}\end{array}\text{(9)}& {\displaystyle }\begin{array}{rl}{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{mw}}\left(\text{VSMOW}\right)=& \left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{phos}}\right(\text{VSMOW})-\mathrm{24.39})\\ & /\mathrm{0.91}\end{array}\end{array}$$

Finally, palaeotemperature estimations from δ¹⁸O_mw are typically approached using a geographically adjusted linear regression, which can vary from precise adjustments (aimed at reducing errors) to broader geographical adjustments that encompass more variability but are less precise (e.g. Pryor et al., 2014; Skrzypek et al., 2011; Tütken et al., 2007). In this work, temperatures were calculated considering the linear regression model relating δ¹⁸O_mw and air temperatures proposed by Pederzani et al. (2021b) based on monthly climatic records (monthly mean δ¹⁸O_mw and monthly mean air temperatures) from western, southern, and central European stations from the Global Network of Isotopes in Precipitation (IAEA/WMO, 2020). Considering current IAEA datasets from northern Iberia, there is a strong positive relationship between δ¹⁸O_mw and annual or monthly temperatures (Moreno et al., 2021). However, it is known that Iberia is under a mixed influence of Atlantic and Mediterranean moisture sources that affects the isotopic composition of rainfall (Araguas-Araguas and Diaz Teijeiro, 2005; García-Alix et al., 2021; Moreno et al., 2021). Given uncertainties in past atmospheric circulation patterns and the limited availability of reference stations, it was deemed most appropriate to select an equation that extends beyond the borders of Iberia and incorporates higher variability. Different correlations were for mean annual (Eq. 10), summer (Eq. 11), and winter (Eq. 12) temperatures (T).

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{mw}}\left(\text{VSMOW}\right)=(\mathrm{0.50}\times T)-\mathrm{13.64}\text{(11)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{mw}}\left(\text{VSMOW}\right)=(\mathrm{0.46}\times T)-\mathrm{14.70}\text{(12)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{mw}}\left(\text{VSMOW}\right)=(\mathrm{0.52}\times T)-\mathrm{11.26}\end{array}$$

Nonetheless, oscillations between glacial and interglacial conditions in the past have influenced global ice volume and sea level fluctuations (Dansgaard, 1964; Shackleton, 1987), impacting seawater oxygen isotope composition and the surface hydrological cycle on a worldwide scale, including δ¹⁸O_mw (Schrag et al., 2002). Prior studies used sea level information to correct δ¹⁸O_mw (e.g. Fernández-García et al., 2019; Schrag et al., 2002). Given the chronological uncertainty in the studied levels, a general correction was applied to δ¹⁸O_mw before temperature estimations, following the approach of Fernández-García et al. (2020). Considering the mean sea level descent for MIS 3 (50 m below present-day sea level) (Chappell and Shackleton, 1986), this may have contributed to a potential increase in the global δ¹⁸O_mw value by ≈ 0.5 ‰, implying a bias in calculated air temperatures of ≈ 1 °C.

Due to the uncertainties incurred when converting stable isotope measurements to palaeotemperature, the final estimations in this work should be regarded as exploratory and as a method of standardization to make results comparable among different sites, species, and other non-isotopic palaeoclimatic records. In these estimations, the associated error from converting δ¹⁸O_phos to mean annual temperature (MAT) is enlarged by the uncertainty derived from the transformation of δ¹⁸O_carb (VPDB) to δ¹⁸O_phos (VSMOW) (see Pryor et al., 2014; Skrzypek et al., 2016 for further discussion). However, Pryor et al. (2014) and Pederzani et al. (2023) concluded that the impact of this conversion is negligible compared to the error propagation in subsequent calibrations used for temperature estimations from δ¹⁸O_phos. These associated errors were quantified following the methodology outlined by Pryor et al. (2014) (Sect. S2).

3.6 Inverse modelling applied to intra-tooth profiles

Intra-tooth profiles frequently provide a time-averaged signal compared to the input isotopic signal (δ¹³C/δ¹⁸O) during enamel formation (Passey et al., 2005b). This signal attenuation is caused by time-averaging effects incurred through the extended nature of amelogenesis and tooth formation and through the sampling strategy. During mineralization, the maturation zone, which is time-averaged, often affects a large portion of the crown height and might affect the temporal resolution of the input signal of the sample taken. To obtain climatically informative seasonal information on the analysed teeth, the inverse modelling method proposed by Passey et al. (2005b) is applied in this work. This method computationally estimates the time-averaging effects of sampling and tooth formation to obtain the original amplitude of the isotopic input signal more accurately, thus, to summer and winter extremes (Sect. S5). This method considers parameters based on the amelogenesis trends of each species and sampling geometry, which are critical for a meaningful interpretation of intra-tooth isotope profiles. The model also estimates the error derived from the sampling uncertainty and the mass spectrometer measurements to evaluate the data's reproducibility and precision. This method was initially developed for continuously growing teeth, taking into account a constant growth rate within a linear maturation model, with a progressive time-averaged increment as sampling advances along the tooth profile. The species studied in this research exhibit non-linear tooth enamel formation, particularly in later-forming molars (Bendrey et al., 2015; Blumenthal et al., 2014; Kohn, 2004; Passey and Cerling, 2002; Zazzo et al., 2012). Although the model mentioned above is not ideal, as it does not take into account non-linear enamel formation and as specific growth parameters for the species included are unknown, it is the best estimation based on the current state of the field and remains widely used (Pederzani et al., 2021a, b, 2023). Flat and less sinusoidal profiles are less suitable for the application of the model, given its inherent assumption of an approximately sinusoidal form. Therefore, we chose not to apply this methodology in the analysis of intra-tooth δ¹³C profiles, and it is recommended to approach the interpretation of model outcomes for non-sinusoidal δ¹⁸O curves with caution. Further details on the application of this method can be found in Sect. S5.

Following Pederzani et al. (2021b), mean annual temperatures (MATs) were deduced from the average of δ¹⁸O_carb values between summer and winter detected in original sinusoidal intra-tooth profiles (Sect. S5). This work shows that comparable results for annual means can be obtained before and after model application, but doing it beforehand avoids the associated errors induced by the inverse model. To maximize data, in non-sinusoidal teeth profiles, MAT was deduced from the average of all points within a tooth. However, this approach is less reliable when complete annual cycles are not recorded. When possible, summer and winter temperature estimations were derived from the obtained δ¹⁸O_carb values after inverse modelling application, aiming to identify the corrected seasonal amplitude, which is dampened in the original δ¹⁸O_carb signal.

3.7 Present-day isotopic and climatic data

Present-day climatic conditions surrounding each site have been considered, allowing an inter-site comparison, essential for comparing this study with other regional and global data. Considering current MATs and MAPs, estimated climatic data are expressed in relative terms as MAT and MAP anomalies (MATA and MAPA, respectively). Present-day summer and winter temperatures were also considered. Present-day temperatures and precipitation values were obtained from the WorldClim dataset v2 (Fick and Hijmans, 2017) (Sect. S2). This dataset includes the average of bioclimatic variables between 1970–2000 in a set of raster files with a spatial resolution every 2.5 min. The exact location of the selected archaeo-palaeontological sites was used, using geographical coordinates in the projection on modern climatic maps with QGIS software.

Present-day δ¹⁸O_mw values from the analysed sites were obtained using the Online Isotopes in Precipitation Calculator (OIPC Version 3.1 (4/2017); Bowen, 2022) based on datasets collected by the Global Network of Isotopes in Precipitation from the IAEA/WMO (Sect. S2).

Table 3 Mean δ¹³C from enamel carbonate (δ¹³C_carb) and diet (δ¹³C_diet), along with δ¹⁸O from enamel carbonate (δ¹⁸O_carb) and meteoric waters (δ¹⁸O_mw), from species on the Vasco-Cantabrian and northeastern Iberia areas. Max: maximum value; Min: minimum value; SD: standard deviation.

4 Results

In northwestern Iberia, specifically in the Vasco-Cantabrian region, the mean δ¹³C_carb values range from −13 ‰ to −8.9 ‰, with a mean value of −11 ‰ (SD = 1.2 ‰) (Tables 2 and 3). Considering species' different enrichment factors, the δ¹³C_carb values were transformed into δ¹³C_diet, resulting in mean values that extend from −27 ‰ to −23.5 ‰ (Fig. 4). It must be considered that average values may reflect slightly different periods or be affected by seasonal bias because different teeth encompass diverse periods, but it has been verified in our teeth that the variations are limited when the seasonal information of the sequential sampling is incorporated (± 0.2; Sect. S2). The carbon isotopic composition varies between species. The bovines have generally higher mean δ¹³C_carb (from −12.4 ‰ to −8.9 ‰) than the horses (from −12.6 ‰ to −11.3 ‰), whereas the red deer fall within the horses' range (from −13 ‰ to −11.3 ‰). Average values of δ¹⁸O_carb in all Vasco-Cantabrian individuals extend between −7.2 ‰ and −3.3 ‰ (mean = −5.5 ‰; SD = 0.8 ‰). When transformed to δ¹⁸O expected from meteoric waters (δ¹⁸O_mw), with species-adapted correlations, the δ¹⁸O_mw values range from −10.6 ‰ to −5.5 ‰. Less clear patterns in δ¹⁸O_carb are observed between bovines and horses, with mean values of −5.7 ‰ and −5.2 ‰, respectively. In northeastern Iberia, the site of Canyars, both species have relatively high δ¹⁸O_carb values that fall inside the range of variation observed in the Cantabria region, between −5.5 ‰ and −3.6 ‰ in bovines and between −4.8 ‰ and −4.4 ‰ in horses.

Figure 3 Distribution of mean carbon (δ¹³C_carb) and oxygen (δ¹⁸O_carb) isotopic values of enamel carbonate by species and archaeological site.

Figure 4 Biplot crossing δ¹³C from diet (δ¹³C_diet) and δ¹⁸O from meteoric waters (δ¹⁸O_mw) by species and archaeological site.

Table 4 Summary of palaeoclimatic estimations, based on δ¹⁸O for temperatures (mean annual temperature, MAT; summer; winter) and on δ¹³C for precipitation (mean annual precipitation, MAP). Summer and winter temperature estimations were obtained from teeth with clear seasonal profiles after modelling, while MAT was averaged between summer and winter before modelling. In profiles with an unclear seasonal shape, MAT was deduced from the original average of all tooth points (values marked in italics). Mean error associated to temperature estimations is 5.1 ± 0.6 (see details in Sect. S2). Seasonality is calculated as the temperature difference between summer and winter.

4.1 Axlor (Mousterian, ca. 80 ka BP–50 ka cal BP)

A total of seven bovine teeth were included from levels III (n = 4), IV (n = 1), and VI (n = 2) of Axlor cave (Pederzani et al., 2023). The mean δ¹³C_carb ranges from −9.9 ‰ to −8.9 ‰ (δ¹³C_diet = −24.5 ‰ to −23.5 ‰), whereas mean δ¹⁸O_carb values are between −6.2 ‰ and −4.8 ‰ (δ¹⁸O_mw = −8.3 ‰ and −6.5 ‰), indicating a range of variation around 1 ‰ and 1.4 ‰, respectively (Figs. 3 and 4). Considering isotopic compositions by levels, mean δ¹³C_carb decreases from level III to level IV, whereas mean δ¹⁸O_carb remains stable through the sequence (Table 2; Sect. S2). A range between 0.3 ‰ and 0.5 ‰ is observed in δ¹³C_carb variation within tooth profiles. Individuals show clear δ¹⁸O sinusoidal profiles, with peaks and troughs and intra-tooth ranges from 2.1 ‰ to 3.4 ‰. The δ¹⁸O_mw after inverse modelling intra-tooth profiles ranges from −9.1 ‰ to −7.35 ‰ (Sects. S4 and S5). Mean annual temperatures (MATs) oscillated between 9.1 and 12.6 °C (MATA = −3.1/+0.4 °C) (Table 4). From sinusoidal profiles, summer temperatures were extracted from peaks, ranging from 15.4 to 23.7 °C, and winter temperatures from troughs provided values ranging from −7 to 10.8 °C. Mean annual precipitation (MAP), extracted from δ¹³C_carb, extends between 204 and 326 mm (MAPA = −843/−721 mm). Based on these estimations, a non-clear climatic trend is observed through these levels.

Figure 5 Evolution of δ¹³C in diet (δ¹³C_diet) and δ¹⁸O in meteoric waters (δ¹⁸O_mw) by archaeological levels in a diachronic order. From right to left: all species, including cervids, bovines and horses. Colours correspond to different chronocultures.

4.2 El Castillo (Mousterian and Transitional Aurignacian, 62.5 ka BP–46.4 ka cal BP)

From El Castillo, this work includes bovine (n = 11), horse (n = 5), and red deer (n = 1) teeth from Mousterian (21 and 20E) and Transitional Aurignacian levels (18B and 18C). The mean δ¹³C_carb values are lower for horses, bovines, and red deer (−13 ‰ to −10.2 ‰) than in other sites: between −12.4 ‰ and −10.2 ‰ for bovines (δ¹³C_diet = −24.6 ‰ to −25.8 ‰) and between −12.6 ‰ and −11.5 ‰ for horses (δ¹³C_diet = −26.3 ‰ to −25.2 ‰) (Fig. 3). The mean δ¹⁸O_carb values extend from −6.8 ‰ and −3.3 ‰. Horses and bovines overlap in their isotopic niche (Fig. 4), mainly due to the notably lower δ¹³C_carb reported by bovines. The mean δ¹³C_carb (−13 ‰) of the single red deer tooth is inside the variation range of bovines and horses but with a lower δ¹⁸O_carb mean value (−6.8 ‰). Considering these isotopic compositions by levels, bovine mean δ¹³C_diet values highly increase the variation range from Mousterian levels (20E and 21A) to Transitional Aurignacian levels (18C and 18B). In contrast, horses increase mean δ¹³C_diet values (Fig. 5). Bovine mean δ¹⁸O_mw values decrease from level 21A to level 18B, while horses from 18B have a large intra-level amplitude.

The mean δ¹⁸O_carb values from horses have a more significant variation (range = 3.3 ‰) than those of bovines (range = 2.2 ‰). All individuals show flat δ¹³C_carb intra-tooth profiles (< 0.4 ‰), except for red deer (1 ‰) (Sect. S4). Intra-tooth δ¹⁸O_carb ranges of individuals are around 1 ‰–2 ‰ for horses and 1 ‰–3 ‰ for bovines. Some of the individuals analysed do not show complete annual cycles. No δ¹⁸O_carb sinusoidal profiles are detected in three teeth; the other six have particularly unclear profiles. After modelling, individual δ¹⁸O_carb ranges oscillated between 2.7 ‰ and 7.4 ‰ (Sect. S5). MATs oscillated between 4.6 °C and 12.6 °C (MATA = −8.8/−0.9 °C), with mean summer temperatures around 20.5 °C and mean winter temperatures around −1.1 °C. MAPs extend between 376 and 784 mm (MAPA = −656/−248 mm) (Table 4). Non-important differences in rainfall estimations based on bovines and equids are noticed, probably because they feed on similar ecological resources. Diachronic trends are unclear along the sequence, but mean annual and winter temperatures from levels 18C and 18B seem slightly lower. MAP estimations oscillated more in the upper levels.

4.3 Labeko Koba (Châtelperronian and Aurignacian, 45.1–36.3 ka cal BP)

This work includes bovine (n = 4), horse (n = 4), and red deer (n = 1) teeth from levels related to the Châtelperronian (IXb inf), Proto-Aurignacian (VII), and Aurignacian (VI, V, and IV). Significant differentiation in mean δ¹³C_carb between bovines and horses is observed, with higher values between −9.3 ‰ and −10.4 ‰ in bovines (δ¹³C_diet = −25 ‰ to −23.8 ‰) than in equids, whose values extend from −12 ‰ to −11.6 ‰ (δ¹³C_diet = −25.8 ‰ to −25.2 ‰) (Fig. 3). These horses' values are within the ranges observed from this species in the region. Red deer have similar δ¹³C_carb values to those of horses (δ¹³C_carb = −12.3 ‰; δ¹³C_diet = −25.5 ‰). Mean δ¹⁸O_carb values are similar between species, from −7.2 ‰ to −4.7 ‰ (δ¹⁸O_mw = −8.5 ‰ to −6.1 ‰). However, bovines have a very high variation within mean δ¹⁸O_carb values (2.1 ‰), also reflected in the intra-tooth profiles. These δ¹⁸O values are lower than in other Vasco-Cantabrian sites, especially for two individuals in levels VII and V (Table 3). Differences in δ¹³C_diet values between bovines and horses result in isotopic niche differentiation between both species (Fig. 4). The red deer niche is placed within the horses' niche. The evolution of niche over time cannot be evaluated by levels due to the limited sample. Considering the isotopic compositions by levels (Fig. 5), both bovines and horses experienced a slight increase in mean δ¹³C_diet from level IX inf to IV, from Châtelperronian to Aurignacian. Mean δ¹⁸O_mw values of bovines decrease from VII to V, whereas horses increase from IXb inf to VI to decrease from VI to IV.

Variability in δ¹³C_carb values in intra-tooth profiles is slightly higher (0.1 ‰–0.7 ‰), especially in bovines (0.3 ‰–0.9 ‰), with more oscillating profiles than the generally flat profiles observed in horses and red deer (Sects. S4 and S5). Intra-tooth profile ranges of δ¹⁸O_carb are also larger within bovines (2 ‰–4 ‰) than in horses (1 ‰–2 ‰). Inverse-modelled individual δ¹⁸O_carb ranges oscillated between 5 ‰–8 ‰ and 2 ‰–4 ‰, respectively. Sinusoidal curves are observed in horses and bovines, but bovine profiles are noisier. The red deer has an extensive δ¹⁸O_carb range (6.3 ‰) from a summer peak to an incomplete winter trough. We detect an inverse relation between δ¹³C_carb and δ¹⁸O_carb at some points in these individual profiles. MATs oscillated between 5.2 and 11.4 °C (MATA = $-\mathrm{5.6}/+\mathrm{1.1}$ °C), with summer temperatures from 14.5 to 27.3 °C and winter temperatures from 1.9 to −4.9 °C. MAP extends between 248 and 521 mm, notably drier than nowadays (MAPA = $-\mathrm{798}/-\mathrm{525}$ mm) (Table 4). Lower rainfall levels and higher seasonal amplitudes are recorded along the sequence, especially in samples from the Proto-Aurignacian level VII. Relevant differences are noticed between MAP estimated from bovines and equids, the first providing more arid conditions.

4.4 Aitzbitarte III interior (Gravettian, 27.9 ka cal BP)

A single bovine individual was analysed from Gravettian level V located in the inner part of the cave. It has a high mean δ¹³C_carb (−9.2 ‰) considering the observed range in bovines from the Vasco-Cantabrian region, whereas the δ¹⁸O_carb mean value (−5.5 ‰) is inside the common δ¹⁸O_carb variation observed (Fig. 3). The mean δ¹³C_diet value of −23.8 ‰ is comparable with Canyars and some individuals from Axlor but different from Labeko Koba and El Castillo individuals. The individual δ¹³C_carb fluctuation is slight (0.3 ‰) (Sects. S4 and S5). These teeth show a not quite sinusoidal profile shape in δ¹⁸O_carb, with an intra-tooth range of around 2.2 ‰. Climatic information is extracted but may only be considered cautiously due to the profile shape and the limited sample size. From the inverse modelled mean δ¹⁸O_mw value (−5.4 ‰), we estimate a MAT of 13 °C (MATA = −0.4 °C), with a summer temperature of 19.7 °C and a winter temperature of −2.9 °C. The MAP estimation reached 235 mm (−1127 mm to nowadays) (Table 4).

4.5 El Otero (Magdalenian, ca. 17.3 ka cal BP)

Two equids and three cervids are included from level IV from El Otero, recently re-dated and chronologically related to the Magdalenian (Marín-Arroyo et al., 2018). The mean δ¹³C_carb values are close, between −11.4 ‰ and −11.3 ‰ for red deer (δ¹³C_diet = −24.4 ‰ and −24.6 ‰) and −11.6 ‰ and −11.3 ‰ for horse (δ¹³C_diet = −25.3 ‰ and −25.3 ‰) (Fig. 3). These δ¹³C values for both species are relatively high concerning other studied samples, especially for cervids (around +1 ‰–2 ‰). Both species have higher δ¹⁸O_carb values concerning the common range of variation observed in the Vasco-Cantabria region, between −5 ‰ and −3.9 ‰ for horses and between −5.1 ‰ and −4.4 ‰ for red deer. When values are transformed to δ¹³C_diet and δ¹⁸O_mw, equid and cervid isotopic niches are separated (Fig. 4). All individuals show δ¹³C_carb intra-tooth profiles of low amplitude (< 0.3 ‰) but especially equids, with an intra-tooth variation around 0.1 ‰ (Sects. S4 and S5). Equids and cervids show δ¹⁸O_carb sinusoidal profiles, with intra-tooth ranges between 1.4 ‰ and 2.4 ‰. Climatic estimations are proposed only for equids, providing MAT estimations from 8.8 °C to 12.6 °C (MATA = $-\mathrm{4.9}/-\mathrm{1}$ °C) and MAP between 400 and 456 mm (MAPA = $-\mathrm{755}/-\mathrm{699}$ mm) (Table 4). A high-temperature seasonality can be seen, with summer temperatures between 19.7 °C and 23.8 °C and winter temperatures from −10.4 °C to −3.1 °C.

4.6 Canyars (Aurignacian, 39.7 ka cal BP)

From the archaeological level I at Canyars, corresponding to the Aurignacian, this work includes bovine (n = 2) and equid (n = 3) teeth. The mean δ¹³C_carb values for bovines are between −9 ‰ and −9.3 ‰ (δ¹³C_diet = −23.6 ‰ and −23.8 ‰), and for horses they are between −10 ‰ and −10.7 ‰ (δ¹³C_diet = −23.7 ‰ and −24.4 ‰) (Fig. 3). In this site, the δ¹³C_carb values for horses are notably higher than in the Vasco-Cantabrian region (around +1 ‰–2 ‰) (Table 3). Both species have relatively high δ¹⁸O_carb values, but they fall inside the range of variation observed in the Vasco-Cantabrian region, between −5.5 ‰ and −3.6 ‰ in bovines and between −4.8 ‰ and −4.4 ‰ in horses. Bovine and equid isotopic niches overlap (Fig. 4), but different responses are seen in mean δ¹⁸O_mw values between the two bovines, with one high mean value but close δ¹³C_diet mean values.

All individuals show flat δ¹³C_carb intra-tooth profiles (< 0.3 ‰ variation). Some individuals analysed do not show δ¹⁸O_carb sinusoidal profiles, with intra-tooth profiles moderately flat and ranging from 1.1 ‰ to 1.6 ‰. We detect an inverse relation between δ¹³C_carb and δ¹⁸O_carb in some points of bovine individual isotopic profiles. MATs oscillated between 9.8 °C and 11.9 °C (MATA = $-\mathrm{5.4}/-\mathrm{3.3}$ °C), with summer temperatures from 16.3 °C to 27.5 °C and winter temperatures from −0.5 °C to 1.8 °C (Table 4). MAP extends between 211 and 316 mm (MAPA = $-\mathrm{431}/-\mathrm{326}$ mm). No substantial differences are noticed in the estimations based on bovines and equids because mean δ¹³C diet values differed relatively little.

5 Discussion

5.1 Diet and ecological niches: carbon ratios

Carbon isotopic ratios are valuable indicators for discerning past animal diets, partially influenced by the physiology of the animal. Considering species trends in the studied sites, bovines have generally higher mean δ¹³C_carb values (from −12.4 ‰ to −8.9 ‰) than horses (from −12.6 ‰ to −11.3 ‰), whereas the red deer fall within the horses' range (from −13 ‰ to 11.3 ‰). In the northeastern site of Canyars, bovines also show higher mean δ¹³C_carb values (−9 ‰ to −9.3 ‰) compared to horses (−10.7 ‰ to −10 ‰). These differentiated isotopic ranges for equids and bovines can be potentially linked to feeding behaviour. Still, these species are expected to present different basal δ¹³C_carb driven by their feeding behaviour and distinct physiological characteristics. Bovines, being ruminants, have been suggested in previous studies to exhibit higher δ¹³C_carb values due to increased methane production (Cerling and Harris, 1999; Tejada-Lara et al., 2018). Therefore, transforming δ¹³C_carb to δ¹³C_diet values using species-specific equations is crucial to mitigate the species-specific impact, particularly when comparing ruminants and non-ruminants. Bovines report δ¹³C_diet values between −27.5 ‰ and −23.5 ‰, and horses report values between −26 ‰ and −25 ‰. These carbon compositions are typical of animals feeding on C3 plants (commonly accepted range between −34 ‰ and −23 ‰), as can be expected from high-latitude ecosystems during the Pleistocene (Bocherens, 2003; Cerling and Harris, 1999; Drucker, 2022).

Environmental factors, such as light exposure, water stress, temperature fluctuations, salinity, and atmospheric CO₂ changes, can influence variations in δ¹³C values in a diet primarily based on C3 plants (Bocherens, 2003; Kohn, 2010). Typically, δ¹³C_diet values below −27 ‰ (δ¹³C_carb = −13 ‰) are associated with animals feeding on C3 vegetation found in closed forested environments, whereas δ¹³C_diet values between −27 ‰ and −23 ‰ are linked to C3 open landscapes, which could include grasslands and steppe areas (Bocherens, 2003). The relatively high δ¹³C_diet observed here points to animals predominantly feeding in open environments. The canopy effect, characterized by a depletion in ¹³C isotopes due to dense tree cover, seems unlikely among the analysed samples, since none of the individuals reported δ¹³C_diet below the standard cut-off of −27 ‰ (Drucker et al., 2008; Kohn, 2010; van der Merwe, 1991). Therefore, in general terms, open mosaic landscapes, ranging from light forests to meadows and grasslands, can be inferred for northwestern Iberia. Given the generally higher δ¹³C_diet values reported by bovines, it is likely that they were foraging in more open environments than horses, and they can be regarded as predominantly grazers. In particular, bovines from El Castillo exhibit distinct feeding behaviour compared to other Vasco-Cantabrian sites, as evidenced by their lower δ¹³C_diet values, indicating a potential preference for browsing and feeding in closer environments, possibly in lightly forested areas. Both extinct aurochs (Bos primigenius) and steppe bison (Bison priscus) are usually classified as grass-dominant mixed feeders during the Pleistocene, although it should be noted that modern European bison (Bison bonasus) could include browsing in their diet (Rivals et al., 2022). For aurochs, a browse-dominated mixed feeding behaviour is also frequently described.

The δ¹³C_diet range in equids also indicates feeding in open environments, suggesting a general mixed feeding pattern for the Vasco-Cantabrian region. However, individuals from northeastern Iberia likely graze in more open environments, as evidenced by their notably higher δ¹³C_diet values compared to the Vasco-Cantabrian region (+1 ‰–2 ‰). Evaluating if other factors contribute to lower δ¹³C_diet values in horses is critical. In the case of equids from the Vasco-Cantabrian region, it should be considered that they were pretreated with a combination of NaClO and acetic acid, which could potentially affect the isotopic values. Samples after organic removal pretreatment can potentially show either higher or lower δ¹³C values and higher δ¹⁸O values based on previous experiments (Pellegrini and Snoeck, 2016; Snoeck and Pellegrini, 2015), with δ¹³C values generally varying below 0.3 ‰. Based on the observation that horses in the Vasco-Cantabrian region present lower δ¹³C_carb values compared to bovines but similar mean δ¹⁸O_carb value ranges, the influence of the pretreatment on our samples is deemed to be limited.

Furthermore, the high variability in δ¹⁸O_carb values at El Castillo and Labeko Koba does not correlate with a significant variation in δ¹³C_carb values. Based on dental wear and stable isotope analysis, Middle and Late Pleistocene horses (Equus ferus) were primarily grazers, although some rare cases have been reported as mixed feeders or browsers, such as at Igue des Rameaux and Schöningen (Rivals et al., 2022). Horse populations from northern and eastern Europe were found to be browsers or mixed feeders, while those from the Mediterranean region tended to be grazers (Rivals et al., 2022).

Finally, the few cervids included in this study exhibit δ¹³C_diet values that frequently overlap with horses, indicating a mixed feeding behaviour that varies from more closed environments in El Castillo to more open habitats in El Otero. During the Pleistocene, the red deer (Cervus elaphus) exhibit a flexible, mixed feeding behaviour, consuming leaves, shrubs, forbs, grass, and sedges, similar to their present-day counterparts (Merceron et al., 2021; Rivals et al., 2022). Today, this species inhabits diverse habitats ranging from steppes to closed temperate forests.

5.2 Seasonality, mobility, and water acquisition: oxygen ratios and intra-tooth profiles

Average values of δ¹⁸O_carb in Vasco-Cantabrian individuals extend between −7.2 ‰ and −3.3 ‰ (Table 3). Even if no clear species patterns in δ¹⁸O_carb are observed, in general, bovines present slightly lower δ¹⁸O_carb values from −7.2 ‰ to 4.8 ‰ than other species; horses have a significant variation from −6.6 ‰ to −3.3 ‰ and red deer from −6.8 ‰ to −4.4 ‰. In Canyars, both species have relatively high δ¹⁸O_carb values that fall inside the variation range observed in the Vasco-Cantabrian region, between −5.5 ‰ and −3.6 ‰ in bovines and between −4.8 ‰ and −4.4 ‰ in horses. Each species shows different δ¹⁸O_carb intra-tooth ranges, with bovines between 1 ‰ and 3 ‰, horses mostly around 1.5 %, and red deer from 1 ‰ to 6 ‰ presenting the higher ranges (Table 3; Sect. S4). After applying inverse modelling to correct the dampening effect (Passey et al., 2005b), the majority of teeth increase the δ¹⁸O_carb intra-tooth range, between 3 ‰ and 8 ‰ for bovines and 2 ‰ and 7 ‰ for horses (Sect. S5). Most bovines from Axlor and Labeko Koba and horses from El Castillo and El Otero exhibit well-defined sinusoidal profiles in their δ¹⁸O_carb and large intra-tooth individual ranges, related to the predominant consumption of water sources that reflect seasonal fluctuations between summer and winter. Although not all samples consistently follow this pattern, specific intra-tooth profiles, particularly those from bovines in El Castillo and Canyars, exhibit sharp profiles with narrow ranges (< 1.5 ‰). This phenomenon was previously reported in the region in preliminary studies conducted at the sites of El Castillo (Jones et al., 2019) and in the Magdalenian levels of the El Mirón Cave (Geiling, 2020).

Non-sinusoidal profiles observed in the data can be attributed to various factors, including sampling techniques, preservation issues, and the inherent variability in the original isotopic signal. Factors related to sampling and methods can be connected to (1) the sampling process (e.g. sampling grooves that are too deep or too distant), (2) the imprecision of the mass spectrometer measurements, (3) uncontrolled effects of sample pretreatments, and (4) diagenetic alterations affecting the carbonate fraction. However, it must be noted that technical reasons, whether related to sampling or pretreatment, do not appear to impact the obtained results significantly. Firstly, this study reproduces the same intra-tooth sampling methods that previously yielded reliable results in similar research (e.g. Pederzani et al., 2023, 2021a). Secondly, non-significant alterations in intra-tooth profiles of pretreated horse samples (El Castillo, Labeko Koba, and El Otero) are noticed in comparison to untreated bovid samples (Sect. S5). Some bovid samples show these non-sinusoidal profiles equally. In sites where both species are analysed, no correlation is observed between δ¹⁸O_carb and δ¹³C_carb. In tooth enamel, diagenetic alterations are generally less pronounced than in bone due to its higher mineral content. However, carbonates within tooth enamel can be more susceptible to diagenesis and recrystallization compared to the phosphate fraction, which contains a more extensive reservoir of oxygen and stronger oxygen bonds (Zazzo et al., 2004; Chenery et al., 2012; Bryant et al., 1996). The carbonate content in our samples, ranging from 3.9 % to 8.9 %, is similar to the proportion found in modern tooth enamel, suggesting no immediate indication of diagenetic alteration. Diagenesis can also be evaluated by comparing the isotopic values of the carbonate and phosphate fractions in a sample, as there is a predictable difference between them. However, phosphate fraction measurements were still unavailable in our study, except at Axlor (Pederzani et al., 2023), where good preservation was attested. Additionally, in the case of diagenetic alteration, we would expect specimens from the same archaeological levels to be affected similarly, which is not the case.

Based on these arguments, it is suggested that the non-sinusoidal δ¹⁸O_carb signal observed in some individuals may not be attributed to poor preservation; instead, it likely reflects the original isotopic signature from water input, which appears to be non-seasonal. Several factors can explain why some teeth do not reflect an evident seasonal fluctuation, which could be related to animals' mobility, the isotopic composition of the water sources, and seasonal buffering within those water sources (Pederzani and Britton, 2019). The main factors considered in our study are (1) the high mobility of the animals analysed among ecosystems with different isotopic baselines due to large migrations; (2) the inland–coastal or short altitudinal movements through the region, which lead to the acquisition of water from sources with different isotopic signal; and (3) the acquisition of water from sources with no clear seasonal signal, such as large bodies of water, rivers, groundwaters, and meltwaters. At mid-latitudes, the temperature effect is currently the dominant factor. However, it is crucial to note that past changes in rainfall density (as the “amount effect”; Dansgaard, 1964) cannot be dismissed from having a more significant role then, particularly during glacial and arid periods. These effects, with their potential to mask temperature oscillations, underscore the urgency and importance of our research in understanding and predicting climate patterns. Furthermore, variability between species and within the same species, even within populations living in the same habitat, is also possible. This can be attributed to multiple factors, from minor differences in foraging and drinking behaviour to slight metabolic and physiological variations, including body size, metabolic rate, breathing rate, moisture content of food, and faeces, among others (Hoppe et al., 2004; Kohn, 1996; Magozzi et al., 2019).

Analyses of nitrogen and sulfur stable isotopes on ungulate bone collagen from Axlor, El Castillo, and Labeko Koba (Jones et al., 2018, 2019; Pederzani et al., 2023) have already revealed large variation ranges linked to the existence of several microenvironments in just a few kilometres within the Vasco-Cantabria region. Long migrations and long hunting distances cannot solely explain these diverse values because of the range of species involved and their likely small-scale movements. In our study, the minimal δ¹³C_carb intra-tooth variation within individuals (< 1 ‰) indicates limited seasonal changes in their feeding behaviour that influenced the carbon isotopic composition (Sect. S4). Therefore, considering the diverse topography of the Vasco-Cantabrian, characterized by steep valleys connecting the Cantabrian Mountains with the Atlantic Ocean through rivers over short distances (30–50 km), the availability in the past of a wide range of water sources in small areas seems highly likely. Certain drinking behaviours can influence δ¹⁸O, as animals may acquire water from various sources, with small streams better reflecting seasonal isotopic oscillations than large lakes or evaporating ponds (see synthesis in Pederzani and Britton, 2019). Systematic consumption of highly buffered water sources can significantly attenuate the final recorded signal. Furthermore, rivers in the region frequently contain meltwater from snow during winter and spring and from water springs.

5.3 Regional trends and ecological niches

This study provides valuable insights despite the limited sample size at each archaeological level. It establishes a baseline of isotopic values for northern Iberia, allowing the evaluation of regional trends. In the northwest, in the Vasco-Cantabrian region, the δ¹³C_carb values obtained oscillated between −13 ‰ and −8.9 ‰, and the δ¹⁸O_carb values were between −7.2 ‰ and −3.3 ‰. These values are within the range expected, considering previous regional studies in ungulates (Carvalho et al., 2022; Jones et al., 2019; Lécuyer et al., 2021; Pederzani et al., 2023). Although oxygen variability trends are less precise, the main factor distinguishing the observed changes over time is the variation in carbon isotopic composition among species and regions. The combination of mean δ¹³C_diet and δ¹⁸O_mw values (Figs. 4 and 5) accentuates disparities in ecological niche overlap between horses and bovines, whereas cervids and horses frequently exhibit shared ecological niches. The dissimilarities between bovines and horses could be attributed to shifts in feeding behaviour, which may be accompanied by ecological and environmental changes, either independently or in parallel.

Comparing the entire dataset and across all sites, the consistently lower δ¹³C_diet values in horses compared to bovids throughout time suggest both animals inhabited open landscapes, with bovines exhibiting a grazer preference and horses showing a mixed feeding diet. Only in the Middle to Upper Palaeolithic transition 18B and 18C levels of El Castillo is an exception observed with lower δ¹³C_diet values in bovines, linked to a higher browser input due to a higher habitat in closer environments, such as open forests, similar to those inhabited by the horses. This generates a niche overlap between horses and bovines, most likely reflecting stable conditions that could support both species in similar ecosystems. Contrarily, in the Châtelperronian and early Aurignacian levels from Labeko Koba, a clear differentiation between horses and bovines is observed, mainly in δ¹³C_diet values, highlighting the occupation of different parts of the landscape by both species. This spatially driven niche separation between species could result from resource competition derived from an unstable climatic period, where species needed to specialize to adapt to the changing conditions. Notable changes are also observed in the δ¹⁸O_carb values from Labeko Koba compared to the older El Castillo and Axlor sites, with bovines exhibiting a higher fluctuation range and the lowest values in the region. These trends are consistent with values observed in bone collagen from previous studies in these sites. During the Middle to Upper Palaeolithic transition in the region, by comparing horses and red deer, a decrease in mean δ¹³C (from −20 ‰ to −21 ‰) and δ¹⁵N (from 6 ‰ to 2.5 ‰) values in bone collagen was observed in contrast to stable red deer mean δ¹³C (Fernández-García et al., 2023; Jones et al., 2018, 2019). This decrease was previously interpreted as niche fractionation, derived from an opening landscape, which drove equids into low-quality pastures compared to cervids. Pollen evidence in the region suggests a prevalence of steppe vegetation and low tree cover for the Châtelperronian and Aurignacian (Iriarte-Chiapusso, 2000).

In the same period in Canyars in the northeastern area, higher mean δ¹³C_diet values are observed in both species (between −23.6 ‰ and −24.4 ‰), indicating a preference for more open landscapes by bovines and equids. The indication of open areas could be linked to the arid climatic conditions associated with Heinrich Event 4 (HE4), which coincides with the formation of the studied level. This predominance of open areas coincides with the presence of typical steppe herbivore species, such as Equus hydruntinus and Coelodonta antiquitatis, the microfauna and pollen taxa, and the data offered by the use-wear analysis on ungulate remains identified at the site (Daura et al., 2013; López-García et al., 2022; Rivals et al., 2017).

Aridity is a plausible explanation for the higher niche partitioning observed in Labeko Koba and the higher δ¹³C_diet values found in Canyars for both species during the Aurignacian. The δ¹³C_diet results of bovines from Aitzbitarte III interior during the Gravettian are consistent with the trend observed in Labeko Koba, where previous studies have already suggested this period to be notably arid and cold (Arrizabalaga et al., 2010). Finally, in the Magdalenian level of El Otero, higher δ¹³C_diet values resemble those observed in Canyars. However, this time, carbon values are related to niche partitioning between horses and red deer. In contrast, higher δ¹⁸O_mw values might indicate warmer conditions but are still associated with open landscapes in the Vasco-Cantabrian area.

5.4 Late Pleistocene climatic evolution in northern Iberia

Carbon and oxygen isotopes were used to estimate quantitative parameters related to past temperatures and precipitation. In the case of oxygen isotopic compositions, an evaluation of environmental water composition can be addressed before approaching temperature estimations. When transformed to δ¹⁸O_mw using species-adapted correlations and correcting bias in seawater δ¹⁸O_mw, the summer δ¹⁸O_mw values obtained from the modelled teeth range from −8.9 ‰ to −2.2 ‰, while the winter values range from −17.1 ‰ to −8.9 ‰. These values can tentatively be compared with the current trends observed in the δ¹⁸O_mw range recorded by the IAEA stations (IAEA/WMO, 2020) in Santander (from −3.5 ‰ in summer to −6.6 ‰ in winter) and in Barcelona (from −2.2 ‰ in summer to −6.3 ‰ in winter) and with the OIPC (Bowen, 2022) estimations for studied locations (from −1 ‰ to −9 ‰) (Sect. S2). As observed in the present, Canyars exhibits mean annual δ¹⁸O_mw values around −8.2 ‰, which is lower than the current δ¹⁸O_mw estimated for this location (−5.4 ‰) but higher than the Labeko Koba mean annual δ¹⁸O_mw (−9.5 ‰). This raises the question of whether the baseline δ¹⁸O_mw differences between Canyars and the other sites can be attributed to Mediterranean rather than Atlantic influence, assuming equivalent air circulation patterns and moisture sources experienced in the past to in the present (Araguas-Araguas and Diaz Teijeiro, 2005; García-Alix et al., 2021; Moreno et al., 2021). However, these comparisons must be approached thoughtfully, considering that moisture fluxes and precipitation trends may have varied significantly during the Pleistocene and the Holocene (Dansgaard, 1964; Shackleton, 1987).

As indicated by the climate reconstructed here, temperatures were colder and precipitation levels were notably lower in the Late Pleistocene period in this region than they are nowadays (Table 4; Sect. S2). From 80 to 50 ka BP, in the Mousterian levels of Axlor, temperatures were slightly colder than today, but older levels showed higher differences between summer and winter temperatures. Rainfall estimations exhibit an unusually arid pattern, possibly affected by bovines predominantly feeding in open areas at that time. This aligns with the impact of basal feeding behaviour on rainfall estimations, as previously advised by Lécuyer et al. (2021). In this case, it is not possible to isolate the effect of diet from environmental interference, but previous studies have highlighted stable climatic conditions at the site (Pederzani et al., 2023). Climatic reconstruction, relying on a compilation of lake sediments from northern Iberia (Moreno et al., 2012), suggests that, from late MIS 4 to 60 ka cal BP, cold but relatively humid conditions predominated, with drier conditions emerging later. Additionally, stalagmites from Ejulve cave in the Iberian range indicate a dry climate until 65.5 ka BP, preceding HE6, followed by more humid conditions afterwards (Pérez-Mejías et al., 2019).

During the late Middle Palaeolithic and Early Aurignacian occupations, the observed shift in the niche configuration of species suggests potential climatic perturbations. There is a decreasing trend in temperature from the Transitional Aurignacian levels in El Castillo (18C and 18B; ca. 47–46 ka cal BP) to the Châtelperronian (X inf; 45.1 ka cal BP) and Early Aurignacian (VII–V; from 40.7 to 36.3 ka cal BP) levels in Labeko Koba. Lower mean annual and winter temperatures are particularly notable at El Castillo and Labeko Koba. Labeko Koba levels exhibit high seasonal amplitude, especially at level VII. Additionally, there is a slight decrease in rainfall and increased fluctuations from the Transitional Aurignacian levels in El Castillo (18B–18C) to the Aurignacian levels in Labeko Koba (VII–V). Previous studies in the northern Iberian region underlined an environmental and ecological shift after GS13/HE5, from 48 to 44 ka cal BP, based on a progressive trend to colder temperatures, an increase in aridity, and open environmental conditions, matching with the late Neanderthal occupations, followed by a population hiatus before the arrival of anatomically modern humans (Fernández-García et al., 2023; Vidal-Cordasco et al., 2022). This episode coincides with the maximum extent of glaciers in this region, as recorded in Lake Enol and Vega de Comeya, and a significant decrease in plant biomass and herbivore abundance around 44 to 38 ka BP (Ballesteros et al., 2020; Jiménez-Sánchez et al., 2013; Ruiz-Fernández et al., 2022). Moreover, previous isotopic analyses in the region pointed to some ecological alterations considering perturbations observed in the δ¹³C and δ¹⁵N of bone collagen (Jones et al., 2018, 2019). This tendency towards increased aridity aligns with observations made in regional lake sediments from northern Iberia between 60 and 23.5 ka cal BP, marked by abrupt climate changes associated with HE5 (Moreno et al., 2012). Supporting this, the marine core MD04-2845 in the northern margin of Iberia reveals a decline in the Atlantic forest and an expansion of steppe and cold grasses from 47 to 40 ka BP (Fourcade et al., 2022).

When comparing the environmental reconstruction of the Aurignacian period between the Vasco-Cantabrian (levels V–IV from Labeko Koba) and the northeastern region (level I from Canyars), which are synchronous to HE4 (39 ka BP), this study reveals notably lower rainfall levels for the latter. This is due to the feeding behaviour observed in animals, mainly in open areas. However, these drier conditions align with the specific climatic conditions expected for this period and support previous findings revealing aridity and the predominance of open landscapes (Daura et al., 2013; Rivals et al., 2017). The temperature data indicate that, at Canyars, colder conditions were experienced, especially during the winter season, compared to the present. However, in comparison to Labeko Koba, Canyars experienced warmer conditions. As explained earlier, the Mediterranean Basin had consistently higher temperatures, even during colder periods. This is consistent with the persistence of Mediterranean open forests in the surroundings, as indicated by other studies (López-García et al., 2013; Rivals et al., 2017). Continuous natural records are lacking in the northeastern Iberian margin. However, the inland stalagmite record from Ejulve cave (Pérez-Mejías et al., 2019) and the sedimentary lacustrine sequence of Cañizar de Villarquemado (González-Sampériz et al., 2020) have identified the most arid intervals during HE5 and HE4. These periods were characterized by steppe vegetation expansions, followed by deciduous woodland expansion. To the south, the Padul sequence agrees with cold and dry conditions alternating with forest recovery (Camuera et al., 2019), as documented in the Alborean Sea (Martrat et al., 2004).

Finally, the sites Aitzbitarte III interior (27.9 ka cal BP) and El Otero (17.3 ka cal BP) provided valuable climatic insights into the Vasco-Cantabrian region during the Upper Palaeolithic, specifically during the Gravettian and Magdalenian, respectively. Considering previous research in the region, the climatic trend reported for the Aurignacian, characterized by colder and more arid conditions, was expected to continue or even intensify during the Gravettian (Fernández-García et al., 2023; Garcia-Ibaibarriaga et al., 2019; Lécuyer et al., 2021). Both sites indicate lower precipitation in this area than today, indicating significant aridity, with ungulates feeding predominantly in open landscapes. However, El Otero's higher mean annual temperatures recorded in Magdalenian horses, with respect to other sites within the Vasco-Cantabrian, are consistent with a climatic amelioration following the Last Glacial Maximum (Jones et al., 2021). MIS 2 is marked by the most extreme glacial conditions, as indicated by NGRIP and marine cores in Iberian margins (Martrat et al., 2004; Sánchez Goñi et al., 2002). However, other regional proxies, such as lake sediment and the stalagmite sequence in Pindal Cave (Moreno et al., 2010), suggest a complex and highly variable climate during MIS 2. These proxies identify the coldest and most arid period within MIS 2 as the interval from 18 to 14 ka cal BP rather than the global Last Glacial Maximum (23 to 19 ka cal BP).

6 Conclusions

This study provides a detailed analysis of the temporal evolution of the environmental and climatic conditions in northern Iberia, spanning from the Middle Palaeolithic to the late Upper Palaeolithic (from GS21 to GS2, ranging from 80 ka BP to 17 ka cal BP). In the Vasco-Cantabrian region, the results reveal a heterogeneous open mosaic landscape, ranging from light forest to meadows and grasslands. This landscape reconstruction is primarily inferred by the feeding locations of the studied animals and, consequently, related to the ecosystems where hominins captured them. Despite shifts in niche configuration observed between equids and bovines, both species typically foraged in open areas, with bovines showing a higher preference for grazing. Only in El Castillo, during the late Mousterian and Transitional Aurignacian levels, did bovines show unusually low δ¹³C_diet related to higher browsing and overlapping with horse isotopic niche. This might indicate a slightly closed mosaic landscape that could sustain both species. In contrast, only horses from Canyars exhibit a preference for grazing behaviour.

Stable climatic conditions are described for the Mousterian in Axlor and El Castillo levels from 80 to 50 ka cal BP. However, some elements indicate environmental perturbations initiated during the Transitional Aurignacian levels of El Castillo, around 48–45 ka BP and after HE5/GS13. After GS12 (44.2–43.3 ka BP), horses and bovines potentially occupied different ecological niches during the Châtelperronian and early Aurignacian levels of Labeko Koba, pointing to a species' environmental specialization, which can be a consequence of competition for food resources during an unstable ecological period. The climatic estimations indicate a temperature shift during this period, with a slight decrease in temperature and evidence of fluctuations in rainfall. Previous environmental studies on the region have underlined ecological stress and increasing aridity from around 42.5 ka cal BP, which may relate to a broader ecosystem decline. When comparing the environmental conditions during the Aurignacian period in the northeast (Canyars) and the northwest (Labeko Koba), the first had higher baseline temperatures but also experienced higher aridity. Animals continued to feed on open landscapes during the Gravettian and Magdalenian levels in the Vasco-Cantabrian region, represented by Aitzbitarte III interior and El Otero. However, there is evidence of a temperature recovery after the LGM at the El Otero.

The results presented here, derived from the first extensive sampling in the Vasco-Cantabrian, establish the basis of future stable isotopic studies on faunal tooth enamel in Iberia. Despite the uncertainties inherent in this work, both δ¹⁸O and δ¹³C contributed to the regional climatic characterization, including the estimation of temperatures and precipitations, and to the seasonality range between summer and winter. The potential influence of pretreatment effects and uncontrolled diagenetic alterations on the enamel carbonate fraction has been assessed. However, complementary diagenetic tests, using new techniques like δ¹⁸O_phos and FTIR analyses, are advised in further works to gain more insights into sample preservation. Ongoing sulfur, hydrogen, and strontium studies will provide additional information on the mobility patterns of animals hunted by Late Pleistocene hominins and, therefore, will help us better understand the ecological and environmental context occupied by Neanderthals and modern humans and their landscape use in this particular region. Finally, a more comprehensive characterization of the baseline oxygen values would also enhance the environmental interpretation of the existing data.