# Human forager response to abrupt climate change at 8.2 ka on the Atlantic coast of Europe

A. García-Escárzaga, I. Gutiérrez-Zugasti, A.B. Marín-Arroyo, R. Fernandes, S. Núñez de La Fuente, D. Cuenca-Solana, E. Iriarte, C. Simões, J. Martín-Chivelet, M.R. González-Morales, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

A. García-Escárzaga, I. Gutiérrez-Zugasti, A.B. Marín-Arroyo, R. Fernandes, S. Núñez de La Fuente, et al.. Human forager response to abrupt climate change at 8.2 ka on the Atlantic coast of Europe. Scientific Reports, 2022, 12 (1), pp.6481. 10.1038/s41598-022-10135-w. hal-03662110

HAL Id: hal-03662110

https://hal.science/hal-03662110

Submitted on 30 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

OPEN

# Human forager response to abrupt climate change at 8.2 ka on the Atlantic coast of Europe

Asier García-Escárzaga<fnref n="1" />,²,³✉, Igor Gutiérrez-Zugasti⁴✉, Ana B. Marín-Arroyo⁵,⁶, Ricardo Fernandes<fnref n="1" />,⁷,⁸, Sara Núñez de la Fuente⁴, David Cuenca-Solana⁴,⁹, Eneko Iriarte¹⁰, Carlos Simões¹¹, Javier Martín-Chivelet¹², Manuel R. González-Morales⁴ & Patrick Roberts<fnref n="1" />,¹³

The cooling and drying associated with the so-called '8.2 ka event' have long been hypothesized as having sweeping implications for human societies in the Early Holocene, including some of the last Mesolithic hunter-gatherers in Atlantic Europe. Nevertheless, detailed 'on-site' records with which the impacts of broader climate changes on human-relevant environments can be explored have been lacking. Here, we reconstruct sea surface temperatures (SST) from $\delta^{18}$O values measured on subfossil topshells *Phorcus lineatus* exploited by the Mesolithic human groups that lived at El Mazo cave (N Spain) between 9 and 7.4 ka. Bayesian modelling of 65 radiocarbon dates, in combination with this $\delta^{18}$O data, provide a high-resolution seasonal record of SST, revealing that colder SST during the 8.2 ka event led to changes in the availability of different shellfish species. Intensification in the exploitation of molluscs by humans indicates demographic growth in these Atlantic coastal settings which acted as refugia during this cold event.

Current global climate warming is having, and will continue to have, widespread consequences for humans. Looking to the past, multiple climatic and environmental changes have long been thought to have shaped human evolution and behaviour<fnref n="1" />⁻³. The Holocene (11.7–0 ka cal BP) is a geological epoch characterized by comparative stable climatic conditions. However, that stability was punctuated by a series of sudden climate changes, particularly during the Early Holocene<fnref n="1" />. Among these, the ‘8.2 ka event’ has been identified as the largest and most abrupt climatic event of the Holocene⁵,⁶. Climate scientists suggest that this ‘event’ was the result of an outburst of glacial meltwater from the Laurentide lakes in North America⁷. The influx of cold water into the Atlantic Ocean led to a reduction of sea surface salinity and a decline of the Atlantic meridional overturning circulation (AMOC), provoking a reduction in sea surface temperatures (SST) across the North Atlantic⁸. The cooling effects of this event have been documented in proxies from the Greenland ice cores⁶ and across Europe⁵,⁹,¹⁰. Short and sharp periods of colder or drier conditions have also been recorded ~ 8.2 ka throughout the northern¹¹⁻¹³ and southern hemispheres¹¹,¹²,¹⁴. Nevertheless, while now a well-established climatic phenomenon, the sweeping impacts of the ‘8.2 ka event’ on different environments often remain assumed rather than proven, and local records of marine or terrestrial conditions available at an appropriate resolution often remain lacking. Furthermore, while temperature or rainfall changes associated with the ‘8.2 ka event’ are often well-established, there are few palaeoclimatic proxies available that have enabled insights into the influence of the ‘8.2 ka event’ on the seasonality of climatic conditions in different parts of the world.

<footnote n="1">¹Department of Archaeology, Max Planck Institute for the Science of Human History, Jena, Germany. ²Departamento de Ciencias Humanas, Universidad de La Rioja, Logroño, Spain. ³Department of Prehistory and Institute of Environmental Science and Technology (ICTA), Universitat Autònoma de Barcelona, Bellaterra, Spain. ⁴Instituto Internacional de Investigaciones Prehistóricas de Cantabria, Universidad de Cantabria, Santander, Gobierno de Cantabria, Santander, Spain. ⁵Grupo de I+D+I EVOADAPTA, Departamento de Ciencias Históricas, Universidad de Cantabria, Santander, Spain. ⁶Department of Archaeology, University of Cambridge, Downing Street, Cambridge CB2 3DZ, UK. ⁷School of Archaeology, University of Oxford, Oxford, UK. ⁸Faculty of Arts, Masaryk University, Brno, Czech Republic. ⁹Centre de Recherche en Archéologie, Archeosciences, Histoire (CReAAH), UMR-6566, Université de Rennes 1, Rennes, France. ¹⁰Laboratorio de Evolución Humana-IsoTOPIK Stable Isotope Laboratory, Universidad de Burgos, Burgos, Spain. ¹¹ICArEBH—Interdisciplinary Center for Archaeology and the Evolution of Human Behaviour, Universidade Do Algarve, Faro, Portugal. ¹²Facultad de Ciencias Geológicas, Universidad Complutense de Madrid & Instituto de Geociencias (CSIC-UCM), Madrid, Spain. ¹³School of Social Sciences, University of Queensland, Queensland, Australia. ¹⁴email: garcia-escarzaga@shh.mpg.de; fernandoigor.gutierrez@unican.es</footnote>
These gaps relating to the impact of this cold event on local environments are problematic given that there have been a number of attempts to link the '8.2 ka event' to changes in past human societies, with cooling and drying being argued to have stimulated behavioural or settlement responses¹⁵,¹⁶. In the Middle East, a 300-year period of drying and cooling has been associated with economic and cultural changes leading to an increase in social stratification and urbanism in the region¹⁷. Aridity at ~8.2 ka has also been correlated with changes in settlement patterns and human subsistence in North and East Africa¹⁸. On the Atlantic façade of Europe, the '8.2 ka event' occurred at a time when some of the last Mesolithic hunter-gatherers in Europe were producing dense shell middens along the coast prior to the arrival of agriculture¹⁹. Given the apparent reliance of human societies on marine resources, it has been hypothesised that this abrupt climate change must have affected the availability of exploited marine taxa and thus the foraging behaviours, settlement patterns, and spatial mobility of prehistoric human populations²⁰. However, as with previous examples, high-resolution local palaeoclimatic records that are directly associated with areas of human occupation have remained generally lacking¹⁶. Moreover, given the known importance of seasonal patterns in subsistence strategies to hunter-gatherers²¹,²², seasonal records of changes in marine and terrestrial conditions are particularly important. Nonetheless, suitable proxies and precise sequence chronologies have remained scarce.

Marine molluscs have very specific environmental requirements and are extremely sensitive to climate change, as revealed by previous studies across Atlantic Europe²³,²⁴. Recent research has demonstrated that the current rising of SST is forcing a rapid expansion of the northern limit of warm-adapted marine species²⁴,²⁵. Similarly, warmer conditions during the Early Holocene have been argued to have led to significant changes in species representation when compared to Pleistocene glacial conditions²⁶. Mollusc species recovered from archaeological deposits thus have high potential to provide temperature information of direct significance to humans occupying particular sites and exploiting marine resources²⁷. Recent developments in the radiocarbon dating of shell carbonate have improved our understanding of the dating of this material, allowing researchers to accurately determine the chronology of climate change decoded from shell remains²⁸. Beyond taxonomic indications, stable oxygen isotope ratios derived from archaeologically-recovered marine mollusc shells (δ¹⁸O₃ₘ) have been demonstrated to be powerful recorders of the seasonal SST variations experienced by that animal in the past²⁹,³⁰, offering unparalleled insights into the seasonality of marine environments relied upon by past human populations. Significantly, the effectiveness of these seasonal palaeothermometers has also been demonstrated on the European Atlantic façade for *Phorcus lineatus* (da Costa, 1778)²¹,³¹,³², a topshell species that comprises a significant proportion of Mesolithic shell middens found in the region. The difficult task, until now, however, has been to find a shell midden with sufficient temporal resolution to enable inference of short and abrupt changes in climate conditions and human behaviour¹⁹.

Here, we report chronological, archaeomalacological and isotopic data from shellfish remains recovered from the site of El Mazo in northern Spain (Fig. 1; Supplementary Text 1; Supplementary Fig. 1). We performed a large number of radiocarbon measurements on different materials, including bone from terrestrial animals, microbotanical remains, and marine mollusc shells to generate a high-resolution chronological sequence of the human use of the cave. Meanwhile, measurements of δ¹⁸O₃ₘ on subfossil topshells *P. lineatus* were used to accurately reconstruct short-term changes in coastal SST throughout the Early Holocene. In addition, distributions of shell assemblages were quantified to determine marine mollusc species representation, diachronic changes in shell size, harvesting areas, and the occurrence of shell size selection by humans. Significantly, the site is unique in the context of the European Atlantic façade given its long stratigraphic sequence (9–7.4 ka cal BP), but especially, the high chronological resolution of each unit. As a result, it offers an important opportunity to decipher, for the first time, the duration and the seasonal impact of the 8.2 ka event on SST in Atlantic Europe, its impacts on littoral resources that were being heavily exploited by humans, and its broader ramifications for the settlement patterns and economies of some of the last dedicated hunter-gatherer groups in Europe prior to the expansion of domesticated plants and animals into the region. The study of ‘on-site’ palaeoenvironmental indicators, at a seasonal resolution, also demonstrates the significant potential of this methodology, which can be applied to multiple sites, to explore the overall, and intra-annual, effects of abrupt climate changes on specific environments where the resolution of archaeological sites allows.

## Results

**Bayesian chronological modelling of El Mazo.** Since 2012, a total of 65 ¹⁴C AMS measurements were obtained for archaeological samples from 23 of the 25 stratigraphic units documented in squares X15 and X16 of the shell midden site of El Mazo (Fig. 1c; Supplementary Text 1; Supplementary Table 1; Supplementary Table 2): 18 on bones from terrestrial mammals, specifically ungulates with anthropogenic modifications; 5 on wood charcoal; and 42 on marine shells. 57 ¹⁴C measurements were produced at the Oxford Radiocarbon Accelerator Unit (ORAU) and 8 at the International Chemical Analyses (ICA) laboratory (Supplementary Table 1; Supplementary Table 2). The Bayesian modelled sequence places the shell midden formation in the Early and Middle Holocene, between 9.0 and 7.4 ka cal BP (Fig. 2; Supplementary Fig. 2; Supplementary Data 1). Bayesian modelling offers a high chronological resolution, allowing for the identification of a unit which formed precisely during the 8.2 ka event (Fig. 2). Three out of thirteen shells used for δ¹⁸O₃ₘ analysis to reconstruct SST, and that were radiocarbon dated, were identified as outliers by the Bayesian model (Supplementary Data 1) and were thus discarded from the palaeoclimatic discussion (Supplementary Text 2).

**Coastal sea surface temperatures (SST) during the early Holocene.** SST derived from δ¹⁸O₃ₘ values indicates important changes in absolute values for both annual maxima and minima, as well as in the SST annual range over time (Fig. 3; Supplementary Data 2). The coldest winter temperature (9.4 °C) and largest seasonality (i.e., difference between winter and summer) (13.8 °C) were recorded at the bottom of the sequence
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Location of the study area (Cantabrian region, northern Spain) and the shell midden site of El Mazo. (b) External view of El Mazo; (c) topographic map of the site showing excavation areas and provenance of the studied samples (shaded squares), and (d) stratigraphy of the inner test pit (squares X15 and X16). Satellite map in the upper-left corner of (a) was created by authors using Adobe® Illustrator software. The original base map was extracted from the maps-for-free web site (https://maps-for-free.com/) (© Collaborators of OpenStreetMap, Open Database License [ODbL]). Map in the upper-right corner of (a) was generated by Alejandro Garcia-Moreno using ArcGIS software. Lower map of (a) was created by Luis C. Teira Mayolini using ArcGIS and Adobe Suite software.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Radiocarbon-based chronology for the midden site sequence of El Mazo cave (Asturias, N Spain) modelled using OxCal v.4.4.2<sup>33,34</sup> and the radiocarbon calibration curves IntCal20<sup>35</sup> and Marine20<sup>36</sup>. Sequence of posterior distributions summarize the chronology of archaeological phases at El Mazo. For each archaeological phase, the OxCal Date command for an undated object was used to summarize the chronology of a phase. The full output of the Bayesian chronological model, showing estimates for start and end of archaeological phases, is provided in Supplementary Data 1.</figcaption></figure>

(unit 108: 8995–8760 cal BP). Minimum and maximum SST increased at unit 107 (8550–8330 cal BP) by 2.1 and 1.4 °C, respectively. Seasonality during this period (13.1 °C) was lower than in the oldest unit. For both units dating to the 8.2 ka event (105: 8315–8185 cal BP; 112C: 8180–8070 cal BP) there was a decrease in the minimum temperature in comparison with the three centuries before and, specifically, a reduction of 2.3 °C in summer temperatures. These were the coldest summer temperatures of the entire sequence (22.3–22.6 °C). Seasonality decreased, showing the smallest range of the entire sequence (11.1–11.7 °C). SST then increased in the upper part of the sequence (101B: 7945–7795 cal BP). During this period, the warmest temperature in winter (12.4 °C) was recorded, reporting even slightly higher winter SSTs than today in the region.

**Marine mollusc species representation and climate changes.** From the study of the shell assemblages we identified a minimum number of 80,221 individuals and 23 marine mollusc species (Supplementary Table 3). The limpets *Patella* spp. and the topshells *P. lineatus* together represent ca. 95% of the total Minimum Number of Individuals (MNI) for all stratigraphic units. The Simpson's Diversity Index obtained for all units (ca. 0.7) confirms that only a few species dominate the taxonomic assemblage. The results also showed changes in the abundance of the warm-adapted topshell *P. lineatus* and the cold-adapted limpet *P. vulgata* over time (Fig. 4a), with four major shifts found throughout the shell midden sequence. The most represented species during the formation of the oldest unit of the shell midden (ca. 8825 cal BP) was the limpet *P. vulgata*, while the topshell *P. lineatus* was the most exploited species during the middle part of the sequence (ca. 8785–8300 cal BP). Subsequently, the number of *P. lineatus* topshells showed a notable decrease, while the limpet *P. vulgata* increased, becoming once again the most exploited species (ca. 8255–8020 cal BP). There was a drastic decrease in *P. vul*
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Estimated sea surface temperature (SST) derived from δ¹⁸O₃ₐshell values obtained from *P. lineatus* shells recovered from the shell midden site of El Mazo cave (Asturias, N Spain). SST was calculated using Eqs. (1) and (2). Grey square shows current SST range in N Spain³¹,³².</figcaption></figure>

gata for the most recent units (ca. 8005–7655 ka cal BP), while *P. lineatus* and the limpet *P. depressa* increased in number, the latter becoming the most exploited limpet.

**Shell size variations.** The dimensions of a total of 20,244 *P. lineatus*, 11,032 *P. vulgata* and 4549 *P. depressa* specimens were measured (Supplementary Table 4). The biometric results showed significant diachronic variations in the average shell size per unit for both limpets and topshells (Fig. 4b; Supplementary Fig. 3). The profiles of the three main species reported similar patterns in shell size variations over time. A clear and abrupt decrease in shell size was observed starting at ca. 8300 cal BP, reaching minimum sizes between ca. 8100 and 8070 cal BP. Following this, sizes began to increase again, reaching similar values to those observed prior to the 8.2 ka event between ca. 7935 and 7870 cal BP. An additional abrupt decrease in shell size was observed for the three taxa between the previous period and ca. 7655 cal BP.

Discussion

The abrupt 8.2 ka cold event has been widely described from Greenland and North Atlantic records⁵,⁶,⁸⁻¹¹. However, its impact in shelf coastal areas is scarcely documented, and current chronological resolution of the majority of marine records is insufficient to precisely determine the timing of this event. Bayesian modelling of a total of 65 radiocarbon dates from El Mazo allowed us to achieve the highest-chronological resolution for a Holocene sequence in an archaeological shell-midden located in the North Atlantic (Fig. 2; Supplementary Fig. 2; Supplementary Data 1), providing a unique opportunity for understanding the impact of this abrupt climate change on a regional marine environment, mollusc species diversity and composition, human behaviour, and the resilience of some of the last littoral foragers in this part of the world.

The 8.2 ka event has been dated in the NGRIP ice core to between 8250 and 8090 cal BP⁶ (Fig. 4f). However, recently, two North Atlantic SST cooling pulses have been identified within the event, as a result of a slowed AMOC⁴⁰ during two outbursts of the proglacial Lakes Agassiz and Ojibway. These outbursts have been observed in marine records (at 8320 and 8202 cal BP, respectively)⁷,⁴¹ and Iberian speleothem isotopic records (at 8350–8340 and 8221–8211 cal BP, respectively)⁴². The stratigraphic sequence at El Mazo covers the interval between 9 and 7.4 ka cal BP and, according to the Bayesian model, includes several units formed during the 8.2 ka event. Among them, units 105 (8315–8185 cal BP) and 112C (8180–8070 cal BP) were selected for reconstructing SST in the northern Iberian coast and determining whether the 8.2 ka event had an impact in the region. The results from δ¹⁸O₃ₐshell values (Figs. 3 and 4d), clearly demonstrate that colder SST existed during the 8.2 ka event. Beyond an average annual SST cooling of –1.8 °C with respect to the previous unit, sequential isotopic analysis revealed the seasonal impact of the 8.2 ka event on coastal SST for the first time. The decline in littoral SST is, mostly, observable in the summer season (–2.3 °C), with a smaller winter decrease (–0.6 °C) across a period of two centuries (Fig. 3). These differential impacts are likely linked to differential seasonal insolation variations during the Early Holocene (Supplementary Text 3; Supplementary Fig. 4)⁴³.

The analysis of mollusc species representation throughout the archaeological sequence at El Mazo has provided significant information on the impact of decreasing SST on mollusc populations during the 8.2 ka event. Results show that the warm-adapted topshell *P. lineatus* considerably decreased in abundance between ca. 8255 and 8020 cal BP, while the cold-adapted *P. vulgata* increased, becoming the most exploited species during this time (Fig. 4a). Warm-adapted species did not recover a similar percentage to the pre-8.2 ka event period until ca. 7935 cal BP (Fig. 4a). Therefore, climatic change had an important impact, not only on the mollusc population
<figure><img src="image_7.png" /><figcaption>**Figure 4.** (a) Relative representation of the four main marine mollusc taxa recovered from El Mazo site over time (*P. lineatus*, *P. vulgata*, *P. depressa*, and *P. ulyssiponensis*). (b) Average shell size variations on *P. lineatus*, *P. vulgata* and *P. depressa* through time. (c) Percentage of *P. vulgata* shells collected in the lower intertidal zone over time deduced from L/H ratios³⁷. (d) Maximum and minimum estimated SST from δ¹⁸O₃₈ values for different stratigraphic units through shell midden sequence. (e) Percentage of arboreal pollen over time in N Spain³⁸. (f) δ¹⁸O from the Greenland NGRIP ice core³⁹. Blue bar shows 8.2 k event duration according to the Greenland NGRIP ice core⁶.</figcaption></figure>

numbers of each species, but also on human foraging behaviour, by conditioning the food supplies available for consumption by Mesolithic societies. In contrast to *P. lineatus*, the warm-adapted limpet *P. depressa* showed a gradual increase over time, including during the 8.2 ka event (Fig. 4a). This is a consequence of the species-specific SST tolerance. While *P. depressa* is a warm-adapted species, this limpet has a higher resistance to colder temperatures than the topshell *P. lineatus*²³, explaining why *P. depressa* continued to thrive during the 8.2 ka event.
A decrease in shell size was also observed at El Mazo during this colder period, something that has been variously linked to either changes in SST³⁷ or a higher human pressure on coastal resources⁴⁴, mostly in the form of demographic increase and additional resource pressure⁴⁵. At El Mazo, a clear reduction in shell size for the three most dominant mollusc species can be seen during the 8.2 ka event. Given that each of these species have well-documented differences in temperature tolerance²³,²⁶, it seems that changes in SST cannot be considered responsible for changes in shell size. The biometric profiles obtained from the three species show a significant correlation between them, particularly between *P. depressa* and *P. vulgata* (R² = 0.75; *p* < 0.0001). In addition, changes in shell sizes induced by variations in shore morphology and sea level rise over time are not expected, since this was very similar throughout the study period (Supplementary Text 4). Regarding the impact of the 8.2 ka event on the abundance of molluscs, it has been suggested that changes in upwelling, nutrients, and marine productivity could influence the abundance of marine resources available on the coast²⁰. However, a recent study has suggested that ΔR values, which can be used as a proxy for upwelling, were not significantly different during the formation of stratigraphic units dated around the 8.2 ka event compared to the ΔR value derived from a unit dated at ca. 8.4 ka²⁸. As a result, the reduction of shell sizes in littoral areas at 8.2 ka does not seem to be linked to natural marine processes. It therefore appears that, while SST influenced the abundance of different species at the site, the main cause for the decrease in their size was related to changing anthropogenic pressures.

Diachronic estimates of SST obtained from δ¹⁸O₃₆ measurements (Figs. 3 and 4d) and from the distribution of the main mollusc species (Fig. 4a) reveal, for the first time, colder surficial water masses between ca. 8255 and ca. 7935 cal BP in northern Iberia. We interpret this as a local response to the 8.2 ka event, consistent with the cooling recorded by global proxies (e.g., δ¹⁸O from Greenland ice cores)⁵,⁶. Based on the current chronological information (Fig. 2), SST decrease in northern Iberia probably began before that observed in the Greenland ice cores (8250 cal BP). Warm-adapted and cold-adapted species abruptly decrease and increase, respectively, from ca. 8255 cal BP, although this changing trend probably started earlier (ca. 8355 cal BP) (Fig. 4a; Supplementary Text 5). Pollen records from El Mazo also show the same trend³⁸. Tree pollen decrease also started ca. 8355 cal BP, with more substantial decreases occurring after 8250 cal BP (Fig. 4e; Supplementary Text 5). This suggests that the two SST cooling pulses, which occurred in the North Atlantic as a consequence of the drainage of the proglacial Lakes Agassiz and Ojibway⁷,⁴¹, also had an impact on Iberian coastal locations, in agreement with speleothem isotopic records for this area⁴² (Supplementary Text 5). On the other hand, the impact of the 8.2 ka event on SST also seems to end later than previously observed in Greenland records (8090 cal BP). Both warm- and cold-adapted mollusc species do not recover to their previous abundance until between ca. 7975 and 7935 cal BP (Supplementary Text 5). Considering the usually rapid responses of these taxa to climate change²³⁻²⁵, the local effect of the 8.2 ka event on oceanographic conditions seems to have continued until after 8 ka cal BP in northern Iberia. This chronological timescale agrees with multiple proxies from Atlantic Europe, which reported a longer impact of this climatic anomaly than has been observed in Greenland⁵. It also appears that coastal environments recovered faster than interior settings following the 8.2 ka event. Studies previously conducted on pollen³⁸,⁴⁶,⁴⁷ and speleothems⁴⁸,⁴⁹ from northern Iberia suggest that terrestrial ecosystems only started to return to pre-8.2 ka states between 7780 cal BP and 7655 cal BP (Fig. 4e, Supplementary Text 5).

Reconstruction of SST from δ¹⁸O₃₆ values revealed that the seawater temperature decreased ca. -2 °C during the 8.2 ka event in northern Iberia. A similar SST decrease has been recorded based on changes in the abundance of planktonic foraminifera species in the southern Landes Plateau (Bay of Biscay)⁸ and in North Atlantic marine locations⁵⁰. Mathematical models have also predicted the SST cooling of coastal waters in northern Iberia⁴¹,⁵¹, suggesting connections between northern Iberia and the Nordic oceanic realm as part of a common driving mechanism linked to the North Atlantic Current⁸. The SST variations estimated from δ¹⁸O₃₆ values reported in this study follow a similar pattern to that seen for the vegetation observed in the El Mazo pollen records³⁸, which showed a decrease in the percentage of arboreal pollen during the 8.2 ka event (Fig. 4e). Meanwhile, the proportion of pine forest during this period, favoured by colder climate conditions, and other meso-mediterranean species, also increased as a result of more arid and colder conditions³⁸. Our marine proxy results also agree with additional palaeoclimatic information previously published for northern Iberia from terrestrial proxies⁴⁶⁻⁴⁹. Research previously conducted on pollen remains⁴⁶,⁴⁷ and speleothems⁴⁸,⁴⁹ recovered from sites located along the N Iberia has also revealed more arid and colder conditions during this period, suggesting that the 8.2 ka event had impacts for both marine and terrestrial environments exploited by humans in this Atlantic littoral area.

Previous studies in western Europe have hotly debated the impact of the 8.2 ka event, and its associated environmental changes, on human behaviour¹⁶,⁵². On the Iberian Peninsula, Bicho et al.²⁰ reported changes in subsistence strategies, human settlement and technology in the coastal areas of central Portugal after the 8.2 ka event. In addition, research conducted in the Ebro Valley basin showed that, although not all of the valley was depopulated, several settlements were abandoned in the central part of the valley due to colder and drier climate conditions during the 8.2 ka event. This abandonment has been associated with a migration of Mesolithic populations towards less climate sensitive areas located in other parts of the Ebro Valley⁵³ and/or along the Mediterranean coast⁵⁴,⁵⁵. In northern Iberia, archaeological sites located in the mountainous inland areas seem to have been occupied with lower intensity during the harsher period of the 8.2 ka event. The site of El Espertín does not report robust evidence of occupation from 8.4 to 8.0⁵⁶. Only after 8.0 ka cal BP does a clear signal of human occupation emerge in the mountain areas located to the south of the Cantabrian region⁵⁶,⁵⁷. Given the higher levels of human predation on littoral resources observed during this time at El Mazo, it seems highly likely that the northern Iberian coast, a less climatically sensitive region than the nearby mountain areas⁴⁶ and Ebro river basin⁵⁴,⁵⁵, provided a refuge for humans moving from inland regions. A similar scenario has also been argued for this region during the Last Glacial Maximum⁵⁸ when northern Iberia became a refugium for forager groups. Population migrations from adjacent areas are also in agreement with the growing presence of geometric microliths seen in northern Iberia from ca. 8.0–8.2 ka cal BP⁵⁹,⁶⁰. Increasing human pressure on coastal environments during the 8.2 ka event is suggested by decreasing mollusc size, as well as evidence for the growing
exploitation of lower and more exposed areas during this period (Fig. 4c), which are continually wave-beaten, being therefore more dangerous than higher intertidal zones⁴⁴. Furthermore, recent palaeo-demographic studies based on summed probability distributions of calibrated radiocarbon dates have recorded a demographic increase in northern Iberia and Atlantic Portugal from 8.2 ka ago, coincident with a population decrease in the Ebro valley⁶¹. This process has been related to an increased use of marine resources in Portugal⁶². Thus, it seems that coastal sites, like El Mazo, provided more stable, reliable, and predictable environments for human foragers in Atlantic Iberia, with populations moving towards the coast from the increasingly unpredictable interior.

The benefits of increased human use of coastal areas during this climatic event are clear. However, importantly, growing human pressures, and populations, along the northern Iberian coast during this part of the Early Holocene do not seem to have led to obvious overexploitation of coastal resources, despite evidence for an overall decline in mollusc size. Larger specimens are continuously and systematically harvested, instead of smaller taxa, throughout the whole sequence of occupation (Supplementary Text 6, Supplementary Table 5). In addition, limpet sizes never decrease, at least at the overall population level, to 20 mm or lower (Supplementary Fig. 3), the minimum size currently established by the regulations in order to guarantee the survival of these species over the long-term in the area (Boletín Oficial de Cantabria, BOC-2017-4672, https://boc.cantabria.es). Similarly, scheduled littoral exploitation patterns year-round, which were conceived to maximize the meat yield extracted from molluscs, continued to be practiced despite greater intensification in the use of available resources during the 8.2 ka event²². Our results therefore suggest an ongoing application of local marine ecological knowledge by some of the last foragers in western Europe, despite major climatic and demographic changes, with foragers continuing to apply sustainable resource management strategies.

This study highlights that taxonomic, sclerochronological, geochemical and chronological analyses of molluscs from archaeological sites have the potential to uncover direct insights into the impacts of climatic changes on human-relevant local marine environments. Not only that, but the high seasonal resolution of the climate proxies employed enables the seasonal impacts of the 8.2 ka event to be directly observed. The detailed chronology of El Mazo, as well as these high-resolution proxies, provide a more refined linkage between changing environmental conditions and human behaviour in this region during the 8.2 ka event. Based on comparison with inland records derived from palynological studies and geochemical analyses of speleothems, it seems that coastal regions provided crucial refugia for growing human populations bolstered by migrations from the interior. The resolution provided by this methodology has major implications for other studies seeking to determine the significance of climate change for marine environments, and human adaptations, in regions, and periods of time, where ocean resources were crucial to our species and its hominin relatives⁴⁴,⁶³,⁶⁴. Furthermore, understanding the implications for marine ecology and human behaviour of colder and drier conditions provoked by climatic events such as that seen 8.2 ka ago in North Atlantic Europe⁵,⁹,¹⁰, can provide detailed clues to the magnitude and nature of future climate changes, lags in their manifestation between different regions and environments, and their ultimate outcomes for human societies. This is especially the case given that the AMOC is currently weakening⁶⁵. As the palaeosciences become increasingly recognised by earth scientists, policy makers and climate prediction scientists as an important, detailed source of data for climate predictions beyond anecdotal comparisons, high-resolution methodologies, like that presented here, will become increasingly important.

## Methods

**Radiocarbon dating and Bayesian modelling.** An outlier Bayesian chronological model for El Mazo was implemented using the OxCal v.4.4.2 software³³,³⁴. Radiocarbon calibrations relied on the IntCal20 calibration curve for the northern hemisphere³⁵ and on the marine calibration curve from Heaton et al.³⁶. Corrections for marine radiocarbon reservoir effects (ΔR) in shell samples from northern Iberia followed the recommendations in García-Escárzaga et al.²⁸. According to that study, a different ΔR value should be used depending on the mollusc species: -114 ± 170 for *P. vulgata* and -257 ± 120 for *P. lineatus*. A stratigraphic model was defined by grouping samples from the same archaeological layer within a phase with start and end boundaries in accordance with OxCal terminology. These were set in a sequence following the available stratigraphic information. The OxCal model code is given in Supplementary Code 1. The Bayesian model was run 4–5 times and results compared to assess reproducibility, which was observed to be present for scales of ca. 5–10 years and consequently all reported dates here were rounded to the nearest lustrum. The obtained chronological sequence was compared with the Greenland Oxygen Isotopic record (NGRIP)⁶,³⁹. We used the OxCal Date command to provide a reference for an undated object within a phase. This is also used to provide a rough estimate of the chronology of each phase in this study. Estimates of the start, end, and duration of each phase, plus dates for each object within the stratigraphic units, are given in Supplementary Data 1.

**Stable oxygen isotope analysis on *Phorcus lineatus*.** A total of 27 shells from the topshell *P. lineatus*, which were recovered from five different stratigraphic units (108, 107, 105, 112C and 101B) and located at different points within the shell midden sequence of El Mazo (Fig. 1d), were subjected to stable oxygen isotope (δ¹⁸O) measurements to reconstruct sea surface temperature (SST) diachronic variations. Four specimens were analysed per stratigraphic unit, except for SU8 108 and 112, where ten and five shells were selected, respectively. According to modern baseline isotopic studies previously published for this species in northern Iberia⁵¹,³², specimens smaller than 17 mm in diameter were selected to avoid growth stoppages, which would preclude the reconstruction of annual minima and/or maxima temperatures. All shells were measured using a digital calipper to the nearest 0.01 mm. Prior to sampling the shell samples for oxygen isotope analysis, the samples were first placed in 30 vol% H₂O₂ for 48 h to remove organic matter from the shells following a frequently-applied methodology²²,³¹,³². The cleaned shells were then air-dried at ambient temperature.
The methods applied to extract micro-samples of calcium carbonate from the aragonite layers of shells for δ¹⁸O analyses were similar to those previously employed for this species²²,³¹,³². A first sample was taken from inside the edge of the shell aperture (Supplementary Fig. 5) to avoid breaking the edge when removing the outer layers (i.e. *periostracum* and calcite layer). The outer layers were removed using a Dremmel microdrill with a 2 mm drill bit. The remaining samples were taken sequentially along the whorl, but from the outside of the shell (Supplementary Fig. 5). All calcium carbonate micro-samples were taken manually from the aragonite layer using a dentist microdrill with a 0.5 mm drill bit coupled to a stereoscopic microscope. The sampling spots were separated by 0.3 mm, and the carbonate recovered weighed approximately 150–200 µg. A total of 45 calcium carbonate micro-samples were taken from each shell. The only exception made was for the shells recovered from SU 108, where a total of 30 samples were taken from the shell edge. This give a total of 1065 calcium carbonate micro-samples taken for oxygen stable isotope analyses.

Micro-samples taken on the ten shells recovered from SU 108 (n = 300) and on two shells from the other four stratigraphic unit (MA.107.1, MA.107.2, MA.105.1, MA.105.2, MA.112.1, MA.112.2, MA.101.1, MA.101.2) (n = 360) were analysed using an IRMS Thermo Scientific MAT 253 coupled to Kiel device at the University of Michigan (USA). The instrumental error was better than ± 0.12‰. Micro-samples taken on the rest of the nine shells analysed herein (MA.107.3, MA.107.4, MA.105.3, MA.105.4, MA.112.3, MA.112.4, MA.112.5, MA.101.3, MA.101.4) (n = 405) were analysed using an IRMS Thermo Scientific MAT 253 coupled to Kiel IV device at the Complutense University of Madrid (Spain). In this case, the instrumental error was better than ± 0.03‰. In both cases, the powder samples were dissolved with concentrated phosphoric acid at 70 °C and isotopic ratios were calibrated against NBS-18 (−23.2‰) and NBS-19 international standards (−2.20‰). The results are reported as δ¹⁸O (‰) relative to the Vienna Pee Dee Belemnite (VPDB) reference standard.

**Calculation of past sea surface temperatures (SST) from δ¹⁸O<sub>shell</sub> values.** Following previous investigations conducted on topshell *P. lineatus* species in northern Iberia³¹,³², the water-aragonite fractionation factor obtained by Kim et al.⁶⁶ for synthetic aragonite was used to estimate the SST during the mollusc life span:

$$1000\ln\alpha = 17.88 * (10^3/T) - 31.14 \quad (1)$$

where T corresponds to SST in Kelvin and α is the fractionation between water and aragonite described by the equation:

$$\alpha = 1000 + \frac{\delta^{18}\text{O}_{\text{shell}}(\text{SMOW}\,\text{‰})}{1000} + \frac{\delta^{18}\text{O}_{\text{water}}(\text{SMOW}\,\text{‰})}{1000} \quad (2)$$

To estimate the δ¹⁸O<sub>water</sub> values, we employed the correction proposed by Fairbanks⁶⁷. Fairbanks⁶⁷ proposed a gradual correction of + 0.011‰ for each metre of sea level rise. In this investigation, the current annual δ¹⁸O<sub>water</sub> value average (0.90‰)³¹,³² was used to estimate the δ¹⁸O<sub>water</sub> value for the past. To calculate the sea level height during the formation of each stratigraphic unit, radiocarbon dates obtained for each stratigraphic unit and previously published sea level rise estimates for northern Iberia through Early Holocene⁶⁸ were employed. Nevertheless, another additional correction should be applied in order to accurately estimate δ¹⁸O<sub>water</sub> values during the past, as the correction developed by Fairbanks did not consider the possibility of short-term changes, including the '8.2 ka event'. Warmer conditions from the Early Holocene led to the collapse of the Laurentide ice saddle and the Agassiz-Ojibway palaeolakes, and the posterior drainage of cold freshwater, with a very low δ¹⁸O signal, into the North Atlantic⁷. Indeed, the 8.2 ka event was caused by a depletion of the thermohaline circulation in the North Atlantic derived from this freshwater mass released into the ocean. This process, radiocarbon dated to around 8.33 cal BP⁷, led to a depletion of δ¹⁸O<sub>water</sub> values in the North Atlantic, with a reduction of between -0.2 and -0.4‰ estimated for northern Iberia⁴¹,⁶⁹. Research conducted in our study area applying δ¹⁸O analyses on speleothems showed negative shifts in δ¹⁸O<sub>precipitation</sub> around this time, which has been explained as a local signal of the lower δ¹⁸O of sea surface waters⁴⁸. Therefore, in this study, δ¹⁸O<sub>water</sub> values for those units affected by the meltwater discharge were estimated by subtracting 0.3‰ from the values estimated using the Fairbanks correction⁶⁵.

**Archaeomalacological analyses.** The archaeological shells studied here were recovered from the shell midden site of El Mazo cave (N Spain), from 23 stratigraphic units located in the excavation squares X15 and X16, covering the entire chronological sequence. Archaeological materials were recovered across a surface area of 0.5 m². The total volume of sediment excavated and analysed between all stratigraphic units was 293.5 L and the remains employed in our study were recovered using mesh screen sizes of 1, 2, and 4 mm. Anatomic and taxonomic identification and quantification were performed in order to establish species distribution and their relative abundance in the different stratigraphic units. A total of 262,547 mollusc remains were classified into different fragmentation categories following the methodology proposed by Gutiérrez-Zugasti⁷⁰,⁷¹. The fragmentation categories employed are different for bivalves, spiralled gastropods and non-spiralled gastropods. For bivalves: valves (complete, COMV; fragmented valve, FRAV; fragment with complete hinge, FCG; and anterior/posterior hinge fragments, AHF/PHF; separated into right and left) and fragments. For spiralled gastropods: complete individual (COMI), fragmented individual (FRAI), apical fragment (API), stoma fragment (STOF), umbilical fragment (UMBF) and fragments (FTS). For non-spiralled gastropods: complete individual (COMI), fragmented individual (FRAI), apical fragment (APIF) and fragments (FTS). This procedure, based on non-repetitive elements, is similar to those proposed by other scholars for other taxa in other littoral areas⁷². The minimum number of individuals (MNI) for each taxon was calculated using the corresponding formulae: (i) Bivalves: COMV + FRAV + FCH + (AHF or PHF, whichever is greater). First, the right and left valves are counted separately, and the MNI is the greater of the two totals. (ii) Spiralled gastropods: COMI + FRAI + [APIF
or (STOF + UMBF), whichever is greater]. (iii) Non-spiralled gastropods: COMI + FRAI + APIF. The species nomenclature applied in this study was that proposed by the World Register of Marine Species (WoRMS, http://www.marinespecies.org/). Taxonomic identification was performed using different reference manuals and the comparative collection of the Institute of Prehistory of the University of Cantabria (IIIPC). Limpet remains that were not assigned to any *Patella* species as a consequence of taphonomic processes, which precluded us correctly classifying every APIF and FTS, were denoted as *Patella* sp. Finally, in order to accurately estimate the number of each limpet species and to avoid the artificial effects caused by taphonomic processes in the estimation of *Patella* taxa, the MNI corresponding to *Patella* sp. was assigned to each limpet species according to its percentage representation in each unit.

A biometric analysis was also performed (i) to identify changes in shell size through time, (ii) to determine the collection areas of the limpet *P. vulgata*, and (iii) to establish whether there was any kind of shell size selection during collection. Measurements were taken using digital callipers to the nearest 0.01 mm and followed the protocols established for gastropods<sup>71</sup>. The collection areas were defined using the equations proposed by Bailey and Craighead<sup>37</sup> for *P. vulgata*. The degree of exposure of the coast where the shells were collected was calculated using the following equation: Length (0.142) – Height (0.06) + Width (0.0489) – 5.328. When the results were lower than –0.15, the individuals were assigned to exposed areas, while results higher than that limit were assigned to sheltered shores. To determine the intertidal zone where *P. vulgata* specimens were collected, length/height (L/H) ratios were calculated from each limpet. Then the control value obtained by Gutiérrez-Zugasti<sup>71</sup> for northern Iberia was applied to determine whether the specimens were harvested in the lower or higher intertidal zone. In order to establish the occurrence of size selection during collection, a normality test (Shapiro-Wilk) was applied under the assumption that non-normal distributions reflect some kind of shell size selection<sup>71</sup>. The statistical tests were performed using the software PAST. Stratigraphic units where the number of measurable specimens was lower than 30 were not considered for the biometric study.

## Data availability

All of the data reported in this article is provided in the Supplementary Information, Tables, Figures, and Data files. Bayesian modelled results were included in Fig. 2, as well as in Supplementary Fig. 2 and Supplementary Data 1. All stable oxygen isotope values obtained from *Phorcus lineatus* (da Costa, 1778) mollusc species, as well as the sea surface temperatures (SST) estimated from these isotopic values, are provided in Supplementary Data 2. Mollusc species representation and shell sizes over time are included in Supplementary Table 3 and Supplementary Fig. 3, respectively. All other data supporting the findings and interpretations of this study are available in existing publications and the Supplementary Information provided alongside this manuscript. The archaeological remains studied and sampled for oxygen isotope analysis are curated in the Archaeology collection at the University of Cantabria (Santander, Spain). These materials are to be returned to Archaeological Museum of Asturias (Oviedo, Spain).

Received: 26 November 2021; Accepted: 1 April 2022

Published online: 02 May 2022

## References

1. Carleton, C. & Collard, M. Recent major themes and research areas in the study of human-environmental interaction in prehistory. *Environ. Archaeol.* **25**, 114–130 (2020).

2. deMenocal, P. B. Climate and human evolution. *Science* **331**, 540–542 (2011).

3. Potts, R. Evolution and environmental change in early human prehistory. *Annu. Rev. Anthropol.* **41**, 151–167 (2012).

4. Mayewski, P. A. et al. Holocene climate variability. *Quatern. Res.* **62**, 243–255 (2004).

5. Rohling, E. J. & Pálike, H. Centennial-scale climate cooling with a sudden cold event around 8,200 years ago. *Nature* **434**, 975 (2005).

6. Thomas, E. R. et al. The 8.2ka event from Greenland ice cores. *Quat. Sci. Rev.* **26**, 70–81 (2007).

7. Lewis, C. F. M., Miller, A. A. L., Levac, E., Piper, D. J. W. & Sonnichsen, G. V. Lake agassiz outburst age and routing by labrador current and the 82 cal ka cold event. *Quat. Int.* **260**, 83–97 (2012).

8. Mary, Y. et al. Changes in Holocene meridional circulation and poleward Atlantic flow: The Bay of Biscay as a nodal point. *Clim. Past* **13**, 201–216 (2017).

9. Prasad, S. et al. The 8.2 ka event: Evidence for seasonal differences and the rate of climate change in western Europe. *Glob. Planet. Change* **67**, 218–226 (2009).

10. Seppä, H. et al. Spatial structure of the 8200 cal yr BP event in Northern Europe. *Clim. Past Discuss.* **3**, 165–195 (2007).

11. Alley, R. B. & Ágústsdóttir, A. M. The 8k event: Cause and consequences of a major Holocene abrupt climate change. *Quatern. Sci. Rev.* **24**, 1123–1149 (2005).

12. Morrill, C. & Jacobsen, R. M. How widespread were climate anomalies 8200 years ago? *Geophys. Res. Lett.* **32**, 2 (2005).

13. Dixit, Y., Hodeld, D. A., Sinha, R. & Petrie, C. A. Abrupt weakening of the Indian summer monsoon at 8.2 kyr B.P. *Earth Planet. Sci. Lett.* **391**, 16–23 (2014).

14. Bustamante, M. G. et al. Holocene changes in monsoon precipitation in the Andes of NE Peru based on δ18O speleothem records. *Quatern. Sci. Rev.* **146**, 274–287 (2016).

15. Roffet-Salque, M. et al. Evidence for the impact of the 8.2-kyBP climate event on Near Eastern early farmers. *Proc. Natl. Acad. Sci.* **115**, 8705–8709 (2018).

16. Wicks, K. & Mithen, S. The impact of the abrupt 8.2 ka cold event on the Mesolithic population of western Scotland: A Bayesian chronological analysis using 'activity events' as a population proxy. *J. Archaeol. Sci.* **45**, 240–269 (2014).

17. van der Plicht, J., Akkermans, P. M. M. G., Nieuwenhuys, O., Kaneda, A. & Russell, A. Tell Sabi Abyad, Syria: Radiocarbon chronology, cultural change, and the 8.2 ka event. *Radiocarbon* **53**, 229–243 (2011).

18. Vermeersch, P. M. et al. Early and middle holocene human occupation of the Egyptian Eastern desert: Sodmein cave. *Afr. Archaeol. Rev.* **32**, 465–503 (2015).

19. Gutiérrez-Zugasti, I. et al. Shell midden research in Atlantic Europe: State of the art, research problems and perspectives for the future. *Quatern. Int.* **239**, 70–85 (2011).
20. Bicho, N., Umbelino, C., Detry, C. & Pereira, T. The emergence of Muge Mesolithic shell middens in central Portugal and the 8200 cal yr BP cold event. *J. Island Coast. Archaeol.* 5, 86–104 (2010).

21. Mannino, M. A., Spiro, B. F. & Thomas, K. D. Sampling shells for seasonality: oxygen isotope analysis on shell carbonates of the inter-tidal gastropod *Monodonta lineata* (da Costa) from populations across its modern range and from a Mesolithic site in southern Britain. *J. Archaeol. Sci.* **30**, 667–679 (2003).

22. García-Escárzaga, A. *et al.* Stable oxygen isotope analysis of *Phorcus lineatus* (da Costa, 1778) as a proxy for foraging seasonality during the Mesolithic in northern Iberia. *Archaeol. Anthropol. Sci.* **11**, 5631–5644 (2019).

23. Crisp, D. The effects of the severe winter of 1962–63 on marine life in Britain. *J. Anim. Ecol.* **33**, 165–210 (1964).

24. Mieszkowska, N., Hawkins, S., Burrows, M. & Kendall, M. Long-term changes in the geographic distribution and population structures of *Osilinus lineatus* (Gastropoda: Trochidae) in Britain and Ireland. *J. Mar. Biol. Assoc. U.K.* **87**, 537–545 (2007).

25. Hawkins, S. J. *et al.* Complex interactions in a rapidly changing world: Responses of rocky shore communities to recent climate change. *Clim. Res.* **37**, 123–133 (2008).

26. Gutiérrez-Zugasti, I., Cuenca-Solana D. Biostratigraphy of shells and climate changes in the Cantabrian region (northern Spain) during the Pleistocene-Holocene transition. In: *Archaeomalacology Shells in the Arcaheological Record. British Archaeological Reports International Series 2666* (eds Szabó K, Dupont C, Dimitrijevic V, Gómez-Castélum L, Serrand N). Archaeopress (2014).

27. Thomas, K. D. Molluscs emergent, Part I: Themes and trends in the scientific investigation of mollusc shells as resources for archaeological research. *J. Archaeol. Sci.* **56**, 133–140 (2015).

28. García-Escárzaga, A. *et al.* Bayesian estimates of marine radiocarbon reservoir effect in northern Iberia during Early and Middle Holocene. *Quatern. Geochronol.* **67**, 101232 (2022).

29. Andrus, C. F. T. Shell midden sclerochronology. *Quatern. Sci. Rev.* **30**, 2892–2905 (2011).

30. Wang, T., Surge, D. & Mithen, S. Seasonal temperature variability of the Neoglacial (3300–2500 BP) and Roman Warm Period (2500–1600 BP) reconstructed from oxygen isotope ratios of limpet shells (*Patella vulgata*), Northwest Scotland. *Palaeogeogr. Palaeoclimatol. Palaeoecol.* **317–318**, 104–113 (2012).

31. Gutiérrez-Zugasti, I., García-Escárzaga, A., Martín-Chivelet, J. & González-Morales, M. R. Determination of sea surface temperatures using oxygen isotope ratios from *Phorcus lineatus* (Da Costa, 1778) in northern Spain: Implications for paleoclimate and archaeological studies. *Holocene* **25**, 1002–1014 (2015).

32. García-Escárzaga, A. *et al.* Growth patterns of the topshell *Phorcus lineatus* (da Costa, 1778) in northern Iberia deduced from shell sclerochronology. *Chem. Geol.* **526**, 49–61 (2019).

33. Bronk, R. C. Bayesian analysis of radiocarbon dates. *Radiocarbon* **51**, 337–360 (2009).

34. Bronk, R. C. Dealing with outliers and offsets in radiocarbon dating. *Radiocarbon* **51**, 1023–1045 (2009).

35. Reimer, P. J. *et al.* The IntCal20 northern hemisphere radiocarbon age calibration curve (0–55 cal kBP). *Radiocarbon* **62**, 725–757 (2020).

36. Heaton, T. J. *et al.* Marine20-the marine radiocarbon age calibration curve (0–55,000 cal BP). *Radiocarbon* **62**, 779–820 (2020).

37. Bailey, G. N. & Craighead, A. S. Late Pleistocene and Holocene coastal paleoeconomies: A reconsideration of the molluscan evidence from Northern Spain. *Geoarchaeol. Int.* **18**, 175–204 (2003).

38. Muñez S. *Dinámicas socio-ecológicas, resiliencia y vulnerabilidad en un paisaje atlántico montañoso: la región cantábrica durante el Holoceno*. Unpublished PhD dissertation, Universidad de Cantabria (2018).

39. Rasmussen, S. O. *et al.* A new Greenland ice core chronology for the last glacial termination. *J. Geophys. Res. Atmos.* **111**, D06102 (2006).

40. Ellison, C. R., Chapman, M. R. & Hall, I. R. Surface and deep ocean interactions during the cold climate event 8200 years ago. *Science* **312**, 1929–1932 (2006).

41. LeGrande, A. *et al.* Consistent simulations of multiple proxy responses to an abrupt climate change event. *Proc. Natl. Acad. Sci. U.S.A.* **103**, 837–842 (2006).

42. Domínguez-Villar, D. *et al.* Oxygen isotope precipitation anomaly in the North Atlantic region during the 8.2 ka event. *Geology* **37**, 1095–1098 (2009).

43. Lorenz, S. J., Kim, J.-H., Rimbu, N., Schneider, R. R. & Lohmann, G. Orbitally driven insolation forcing on Holocene climate trends: Evidence from alkenone data and climate modeling. *Paleoceanography* **21**, 2 (2006).

44. Gutiérrez-Zugasti, I. Coastal resource intensification across the Pleistocene-Holocene transition in Northern Spain: Evidence from shell size and age distributions of marine gastropods. *Quatern. Int.* **244**, 54–66 (2011).

45. Marín-Arroyo, A. B. Human response to Holocene warming on the Cantabrian Coast (northern Spain): An unexpected outcome. *Quatern. Sci. Rev.* **81**, 1–11 (2013).

46. Muñoz-Sobrino, C., Ramil-Rego, P., Gómez-Orellana, L. & Díaz Varela, R. A. Palynological data on major Holocene climatic events in NW Iberia. *Boreas* **34**, 381–400 (2005).

47. Moreno, A. *et al.* Revealing the last 13,500 years of environmental history from the multiproxy record of a mountain lake (Lago Enol, northern Iberian Peninsula). *J. Paleolimnol.* **46**, 327–349 (2011).

48. Smith, A. C. *et al.* North Atlantic forcing of moisture delivery to Europe throughout the Holocene. *Sci. Rep.* **6**, 24745 (2016).

49. Rossi, C., Bajo, P., Lozano, R. P. & Hellstrom, J. Younger Dryas to Early Holocene paleoclimate in Cantabria (N Spain): Constraints from speleothem Mg, annual fluorescence banding and stable isotope records. *Quatern. Sci. Rev.* **192**, 71–85 (2018).

50. Hald, M. *et al.* Variations in temperature and extent of Atlantic Water in the northern North Atlantic during the Holocene. *Quatern. Sci. Rev.* **26**, 3423–3440 (2007).

51. Matero, I. S. O., Gregoire, L. J., Ivanovic, R. F., Tindall, J. C. & Haywood, A. M. The 8.2 ka cooling event caused by Laurentide ice saddle collapse. *Earth Planet. Sci. Lett.* **473**, 205–214 (2017).

52. Griffiths, S. & Robinson, E. The 8.2 ka BP Holocene climate change event and human population resilience in northwest Atlantic Europe. *Quatern. Int.* **465**, 251–257 (2018).

53. Alday, A. *et al.* The silence of the layers: Archaeological site visibility in the Pleistocene-Holocene transition at the Ebro Basin. *Quatern. Sci. Rev.* **184**, 85–106 (2018).

54. González-Sampériz, P. *et al.* Patterns of human occupation during the early Holocene in the Central Ebro Basin (NE Spain) in response to the 8.2 ka climatic event. *Quatern. Res.* **71**, 121–132 (2009).

55. García-Martínez de Lagrán, I. *et al.* 8.2 ka BP paleoclimatic event and the Ebro Valley Mesolithic groups: Preliminary data from Artusia rock shelter (Unzué, Navarra, Spain). *Quatern. Int.* **403**, 151–173 (2016).

56. Neira Campos, A., Fuertes Prieto, N. & Herrero, A. D. The Mesolithic with geometrics south of the 'Picos de Europa' (Northern Iberian Peninsula): The main characteristics of the lithic industry and raw material procurement. *Quatern. Int.* **402**, 90–99 (2016).

57. Vidal-Encinas, J. M. & Prada-Marcos, M. E. *Los hombres mesolíticos de la cueva de La Braña-Arintero (Valdelugueros, León)*. Jutan de Castillo y León (2010).

58. Jones, J. R., Marín-Arroyo, A. B., Straus, L. G. & Richards, M. P. Adaptability, resilience and environmental buffering in European Refugia during the Late Pleistocene: Insights from La Riera Cave (Asturias, Cantabria, Spain). *Sci. Rep.* **10**, 1217 (2020).

59. Arias Cabal, P. & Fano Martínez, M. Á. Mesolítico Geométrico o Mesolítico con geométricos? El caso de la región Cantábrica. In *El Mesolítico Geométrico en la Península Ibérica* (eds Utrilla, P. & Montes, L.) (Universidad de Zaragoza, 2009).

60. Fuertes-Prieto N, Risseto J, Gutiérrez-Zugasti I, Cuenca-Solana D, González-Morales MR. New perspectives on Mesolithic technology in northern Iberia: data from El Mazo shell midden site (Asturias, Spain). In: *Foraging Assemblages: Papers Presented at the Ninth International Conference on the Mesolithic in Europe, Belgrade 2015* (eds Boric D, Antonovic D, Mihailovic B) (2020).
61. Fernández-López de Pablo, J. et al. Palaeodemographic modelling supports a population bottleneck during the Pleistocene-Holocene transition in Iberia. *Nat. Commun.* **10**, 1872 (2019).

62. McLaughlin, T. R., Gómez-Puche, M., Cascalheira, J., Bicho, N. & Fernández-López de Pablo, J. Late glacial and early Holocene human demographic responses to climatic and environmental change in Atlantic Iberia. *Philos. Trans. R. Soc. B* **376**, 20190724 (2020).

63. Crowther, A. et al. Coastal subsistence, maritime trade, and the colonization of small offshore islands in eastern African prehistory. *J. Island Coast. Archaeol.* **11**, 211–237 (2016).

64. King, C. L. et al. Marine resource reliance in the human populations of the Atacama Desert, northern Chile—A view from prehistory. *Quatern. Sci. Rev.* **182**, 163–174 (2018).

65. Caesar, L., Rahmstorf, S., Robinson, A., Feulner, G. & Saba, V. Observed fingerprint of a weakening Atlantic Ocean overturning circulation. *Nature* **556**, 191 (2018).

66. Kim, S. T., O'Neil, J. R., Hillaire-Marcel, C. & Mucci, A. Oxygen isotope fractionation between synthetic aragonite and water: Influence of temperature and Mg²⁺ concentration. *Geochim. Cosmochim. Acta* **71**, 4704–4715 (2007).

67. Fairbanks, R. G. A 17.000-year glacio-eustatic sea level record: influence of glacial melting rates on the Younger Dryas event and deep-ocean circulation. *Nature* **342**, 637–642 (1989).

68. Leorri, E., Cearreta, A. & Milne, G. Field observations and modelling of Holocene sea-level changes in the southern Bay of Biscay: Implication for understanding current rates of relative sea-level change and vertical land motion along the Atlantic coast of SW Europe. *Quatern. Sci. Rev.* **42**, 59–73 (2012).

69. Hoffman, J. S. et al. Linking the 8.2 ka event and its freshwater forcing in the Labrador Sea. *Geophys. Res. Lett.* **39**, 2 (2012).

70. Gutiérrez-Zugasti, I. Shell fragmentation as a tool for quantification and identification of taphonomic processes in archaeomagnetic analysis: The case of the Cantabrian Region (Northern Spain). *Archaeometry* **53**, 614–630 (2011).

71. Gutiérrez, Z. I. *La explotación de moluscos y otros recursos litorales en la región cantábrica durante el Pleistoceno final y el Holoceno inicial* (Publican, 2009).

72. Harris, M., Weisler, M. & Faulkner, P. A refined protocol for calculating MNI in archaeological molluscan shell assemblages: A Marshall Islands case study. *J. Archaeol. Sci.* **57**, 168–179 (2015).

## Acknowledgements

This research was performed as part of the projects HAR2016-75605-R and HAR2017-86262-P, funded by the Spanish Ministry of Economy and Competitiveness, MINECO. During the development of this research AGE was funded by the University of Cantabria through a predoctoral grant (no code available), by the Basque Country Postdoctoral Programme (Grant number POS_2020_2_0032), by the University of La Rioja through a posdoctoral grant (no code available) and he is currently supported by the Catalonia Postdoctoral Programme through a Beatriu de Pinós fellowship (grant number 2020 BP 00240). ABMA was funded by the SUBSILIENCE ERC-CoG project (Grant agreement No. 818299). This study has also been supported by the Prehistoric Research Consolidated Group of the Basque Country University (IT1223-19), funded by the Basque Country Government. We thank the Max Planck Institute for Science of Human History (MPI-SHH), University of Cantabria (UC), Instituto Internacional de Investigaciones Prehistóricas de Cantabria (IIIPC), Grupo de Ingeniería Fotónica (GIF), University of the Basque Country (UPV/EHU), University of La Rioja (UR) and Autonomous University of Barcelona (UAB) for providing support. We are also grateful to Emma Loftus (McDonald Institute-University of Cambridge) for her initial advice and help with the Bayesian modelling. We would also like to thank José Miguel López-Higuera (GIF), Adolfo Cobo (GIF), Lucía Agudo-Pérez (EVOADAPTA), Roberto Suarez-Revilla (IIIPC), Alvaro Arrizabalaga (UPV/EHU), Miguel Ángel Fano Martínez (UR), Alejandro García Moreno (IIIPC) and Luis C. Teira Mayolini (IIIPC) for their scientific support.

## Author contributions

A.G.E., I.G.Z. and M.R.G.M. designed the isotopic and archaeomalacological project. A.G.E., I.G.Z., D.C.S. and M.R.G.M. excavated El Mazo cave. A.G.E., A.B.M.A. and R.F. carried out the Bayesian modelling. S.N.F. performed palynological investigation. A.G.E., I.G.Z., A.B.M.A., R.F., S.N.F., D.C.S., E.I., C.S., J.M.C. and M.R.G.M. performed the data analysis. All authors interpreted the data. All authors wrote and provided comments on the manuscript.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-022-10135-w.

**Correspondence** and requests for materials should be addressed to A.G.-E. or I.G.-Z.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the
Creative Commons licence, and indicate if changes were made. The images or other third party material in this
article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the
material. If material is not included in the article's Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022