Christophe Lécuyer
email: christophe.lecuyer@univ-lyon1.fr

Claude Hillaire-Marcel

Ariane Burke

Marie-Anne Julien

Jean-François Hélie

Jean- François Hélie

Temperature and precipitation regime in LGM human refugia of southwestern Europe inferred from δ13C and δ18O of large mammal remains

Keywords: Last Glacial, temperature, precipitation, stable isotopes, apatite, mammals, Western Europe

à 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.

Introduction

The Last Glacial, also named Weichselian in Northwest Europe, covers a time period from 115 ka to 11.7 ka. It was characterized by rapid, large amplitude climatic fluctuations recorded in the Greenland ice cores [START_REF] Dansgaard | Evidence for general instability of past climate from a 250-kar ice-core record[END_REF][START_REF] Grootes | Comparison of oxygen isotope records from the GISP2 and GRIP Greenland ice cores[END_REF][START_REF] Ngrip Members | High-resolution record of Northern Hemisphere climate extending into the last interglacial period[END_REF][START_REF] Rasmussen | A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy[END_REF] and North Atlantic marine sediments [START_REF] Heinrich | Origin and consequences of cyclic ice rafting in the Northeast Atlantic ocean during the past 130,000 years[END_REF][START_REF] Bond | Evidence for massive discharges of ice-bergs into the North Atlantic Ocean during the last glacial period[END_REF]Sanchez Goñi et al., 2008). During this period, global climate conditions oscillated between stadial and interstadial conditions, sometimes punctuated by abrupt climate events known as Heinrich events. These highly contrasted climate regimes influenced faunal communities, including humans, through the development of new migration pathways, allopatric speciation in refugia, and enhanced rates of extinction during the coldest and driest stages of the Last Glacial.

The Last Glacial Maximum (LGM), conventionally dated between 23 and 19 ka ago [START_REF] Margo Project Members | Constraints on the magnitude and patterns of ocean cooling at the Last Glacial Maximum[END_REF][START_REF] Mix | Environmental processes of the ice age: land, oceans, glaciers (EPILOG)[END_REF], is defined as the most recent period of maximum global ice volume [START_REF] Mix | Environmental processes of the ice age: land, oceans, glaciers (EPILOG)[END_REF]. It featured lowered greenhouse gas (CH4, CO2) concentrations, a mean sea level ~120 m lower than today and generally cold and dry climatic conditions resulting in significant modifications of both floral and faunal assemblages [START_REF] Lister | The impact of Quaternary Ice Ages on mammalian evolution[END_REF][START_REF] Sommer | Glacial refugia of mammals in Europe: evidence from fossil records[END_REF][START_REF] Clark | The last glacial maximum[END_REF][START_REF] Binney | Vegetation of Eurasia from the last glacial maximum to present: Key biogeographic patterns[END_REF][START_REF] Cao | Attribution of the Last Glacial Maximum climate formation[END_REF].

During this interval, and in response to global tropospheric cooling, the Northern Hemisphere experienced a southward expansion of continental ice sheets, which reached as far as 52°N in western Europe. Various models suggest a decrease in global mean surface air temperatures during the LGM ranging from 4°C to 10°C vs pre-industrial (PI) climatic conditions [START_REF] Masson-Delmotte | Past and future polar amplification of climate change: climate model intercomparisons and ice-core constraints[END_REF][START_REF] Schneider Von Deimling | How cold was the last glacial maximum[END_REF][START_REF] Annan | A new global reconstruction of temperature changes at the Last Glacial Maximum[END_REF].

The main causes of this decrease are thought to be the combined effect of lower atmospheric CO2 values (with a partial pressure as low as 190 ppm; [START_REF] Bouttes | Last Glacial Maximum CO2 and δ 13 C successfully reconciled[END_REF] and ice-sheet forcing, which would jointly account for about 75% of the total LGM cooling (Schneider von [START_REF] Schneider Von Deimling | How cold was the last glacial maximum[END_REF]. Additional cooling effects are related to changes in atmospheric dust content and vegetation cover. In response to this maximal extent of continental ice (glaciers and ice caps), mean global sea level dropped by 120±10 m [START_REF] Lambeck | Sea level and global ice volumes from the Last Glacial Maximum to the Holocene[END_REF].

Faunal associations described in the archaeological sites of southwestern Europe are generally compatible with cold and dry, open environments, such as tundra and dry steppe, but there are important regional differences as suggested by [START_REF] Yravedra | Essai sur la biodiversité des associations de grands mammifères à la fin du Pléistocène dans le Sud-Ouest de l'Europe[END_REF] who defined several biogeographic zones in southwest Europe. During the LGM, the vegetation in Europe was characterized by an expansion of tundra and steppe. In Southwestern France and in the northern part of the Iberian Peninsula, cold temperate fauna (horse, reindeer, ibex) dominated vertebrate assemblages during the LGM, but co-existed with cold-adapted fauna such as mammoth, woolly rhinoceros and reindeer [START_REF] Yravedra | Essai sur la biodiversité des associations de grands mammifères à la fin du Pléistocène dans le Sud-Ouest de l'Europe[END_REF][START_REF] Sommer | Glacial refugia of mammals in Europe: evidence from fossil records[END_REF]Álvarez-Lao and García, 2011;Álvarez-Lao and García, 2011). Recent models suggest that the most suitable habitats for humans during the LGM were concentrated in southwest France and the Iberian Peninsula, excluding the interiors, as also demonstrated by the distribution of archaeological sites (e.g., [START_REF] Ludwig | The Last Glacial Maximum and Heinrich event I on the Iberian Peninsula: A regional climate modelling study for understanding human settlement patterns[END_REF][START_REF] Burke | Risky business: The impact of climate and climate variability on human population dynamics in Western Europe during the Last Glacial Maximum[END_REF][START_REF] Wren | Habitat suitability and the genetic structure of human populations during the Last Glacial Maximum (LGM) in Western Europe[END_REF]. Modern humans who had arrived in Europe ≈ 45,000 years ago survived the LGM in the SW Europe refugia [START_REF] Fu | The genetic history of ice age Europe[END_REF]. A post-ice age expansion followed, with their descendants spreading out in Iberia, after ≈ 19,000 years ago, then defining the Magdalenian culture [START_REF] Fu | The genetic history of ice age Europe[END_REF].

The LGM in Southwest Europe illustrates a major faunal and vegetation partitioning in biogeographic provinces, which constrained the human strategy of land occupation. Whether Mean Annual Precipitation (MAP) or Mean Annual Air Temperature (MAAT) were limiting factors in this strategy remains to be fully elucidated. In thus study, special attention is paid to sites from coastal areas under the influence of the Atlantic Ocean or the Mediterranean Sea.

They represent three distinct biogeographic zones, following the classification of [START_REF] Yravedra | Essai sur la biodiversité des associations de grands mammifères à la fin du Pléistocène dans le Sud-Ouest de l'Europe[END_REF] for southwestern Europe: southwest France (zone VI), Cantabria/Asturias (zone IV) and the Mediterranean coast (zone IIIa).

MAAT and MAP may be reconstructed using oxygen and carbon isotope compositions of the apatite carbonate from vertebrate skeletons. As highlighted by [START_REF] Pederzani | Oxygen isotopes in bioarchaeology: Principles and applications, challenges and opportunities[END_REF], in particular the oxygen isotope analyses of skeletal remains (δ 18 O) are a powerful tool for exploring past human-environment interactions as well as to reconstruct paleoclimatic parameters. An approach not yet used for documenting the LGM refugia of SW Europe will fill this gap using large mammal remains (teeth) from the above defined human occupation areas. We also intend to demonstrate that carbon isotope (δ 13 C) from such remains may inform specifically on paleoprecipitation.

The basic principles to interpret δ 13 C as MAP and δ 18 O as MAAT are summarized as follows:

The most common strategy consists in estimating air temperature from the oxygen isotope composition of the carbonate or the phosphate groups of the tooth enamel of the vertebrate apatite. During tooth growth, the oxygen isotope composition of apatite carbonate (δ 18 O c ) records body temperature (a near constant close to 37°C) and body water composition (δ 18 O bw ) [START_REF] Koch | Oxygen isotope variation in the tusks of extinct proboscideans: a measure of season of death and seasonality[END_REF][START_REF] Bryant | Oxygen isotope partitioning between phosphate and carbonate in mammalian apatite[END_REF][START_REF] Chenery | The oxygen isotope relationship between the phosphate and structural carbonate fractions of human bioapatite[END_REF]. The oxygen isotope composition of body water reflects that of meteoric water [START_REF] Bryant | A model of oxygen isotope fractionation in body water of large mammals[END_REF][START_REF] Kohn | Predicting animal δ 18 O: accounting for diet and physiological adaptation[END_REF][START_REF] Langlois | Box-modeling of bone and tooth phosphate oxygen isotope compositions as a function of environmental and physiological parameters[END_REF] and depends principally on the composition of oceanic moisture source(s), the length and precipitation history of humid air mass trajectories, evapotranspirative processes along these trajectories, air temperature (MAAT) during precipitation, the relative amount of precipitation and altitude relative to sea level [START_REF] Dansgaard | Stable isotopes in precipitation[END_REF][START_REF] Rozanski | Relation between long-term trends of oxygen-18 isotope composition of precipitation and climate[END_REF][START_REF] Fricke | The correlation between 18 O/ 16 O ratios of meteoric water and surface temperature: its use in investigating terrestrial climate change over geologic time[END_REF]. At mid-and high-latitudes (tropical bands excluded), MAAT of areas close to oceanic moisture sources are linearly related to the weighted mean annual δ 18 O mw [START_REF] Dansgaard | Stable isotopes in precipitation[END_REF][START_REF] Yurtsever | Worldwide survey of stable isotopes in precipitation[END_REF][START_REF] Grafenstein Von | A 200-year mid-European air temperature record preserved in lake sediments: An extension of the δ 18 Op-air temperature relation into the past[END_REF][START_REF] Fricke | The correlation between 18 O/ 16 O ratios of meteoric water and surface temperature: its use in investigating terrestrial climate change over geologic time[END_REF][START_REF] Skrzypek | How cold was it for Neanderthals moving to Central Europe during warm phases of the last glaciation?[END_REF][START_REF] Lécuyer | Water on Earth[END_REF]. This property of the water cycle has already been used for the reconstruction of past air temperatures from the oxygen isotope composition of vertebrate tooth enamel from large mammals (e.g., [START_REF] Ayliffe | The preservation of glacial-interglacial climatic signatures in the oxygen isotopes of elephant skeletal phosphate[END_REF][START_REF] Bryant | Oxygen isotopic composition of fossil horse tooth phosphate as a record of continental paleoclimate[END_REF]Tütken et al., 2007;[START_REF] Bernard | Pleistocene seasonal temperature variations recorded in the δ 18 O of Bison priscus teeth[END_REF][START_REF] Fabre | Late Pleistocene climatic change in the French Jura (Gigny) recorded in the δ 18 O of phosphate from ungulate tooth enamel[END_REF].

The reconstruction of MAP from the carbon isotope composition of apatite carbonate was made in several steps. In the first step, the δ 13 C values of plants (δ 13 Cpl) can be estimated from tooth enamel carbonate (δ 13 Ccarb) taking into account the isotopic fractionation that takes place between diet-derived carbon ingested by the herbivores and the carbon mineralized as carbonate present in the apatite crystal lattice of their teeth [START_REF] Tejada-Lara | Body mass predicts isotope enrichment in herbivorous mammals[END_REF]. In a second step, MAP values are calculated using δ 13 Cpl and a transfer function based on plant leaf discrimination factors [START_REF] Kohn | Carbon isotope compositions of terrestrial C3 plants as indicators of (paleo) ecology and (paleo) climate[END_REF] that may be estimated in the case of animal diets strictly based on Calvin Cycle (C3) plants. Atmospheric pCO2-changes and reduced biomass production by plants during glacials [START_REF] Gerhart | Plant responses to low [CO 2 ] of the past[END_REF] may have led to some global relative increase in the abundance of C4 plants, thus to some possible dietary changes, but C3/C4 plant ratios within given biomes remained broadly stable if biome distributions differed (Ehleringer, 2005). Moreover, C4 plants amount to only about 1% in the European vegetation (Collins and Jones 1985). They are slightly more abundant in herb-dominated biomes (see their taxonomic distribution in [START_REF] Pyankov | European plants with C4 photosynthesis: geographical and taxonomic distribution and relations to climate parameters[END_REF]. Now, if such biomes experienced a ~ 10% relative increase over other biomes during the LGM [START_REF] Binney | Vegetation of Eurasia from the last glacial maximum to present: Key biogeographic patterns[END_REF], their content in C4 plants is unlikely to have increased as these plants are strongly sensitive to temperatures [START_REF] Pyankov | European plants with C4 photosynthesis: geographical and taxonomic distribution and relations to climate parameters[END_REF], adding to the equivocal response of the C4/C3 plant ratio of Poaceae, one of their major components, to pCO2 changes [START_REF] Wand | Responses of wild C4 and C3 grass (Poaceae) species to elevated atmospheric CO2 concentration: a meta-analytic test of current theories and perceptions[END_REF]. One may thus assume a nearly exclusive C3-plant vegetation during the LGM interval whatever the changes in biome distribution.

The goal of the present study is thus to build paleoclimate reconstructions for several archaeological sites of the LGM (Figure 1) using the above isotopic approaches. The specific sampled sites are Laugerie-Haute (Dordogne, France), El Cierro (Asturias Province), Hornos de la Peña, El Pendo and El Ruso (Cantabria Province), Nerja (Málaga Province) and Malladetes (Valencia Province). The present study provides the first quantitative information about LGM-MAAT and -MAP ranges of these sites. They will be discussed in the light of proxies based on pollen assemblages (continental climate) and foraminiferal-δ 18 O data (marine moisture sources), climatic model outputs, past topography, as well as in that of reconstructed vegetation cover and faunal communities. Century the site's stratigraphic integrity has come under close scrutiny [START_REF] Texier | Histoire ge ólogique de sites pre historiques classiques du Pe ŕigord : une vision actualise é : la Micoque, la grotte Vaufrey, le Pech de l'Aze ́ I et II, la Ferrassie, l'abri Castanet, le Flageolet[END_REF][START_REF] Delpech | Biostratigraphie des niveaux solutréens de Laugerie-Haute (Les Eyzies, Dordogne, France). Implications archéologiques[END_REF][START_REF] Verpoorte | Improving the chronological framework for Laugerie-Haute ouest (dordogne, France)[END_REF]. The occupation sequence of the site spans the Aurignacian and the Solutrean. For this study, level 4-8, attributed to Magdalenian II, level 10-16 (F. Bordes collection) attributed to Magdalenian I and level 18-20, attributed to Magdalenian 0, were sampled. Biostratigraphically, levels attributed to Magdalenian I-III form a unit, with faunal assemblages dominated by reindeer, while the Magdalenian 0 level, which is dominated by horse, is distinct [START_REF] Delpech | Biostratigraphie des niveaux solutréens de Laugerie-Haute (Les Eyzies, Dordogne, France). Implications archéologiques[END_REF]. The Magdalenian levels are bracketed by ages of ~ 23 and ~21 14 C cal. kyr [START_REF] Delpech | Biostratigraphie des niveaux solutréens de Laugerie-Haute (Les Eyzies, Dordogne, France). Implications archéologiques[END_REF] and can thus be assigned to the LGM. Four sites were sampled in northern Spain, including: El Ruso, Hornos de la Peña and El Pendo in Cantabria, and El Cierro in Asturias (Fig. 1). El Ruso is a small cave located about 6 km from the modern coastline, near El Pendo and El Juyo (unsampled). The cave contains a sequence of occupation that spans from the Mousterian to the Solutrean Level 3, the Upper Solutrean level sampled here overlies a Solutrean level (IV) dated 20-19 14 C cal. kyr (Yravedra, 2010a).

At El Pendo level 4, the samples are labelled "level 4" and attributed to the Upper Solutrean. If this corresponds to level IV (Eschegaray excavations) this level is attributed to the Late Aurignacian, although it is intercalated between Gravettian and Magdalenian levels [START_REF] Hoyos Gómez | Nuevas aportaciones sobre la estratigrafía y sedimentología de los depósitos del Paleolítico Superior de la Cueva de El Pendo (Santander): sus implicaciones[END_REF]. Unfortunately, the stratigraphic integrity of this site is in some doubt (de los [START_REF] De Los Terreros | Las estrategias de subsistencia en la región central de la cornisa cantábrica ¿Continuidad o ruptura? Nivel cero: revista del grupo arqueológico[END_REF]Barquin, 2014). Consequently, results from this site must be interpreted with caution.

At Hornos de la Peña, we sampled level B. This level is dated 23-21 14 C cal. kyr and is attributed to the Aurignacian in its lower part, and the Solutrean in its upper part [START_REF] Breuil | Les premiers travaux de l'Institut de Paléontologie Humaine[END_REF]Obermeier, 1912 cited in Tejero et al., 2008;Yravedra, 2010b). The faunal assemblage at Hornos de la Peña is dominated by horse [START_REF] Tejero | Arte mueble en el Auriñaciense cantábrico. Nuevas aportaciones a la contextualización del frontal grabado de la Cueva de Hornos de la Peña (San Felices de Buelna, Cantabria)[END_REF]. At El Cierro (Asturias), another deeply stratified cave site excavated by several different teams over the years, a revised chronostratigraphy and several AMS dates are available (Alvarez-Fernandez et al., 2016). We sampled levels G (Magdalenian), G1 and H (Solutrean) which form part of the earliest occupation sequence of the cave. AMS dates obtained for these levels indicate that the site was occupied during the LGM ( 14 C cal. ages of 20-19 kyr) (Alvarez-Fernandez et al., 2016).

Two sites on the Mediterranean seaboard of Spain were sampled: Nerja and Malladetes (aka Mallaetes) (Fig. 1). Nerja, the southernmost sample site, is a deeply stratified cave located on the coastal plain of eastern Andalusia. The site contains one of the longest occupation sequences in Mediterranean Spain, spanning the Gravettian to the Copper Age [START_REF] Tortosa | Cueva de Nerja (Málaga, Espana) y los inicios del solutrense en Andalucia[END_REF][START_REF] Pardo | 70 fechas para una cueva. revisión crítica de 70 dataciones c14 del Pleistoceno superior y Holoceno de la cueva de Nerja (Málaga, Andalucía, España). Espacio, Tiempo[END_REF]. Paleoclimate proxies and a large number of recently reviewed radiocarbon dates provide a robust chronostratigraphic frame for the sequence (Jordá [START_REF] Pardo | 70 fechas para una cueva. revisión crítica de 70 dataciones c14 del Pleistoceno superior y Holoceno de la cueva de Nerja (Málaga, Andalucía, España). Espacio, Tiempo[END_REF]. Unit 8 (vestibule) was sampled for this study. This unit contains Solutrean material and is generally considered to date to the

LGM, based on paleoclimate proxies and conventional and AMS radiocarbon dates yielding 14 C cal. ages comprised between 23 and 21 kyr [START_REF] Vallejo | La fauna de mamíferos del Solutrense en la Cueva de Nerja[END_REF][START_REF] Tortosa | Cueva de Nerja (Málaga, Espana) y los inicios del solutrense en Andalucia[END_REF], although a much younger and possibly unrepresentative AMS date (≈ 14.4 14 C cal. kyr) has also been reported [START_REF] Pardo | 70 fechas para una cueva. revisión crítica de 70 dataciones c14 del Pleistoceno superior y Holoceno de la cueva de Nerja (Málaga, Andalucía, España). Espacio, Tiempo[END_REF].

Malladetes (or Mallaetes) is located in the central Mediterranean region of Spain, in

Valencia. The levels sampled for this study include levels 7 (entrada D), 8 and 9 (entrada D) which should be associated with the Gravettian (de la Peña Alonso, 2013). An AMS 14 C cal. age of 30-26 kyr was obtained from modern human remains associated with unspecified Gravettian levels in sector E, likely corresponding to levels VII to XI [START_REF] Fortea Perez | La Cueva de Les Mallaetes y los Problemas del Paleolítico Superior del Mediterráneo Español[END_REF][START_REF] Arsuaga | The Gravettian occipital bone from the site of Malladetes (Barx, Valencia, Spain)[END_REF], indicate that our samples probably predate the LGM and level 8 at Nerja [START_REF] Fumanal | Schéma paléoclimatique et chrono-stratigraphique d'une séquence du Paléolithique supérieur de la région de Valence (Espagne)[END_REF][START_REF] Tortosa | Cueva de Nerja (Málaga, Espana) y los inicios del solutrense en Andalucia[END_REF].

Sample collection

The samples recovered from the seven archaeological sites described above (Table 1) include mammal remains assigned to Cervidae (Cervus elaphus, Rangifer tarandus), Equidae (Equus sp.) and Caprinae (Capra pyreneica, Rupicapra rupicapra, Saiga tatarica). The cervids (Cervus elaphus) and caprids (Capra pyreneica) were sampled as a priority, but wherever these taxa were poorly represented, we sampled alternate species including smaller cervids and small bovids. At Laugerie Haute Est and Hornos de la Peña we sampled horse remains because of their abundance at both sites.

Habitat and diet of Cervidae, Equidae and Caprinae

Cervus elaphus (red deer) is a large mammal observed in semi-open habitats, often found in forest margins and clearings, while it also occupied open environments such as steppes during the Last Glacial [START_REF] Straus | On the habitat and diet of Cervus elaphus[END_REF][START_REF] Bugalho | The composition of the diet of red deer (Cervus elaphus) in a Mediterranean environment: a case of summer nutritional constraint?[END_REF]. It tolerates a wide range of temperature, as long as the habitat is not too dry. Its diet is highly eclectic and seasonally variable, although it shows a 60-70% herbaceous source. Cervus elaphus mainly browses leaves, conifer needles, twigs, nuts, shoots and bark of trees in winter, and grazes more on grasses, herbs, ferns, mosses, mushrooms, sedges or even rushes, during the summer and intermediate seasons.

Rangifer tarandus (reindeer) lives in subarctic and arctic regions (taiga and tundra), however some sub-species such as R. tarandus tarandus (woodland caribou) occupy the boreal forests [START_REF] Ihl | Habitat and diet selection by muskoxen and reindeer in western Alaska[END_REF][START_REF] Thompson | Factors influencing the seasonal diet selection by woodland caribou (Rangifer tarandus tarandus) in boreal forests in Ontario[END_REF] In winter, reindeers mainly feed on lichens which they access by digging holes in the snow cover. In the mild seasons, their diet is mainly based on grasses, green leaves of shrubs, tree bark and sedges as well as tree shoots such as larch, willow and birch. Mushrooms are also consumed during the late summer and early autumn.

Equus caballus (horse) was widely distributed in Eurasia and the Americas during the late Pleistocene. Equids are adapted to open habitats and are found in steppe, tundra and taiga habitats, as well as in more temperate grasslands. They are non-selective grazers and tolerate a wide range of temperatures [START_REF] Berger | Wild horses of the Great Basin: Social Competition and Population Size[END_REF].

Capra pyreneica (ibex) is a montane species living exclusively in open environments with high relief (cliffs and stone ravines; pastures in spring) at altitudes ranging from about 500 m to about 3500 m as observed in the Alps [START_REF] Villaret | Sexual segregation of habitat by the alpine ibex in the French Alps[END_REF]. Its diet is essentially composed of herbaceous plants.

Rupicapra rupricapra (chamois) feeds on herbs and lives in mountains meadows and on the edge of forests. Chamois are well-adapted to rough terrain and tolerate a large range of temperatures as low as -20°C in winter and up to 25°C in summer [START_REF] Christie | The sensitivity of chamois and red deer to temperature fluctuations[END_REF][START_REF] Morgia | Feeding habits, forage selection, and diet overlap in Alpine chamois (Rupicapra rupicapra L.) and domestic sheep[END_REF] although they are actually distributed in a wide range of temperate conditions.

Saiga tatarica (Eurasian antelope) feeds on salt-coated herbs, some grasses (Poaceae), Euphorbiaceae and Artemisia. This antelope lives in hilly areas and is considered an indicator of a dry climate (steppe environment) with highly-seasonally contrasted temperatures [START_REF] Bekenov | The ecology and management of the saiga antelope in Kazakhstan[END_REF].

Among these large mammals, obligate to semi-obligate drinkers are horse (strict obligate drinker), red deer (drinks large amounts of water during the warm season), reindeer (free water in summer, snow during the cold season) whilst semi-obligate to non-obligate drinkers are ibex, chamois and antelope that drink free water occasionally [START_REF] Macdonald | The Encyclopedia of Mammals[END_REF].

Sample preparation and isotopic analysis

The tooth samples were prepared and analyzed at the Stable Isotopes Laboratory of the Geotop research center (Université du Québec à Montréal). Prior to any treatment, tooth crowns were mechanically cleaned to remove any trace of remaining glue or altered hard tissues. They were then immersed in an ultrasonic bath in order to eliminate potential remains of sedimentary particles. From 10 mg to 20 mg of enamel powder samples were collected using a micro-drill equipped with a diamond-studded drill bit. Grooves were drilled perpendicularly to the tooth growth axis from the cervix to the apex of the tooth crown.

Following the protocol used by [START_REF] Vigeant | Dietary habits in New France during the 17th and 18th centuries: An isotopic perspective[END_REF], organic matter was removed from tooth enamel apatite by soaking an aliquot of the untreated sample in a 2% sodium hypochlorite (NaOCl) solution for 24 h. Then, the samples were rinsed with large amounts of water, then dried before being treated with 1 M acetic acid-Ca acetate buffer (pH = 4.75) for 24h at 20°C to remove secondary carbonate phases [START_REF] Bocherens | Isotopic biogeochemistry ( 13 C, 18 O) of mammalian enamel from African Pleistocene hominid sites[END_REF], and finally washed with distilled water and dried at ambient temperature. For each sample, about 1 mg of enamel powder was analyzed for δ 13 C-values (‰ VPDB) and δ 18 O-values (‰VPDB) of the carbonate using a Micromass Isoprime™ isotopic ratio mass spectrometer operated in dual inlet mode and coupled to a Multicarb™ system. Overall analytical uncertainties (±1σ) for δ 13 C and δ 18 O were close to 0.05‰ and 0.1‰, respectively, based on statistics from replicate measurement of samples, compared to ±0.02 and ±0.03‰, respectively, on reference carbonate materials (e.g., NBS 19, IAEA-603 and the home reference carbonate), thus suggesting some isotopic heterogeneity in the ground carbonate samples.

Results

The dataset is characterized by δ 13 C and δ 18 O values ranging respectively from -12.89‰ to -8.93‰ (mean = -10.42±0.76‰) and -8.65 to -0.14‰ (mean = -5.19±1.99‰). An examination of data compiled in Table 1 reveals that they form two distinct clusters corresponding to sites located on the Atlantic side and along the Mediterranean coast (Figure 1). As a first approach, pooling the taxa in each geographic cluster reveals that the δ 18 O values of mammal remains from Nerja and Malladetes have much higher δ 18 O values (-1.45±0.74) than those (-6.03±0.94) from the sites located on the Atlantic side (Table 1). The difference is not so striking for δ 13 C values which average -10.60±0.70‰ in the Atlantic sites, i.e., slightly but significantly lower than that the mean value of Mediterranean sites, which is -9.65±0.44‰ (Welch t-test: p = 4.4x10 -7 ; H0 = means are equal).

Discrimination between the datasets increases further when distinguishing Cervidae from

Caprinae in the cluster of Atlantic sites, which removes a sampling bias as only Caprinae were sampled in the Mediterranean sites. In the case of Cervidae, the distribution of data satisfies the statistical criteria required by a normal distribution when using Shapiro-Wilk test (p = 0.30 and W = 0.963 in the 95% critical range [0.937; 1.000] for δ 13 C values; p = 0.15 and W = 0.953 in the 95% critical range [0.937; 1.000] for δ 18 O values). Similar conclusions are obtained for the Caprinae when using Shapiro-Wilk test (p = 0.50 and W = 0.964 in the 95% critical range [0.920; 1.000]) for δ 13 C-values; p = 0.37 and W = 0.957 in the 95% critical range [0.920; 1.000] for δ 18 O-values).

Cervidae have δ 13 C values ranging from -12.89‰ to -9.47‰ (mean = -10.96±0.66‰) and δ 18 O values ranging from -8.65‰ to -3.67‰ (mean = -5.96±1.03‰) in the Atlantic sites.

These data statistically differ for δ 13 C values (Welch t-test: p = 1.03x10 -6 ; H0 = mean are equal), but not for δ 18 O values (Welch t-test: p = 0.71; H0 = mean are equal) of Caprinae coming from the same sites (mean δ 13 C = -10.12±0.49‰ and mean δ 18 O = -6.06±0.88‰).

Caprinae from the Atlantic sites have slightly lower δ 13 C-values (-10.12±0.49‰ against -9.65±0.44‰) and much lower δ 18 O-values (-6.06±0.88‰ against -1.45±0.74‰), than those originating from Mediterranean sites with δ 13 C-values of -9.65±0.44‰ and δ 18 O -values of -1.45±0.74‰ (Welch t-test for δ 13 C: p = 3.6x10 -3 ; H0 = mean are equal; Welch t-test for δ 18 O: p = 1.7x10 -16 ; H0 = mean are equal).

Few samples are available for Equidae: 2 samples from Hornos de la Peña and 2 from Laugerie-Haute, have δ 13 C and δ 18 O values in the range of values from other herbivorous mammals at sites located on the Atlantic side. It is worth pointing out that no temporal trend is observed for either δ 13 C-or δ 18 O-values of both Cervidae and Caprinae in the Atlantic sites (Table 1). Nevertheless, whatever their age, Caprinae from this study are characterized by lower δ 18 O values and higher δ 13 C values than Cervidae. In addition, statistics that show a normal distribution with relatively small standard deviations are consistent with inter-annual variations of considered climatic parameters.

Discussion

Preservation and meaning of the isotopic compositions

Assessing the potential blurring effect of diagenetic processes represents the greatest challenge for geochemists hoping to reconstruct environmental parameters based on isotopic proxies. Biogenic apatites have undergone diagenesis, at least during the first step of the mineral transformation commonly referred to as "early diagenesis", which follows the decay of the organic matter. Changes in the content of major and trace elements of biogenic apatites do not necessarily mean that their initial stable C-and O-isotope compositions have been fully reset or even partly modified [START_REF] Goedert | Euryhaline ecology of early tetrapods revealed by stable isotopes[END_REF]. The acquisition of secondary isotopic signatures mainly takes place during the first stages of burial, when processes of mineral dissolution and recrystallization under microbially-mediated conditions may alter the skeletal tissues [START_REF] Kolodny | Dinosaur bones: fossils or pseudomorphs? The pitfalls of physiology reconstruction from apatitic fossils[END_REF][START_REF] Trueman | Diagenetic effects on the oxygen isotope composition of bones of dinosaurs and other vertebrates recovered from terrestrial and marine sediments[END_REF][START_REF] Zazzo | Experimentally-controlled carbon and oxygen isotope exchange between bioapatites and water under inorganic and microbially-mediated conditions[END_REF]. Enamel is considered a more suitable material for paleoenvironmental or paleoclimatic reconstructions because its apatite carbonate fraction has been shown to be generally resistant to isotopic alteration during early diagenesis due to its highly dense and packed crystallites [START_REF] Zazzo | Experimentally-controlled carbon and oxygen isotope exchange between bioapatites and water under inorganic and microbially-mediated conditions[END_REF].

Nonetheless, the preservation of pristine isotopic records needs to be assessed before attempting paleoenvironmental interpretations of the oxygen and carbon isotope compositions of fossil apatites. Whereas there is no method for determining whether the original isotopic compositions are preserved, several tests may help assess the preservation state of the initial isotopic signatures. In the present study, three complementary approaches are used: i) Although δ 13 C-and δ 18 O-values may covariate in the case of specific climate variations, usually detected through geographic or temporal trends (e.g., a combination of higher precipitation and lower temperatures), a linear correlation generally indicates a mixing relationship between two endmembers: the original isotopic signature and a contaminant, such as soil carbonate. Our data do not reveal any such mixing (Atlantic sites {2 outliers}: n = 61; R 2 = 10 -3 ; p = 0.76 for H0: no correlation; Mediterranean sites: n = 13; R 2 = 3x10 -2 ; p = 0.57 for H0: no correlation).

ii) Loose clusters and scrambling of isotopic compositions should affect all tooth material independently of the taxa. Here, the systematic differences in δ 13 C-and δ 18 O-values of Cervidae vs Caprinae supports the hypothesis of a preserved isotopic imprint of their respective diets and habitats.

Reconstitutions of MAP

Method: principles and limits

The δ 13 C of consumed food by vertebrates mainly controls the δ 13 C of their skeleton [START_REF] Ambrose | Experimental evidence for the relationship of the carbon isotope ratios of whole diet and dietary protein to those of bone collagen and carbonate[END_REF]Ballasse et al., 1999;[START_REF] Szostek | Chemical signals and reconstruction of life strategies from ancient human bones and teeth-problems and perspectives[END_REF]. Carbon isotope ratios of apatite carbonate may be considered as a "dirty window" that opens on a combination of factors: the diet, the degree of openness of the environment (exposure to light) and the amount of precipitation. The reconstruction of MAP is thus done in two steps.

The first step consists in converting the δ 13 C of carbonate (δ 13 Ccarb) into the δ 13 C of diet (δ 13 Cdiet), and the second one in converting the δ 13 Cdiet into MAP values (mm.yr -1 ). As a MAP proxy, δ 13 Ccarb is basically tarnished with several important sources of uncertainty. Therefore, we consider that only broad estimates may be obtained. The commonly used isotopic δ 13 Ccarb -δ 13 Cdiet offset value (∆ 13 Cpl) of 14.1‰ [START_REF] Cerling | Carbon isotope fractionation between diet and bioapatite in ungulate mammals and implications for ecological and paleoecological studies[END_REF] cannot be used in this study as it varies from one species to another. More recently, Tejada-Lara et al. ( 2018) proposed a series of equations that consider the specificities in the digestive system as well as the body mass of different herbivorous mammals. Calculations of δ 13 Cdiet of the study taxa in Table 1 follow this model. It must be kept in mind that animal diets do not necessarily reflect the composition of the vegetation cover; the stricter an animal's food preferences, the larger the bias between the isotopic composition of the diet and that of the vegetation cover.

The second step consists in estimating MAP from the δ 13 C of vegetation, here assumed to be close to δ 13 Cdiet. The empirical relationship proposed by [START_REF] Rey | Late Miocene climatic and environmental variations in northern Greece inferred from stable isotope compositions (δ 18 O, δ 13 C) of equid teeth apatite[END_REF], which we will use below, is valid if the vegetation cover was composed of C3 plants. This is still an open debate, however, and conflicting reconstructions of the LGM paleovegetation of western Europe from models vs pollen data exist in the literature (e.g., [START_REF] Elenga | Pollen-based biome reconstruction for southern Europe and Africa 18,000 yr BP[END_REF][START_REF] Leng | Isotopes in Palaeoenvironmental Research[END_REF][START_REF] Woillez | Impact of CO2 and climate on the Last Glacial Maximum vegetation: results from the ORCHIDEE/IPSL models[END_REF][START_REF] Hatté | Excursions to C-4 vegetation recorded in the Upper Pleistocene loess of Surduk (Northern Serbia): an organic isotope geochemistry study[END_REF][START_REF] Binney | Vegetation of Eurasia from the last glacial maximum to present: Key biogeographic patterns[END_REF]; see a discussion and an explanatory hypothesis in Kaplan et al., 2016). Thus, our estimates of MAP from δ 13 C-data must be seen as a first approximation that needs to be confronted with estimates from other approaches. The method established by [START_REF] Rey | Late Miocene climatic and environmental variations in northern Greece inferred from stable isotope compositions (δ 18 O, δ 13 C) of equid teeth apatite[END_REF], using a least square regression from the dataset of [START_REF] Kohn | Carbon isotope compositions of terrestrial C3 plants as indicators of (paleo) ecology and (paleo) climate[END_REF], is based on relationship (1) below, established between plant leaf discrimination (∆ ) as defined above and mean annual precipitations (MAP):

+ 300 = 0.092 ±0.004 ∆ + 1.148 ±0.074 (1) with ∆ = !" $ %& ' ( ) %& ***
and MAP in mm yr -1 + ,-. values of atmospheric CO2 have been measured in air bubbles trapped within the Vostok ice cores between 20 and 43 ka [START_REF] Leuenberger | Carbon isotope composition of atmospheric CO2 during the last ice age from an Antarctic ice core[END_REF][START_REF] Schmitt | Carbon isotope constraints on the deglacial CO2 rise from ice cores[END_REF]. For the LGM, this value has been estimated to be close to -7‰ (VPDB). Calculations of MAP are given in Table 1. However, this equation does not consider the pCO2 effect on δ 13 C values of the vegetation [START_REF] Farquhar | On the relationship between carbon isotope discrimination and the intercellular carbon dioxide concentration in leaves[END_REF][START_REF] Schubert | The effect of atmospheric CO2 concentration on carbon isotope fractionation in C3 land plants[END_REF]. [START_REF] Hare | Atmospheric CO2 effect on stable carbon isotope composition of terrestrial fossil archives[END_REF] related the δ 13 C increase by ≈ 1‰ of plants and faunal collagen to the low pCO2 values (≈ 180 ppm) of the LGM relatively to those (≈ 300 ppm) of the post deglacial rise that took place about 15 kyr ago [START_REF] Peacock | A possible sequence of events for the generalized glacial-interglacial cycle[END_REF]. Consequently, it is most likely that our calculated MAP values according to equation ( 1) are underestimated by about 150 mm y -1 . Final estimates from the model will thus be corrected using this offset.

The validity and robustness of this empirical approach were tested successfully in a study [START_REF] Prud'homme | δ 13 C signal of earthworm calcite granules: A new proxy for palaeoprecipitation reconstructions during the Last Glacial in western Europe[END_REF]

Interpretation

The reconstructed δ 13 C-values of animal diet (δ 13 C diet ) suggest that they mainly fed on C3 plants as revealed by a mean δ 13 C diet of -23.5±1.0‰ for the whole data base (Table 1), thus supporting the scenario of a C3-plant dominated vegetational cover (compare with [START_REF] Hare | Atmospheric CO2 effect on stable carbon isotope composition of terrestrial fossil archives[END_REF]. The C3 Calvin-Benson pathway (~75 wt% of the present-day vegetal biomass including trees, temperate grasses, some forbs, shrubs) is the most discriminating regarding the heavy carbon isotopes, thus leading to low δ 13 C values ranging from -34‰ to -22‰ and averaging -26‰ [START_REF] Park | Carbon isotope fractionation during photosynthesis[END_REF][START_REF] Park | Metabolic fractionation of C 13 & C 12 in plants[END_REF][START_REF] Smith | Two categories of 13 C/ 12 C ratios for higher plants[END_REF][START_REF] O'leary | Carbon isotope fractionation in plants[END_REF].

Variations within this range of δ 13 C values result from external factors such as air temperature, water stress, exposure to light, δ 13 C of atmospheric CO 2 , pCO 2 and rate of CO 2 recycling (O' Leary, 1981). Worth mentioning is the fact that a δ 13 C diet of -23.5±1.0‰ could as well be compatible with a diet of up to ~ 20% of C4 plants, as illustrated for grazers from modern steppe environments [START_REF] Auerswald | Large regional-scale variation in C3/C4 distribution pattern of Inner Mongolia steppe is revealed by grazer wool carbon isotope composition[END_REF]. Interestingly, this study also demonstrates that despite a large overall variability of the relative abundance of C4-vs C3plants, no preferential grazing is observed.

Nevertheless, on the Atlantic side of the Iberian Peninsula, Caprinae have higher δ 13 C diet (mean = -22.8±0.5‰) than Cervidae and Equidae (mean for both = -24.3±0.6‰), all populations being characterized by normal distributions (Caprinae: p = 0.11 and W = 0.953 in the 95% critical range [0.944; 1.000]; Cervidae and Equidae: p = 0.19 and W = 0.960 in the 95% critical range [0.942; 1.000]). Such distinct carbon isotope compositions reflect diet differences related to distinct habitats or grazing habits (Figure 2). This observation is in agreement with [START_REF] Jones | Adaptability, resilience and environmental buffering in European Refugia during the Late Pleistocene: Insights from La Riera Cave (Asturias, Cantabria, Spain)[END_REF] who concluded that ibex adapted their niche during the

LGM.

Among Cervidae, C. elaphus have lower δ 13 C diet values (mean = -24.4±0.6‰) than T.

rangifer (mean = -23.2±0.36‰), even though only three samples were available of C.

elaphus. A higher δ 13 C diet for reindeers, relative to other animals, has been explained by their specific diet mainly based on lichens, especially consumed during the cold season [START_REF] Park | Carbon isotope fractionation during photosynthesis[END_REF][START_REF] Máguas | Spatial variation in carbon isotope discrimination across the thalli of several lichen species[END_REF][START_REF] Ben-David | Utility of stable isotope analysis in studying foraging ecology of herbivores: examples from moose and caribou[END_REF] Based on the estimated δ 13 C diet values, MAP values were calculated in the range ~ 100 to 400 mm.yr -1 for Atlantic sites (100 to 200 mm.yr -1 recorded by Caprinae and 200 to 400 mm.yr -1 recorded by Cervidae and Equidae), while they are close to 100 mm.yr -1 for the two Mediterranean sites (Table 2). The present-day reference values (Table 2) extracted from meteorological records (IAEA/WMO database) correspond to cities in the vicinity of the studied archaeological sites, specifically Almeria (202±66 mm.yr -1 ), for the Nerja site, Valencia (532±173 mm.yr -1 ), for the Malladetes site , Santader (1050±140 mm.yr -1 ), for the sites located in Asturias and Cantabria (Spain), and Cestas-Pierroton ( 937±192 mm.yr -1 ), for the Laugerie-Haute site. The Caprinae record of drier environmental conditions than those inferred from Cervidae and Equidae, suggests differences in diets and/or seasonal/altitudinal grazing habitats between some of the representative herbivorous species sampled.

Diet preferences between C3 and C4 may cause δ 13 C differences by several ‰ between sheep and goat, for example, as shown by [START_REF] Balasse | Distinguishing sheep and goats using dental morphology and stable carbon isotopes in C4 grassland environments[END_REF] for C4 grassland environments of the Central Rift Valley of Kenya. However, C4 was most likely absent or very scarce in the studied LGM-environments (Collins and Jones, 1985). Nevertheless, during the LGM, precipitation followed a pattern similar to the modern one (Figure 3), with sites under Atlantic vs Mediterranean influence defining two distinct clusters.

The Mediterranean cluster is significantly drier than the Atlantic one, but both depict much drier conditions than at present. LGM precipitation estimates obtained here broadly agree with pollen-based reconstructions (e.g., [START_REF] Peyron | Climatic reconstruction in Europe for 18,000 yr BP from pollen data[END_REF], but disagree with the wetter conditions reconstructed from climate model experiments (e.g., [START_REF] Ludwig | Regional atmospheric circulation over Europe during the Last Glacial Maximum and its links to precipitation[END_REF][START_REF] Arpe | A comparison of climate simulations for the last glacial maximum with three different versions of the ECHAM model and implications for summer-green tree refugia[END_REF].

Reconstitutions of MAAT

Method: principles and limits

Oxygen isotope ratios of apatite carbonate may be seen as a "frosted glass window" open on the composition of dietary sources of water and local meteoric water, ultimately used as a proxy for air temperatures (e.g. [START_REF] Fricke | Intra-tooth variations in δ 18O (PO4) of mammalian tooth enamel as a record of seasonal variations in continental climate variables[END_REF][START_REF] Daux | Oxygen isotope composition of human teeth and the record of climate changes in France (Lorraine) during the last 1700 years[END_REF]. The δ 18 O of consumed water by vertebrates mainly controls the δ 18 Oc of their skeleton [START_REF] Longinelli | Oxygen isotopes in mammal bone phosphate: a new tool for paleohydrological and paleoclimatological research[END_REF][START_REF] Kohn | Predicting animal δ 18 O: accounting for diet and physiological adaptation[END_REF]. The oxygen isotope composition of leaf water can somewhat blur the relationship between the δ 18 Oc of apatite carbonate and the δ 18 O of drinking water [START_REF] Ayliffe | Oxygen isotope composition of the bone phosphate of Australian kangaroos: potential as a palaeoenvironmental recorder[END_REF][START_REF] Cormie | Relationship between the hydrogen and oxygen isotopes of deer bone and their use in the estimation of relative humidity[END_REF][START_REF] Luz | Oxygen isotope variations in phosphate of deer bones[END_REF][START_REF] Pederzani | Oxygen isotopes in bioarchaeology: Principles and applications, challenges and opportunities[END_REF]. It constitutes, however, a good indicator of aridity (Levin et al., 2016;[START_REF] Blumenthal | Aridity and hominin environments[END_REF].

In this approach, the reconstitution of MAAT is done in four steps.

First, oxygen isotope compositions of apatite carbonate (δ 18 Oc) measured on the VPDB scale [START_REF] Coplen | Comparison of stable isotope reference samples[END_REF]) must be converted into those of their phosphate counterparts (δ 18 Op) before being used to estimate the composition of drinking water on the basis of the equations available in the literature. This step is achieved by using equation ( 2) that was determined by [START_REF] Iacumin | Oxygen isotope analyses of coexisting carbonate and phosphate in biogenic apatite: a way to monitor diagenetic alteration of bone phosphate?[END_REF] for modern herbivorous mammals. Unfortunately, this step is a source of large uncertainties (δ 18 Op: standard error of estimate = 0.73). We computed Equation (2) using the original data set published in [START_REF] Iacumin | Oxygen isotope analyses of coexisting carbonate and phosphate in biogenic apatite: a way to monitor diagenetic alteration of bone phosphate?[END_REF] in order to tighten uncertainties as follows:

δ 18 Op = 0.975(±0.031) δ 18 Oc -8.43(±0.81) (n = 17; R 2 = 0.985)

The second step consists in calculating the oxygen isotope composition of water consumed (free water and water derived from the food such as leaves and herbs) by the animals. Several oxygen isotope fractionation equations between apatite phosphate (δ 18 Op) and water (δ 18 Ow) have been published for large herbivorous mammals, covering all taxa investigated in the present study. It is worthy to note that animal physiology is integrated in the oxygen isotope fractionation equations [START_REF] Longinelli | Oxygen isotopes in mammal bone phosphate: a new tool for paleohydrological and paleoclimatological research[END_REF][START_REF] Luz | Oxygen isotope variations in phosphate of biogenic apatites, IV. Mammal teeth and bones[END_REF][START_REF] Kohn | Predicting animal δ 18 O: accounting for diet and physiological adaptation[END_REF][START_REF] Langlois | Box-modeling of bone and tooth phosphate oxygen isotope compositions as a function of environmental and physiological parameters[END_REF] that have been established on present-day fauna, either the same species or close relatives. The relevant equations have also been recalculated to set the δ 18 Ow variable as the estimator (Y-axis) in the regression line along with the associated statistics, as follows:

-Deer (Cervus elaphus) (D'Angella and [START_REF] D'angela | Oxygen isotopes in living mammal's bone phosphate: further results[END_REF]):

δ 18 Ow = 0.88(±0.03)δ 18 Op -22.50(±0.52) (n = 7; R 2 = 0.99)

-Reindeer (Rangifer tarandus) [START_REF] Iacumin | Relationship between δ 18 O values for skeletal apatite from reindeer and foxes and yearly mean δ 18 O values of environmental water[END_REF]:

δ 18 Ow = 1.21(±0.01)δ 18 Op -26.41(±1.31) (n = 34; R 2 = 0.74) (4)
Note that 2 outliers were removed from the original database for the setting of equation 4, as recommended by [START_REF] Iacumin | Relationship between δ 18 O values for skeletal apatite from reindeer and foxes and yearly mean δ 18 O values of environmental water[END_REF].

-Horse (Equus burchelli, E. zebra, E. caballus) [START_REF] Bryant | Oxygen isotopic composition of fossil horse tooth phosphate as a record of continental paleoclimate[END_REF]Sanchez Chillon et al., 1994;[START_REF] Delgado-Huertas | Oxygen isotope variations of phosphate in mammalian bone and tooth enamel[END_REF]:

δ
A major caveat here is the possibility of distinct δ 18 Ow vs δ 18 Op relationships during the cold LGM, due to reduced evapotranspiration (e.g., [START_REF] Kohn | Predicting animal δ 18 O: accounting for diet and physiological adaptation[END_REF]. Isotope enrichment from precipitation to plant-water and finally to the animal body-water could have prevailed. Thus, under LGM conditions, the δ 18 Ow-value calculated with the above equations must be seen as a maximum value. (1998), MAAT depicted large amplitude oscillations between <5°C and up to 15°C, with a mean value of about 7 to 9°C (e.g., [START_REF] Kuhlemann | Regional synthesis of Mediterranean atmospheric circulation during the Last Glacial Maximum[END_REF]. Comparatively, surface waters from the western Mediterranean were not only warmer but enriched in 18 O. Planktic foraminifera [START_REF] Hayes | Glacial Mediterranean sea surface temperatures based on planktonic foraminiferal assemblages[END_REF] and biomarkers [START_REF] Rodrigo-Gámiz | Sea surface temperature variations in the western Mediterranean Sea over the last 20 kar: A dual-organic proxy (UK′37and LDI) approach[END_REF], both suggesting a mean LGM sea-surface temperature (SST) of ~ 13.5°C, whereas δ 18 Ocarb values reported from this sector of the Mediterranean during the same interval range from +3.5±0.5‰ vs VPDB in planktic foraminifera [START_REF] Thunell | Glacial-Holocene salinity changes in the Mediterranean Sea: hydrographic and depositional effects[END_REF] to +3.9±0.2 ‰ in shallow mollusk shells (Fergusson et al., 2011). Combining these values, one may set the

LGM δ 18 Osw-value at about +3‰ vs VSMOW in the western Mediterranean, to be compared with a modern value of ~ +1.3‰ (Goddard Institute data base; Schmidt et al., 1999; see also [START_REF] Benetti | Stable isotopes in surface waters of the Atlantic Ocean: Indicators of ocean-atmosphere water fluxes and oceanic mixing processes[END_REF]. Potential moisture generated during the evaporation of the western Mediterranean Sea, which had a δ 18 O higher than the Atlantic Ocean by ~ +2‰, may thus have played some role during the LGM. However, climate models suggest a southward shift of the North Atlantic jet stream during the LGM, governing the relatively abundant winter precipitation in comparison with summer rains which were somewhat controlled by local processes (e.g., [START_REF] Beghin | What drives LGM precipitation over the western Mediterranean? A study focused on the Iberian Peninsula and northern Morocco[END_REF]. Thus, the most likely LGM moisture source, for both the We used the +0.8‰ offset between modern vs LGM moisture sources as a base line for both regions and subtracted it from the δ 18 Ow-values derived from equations (3) to (6) (see Table 1). Nevertheless, we also increased this offset to +1.2‰ for the Mediterranean sites, a value that would account for about 20% of Mediterranean moisture source supplies on a yearly basis, based on the ~ 2‰ offset in δ 18 Osw-values between the Bay of Biscay and the Western Mediterranean. Again, this seems a reasonable estimate considering the regional atmospheric circulation in SW Europe during the LGM (cf. [START_REF] Ludwig | Regional atmospheric circulation over Europe during the Last Glacial Maximum and its links to precipitation[END_REF]; see also [START_REF] Kuhlemann | Regional synthesis of Mediterranean atmospheric circulation during the Last Glacial Maximum[END_REF][START_REF] Beghin | What drives LGM precipitation over the western Mediterranean? A study focused on the Iberian Peninsula and northern Morocco[END_REF].

The fourth step aims at converting δ 18 Ow to MAAT. Here again, one may question whether or not present-day relationships fully apply to the LGM. Firstly, all of the sites of the present study are close to coastlines, thus, they are not affected by continental Rayleighevolution of humid air-masses or by major altitudinal effects (e.g., [START_REF] Gat | Observed isotope effects in precipitations[END_REF], a process well illustrated by [START_REF] Rozanski | Factors controlling stable isotope composition of European precipitation[END_REF] for inland precipitation patterns in Europe. Isotopic compositions of precipitation at the study sites should then follow a relatively linear relationship with mean annual temperatures, as observed at global scale, despite some latitudinal/regional variations in the slope and intercept of the local regression lines (e.g., [START_REF] Lee | Water isotopes during the Last Glacial Maximum: New general circulation model calculations[END_REF]. The slope of a MAT-δ 18 Omw linear relationship is mainly temperaturedependent at mid-and high-latitudes. This point is firstly well illustrated by the existence of a strongly positive linear correlation (R 2 = 0.83) between the dew point (Td) and the oxygen isotope composition of meteoric water (Figure 4). Td, expressed in °C, represents the temperature at which the actual water vapour pressure equals the saturated water vapour pressure (RH = 100%). This climatic parameter plays a key role in the distribution of δ 18 Omw over the continental masses located in the mid-latitudes [START_REF] Lécuyer | Geographic variations in the slope of the δ 2 H-δ 18 O meteoric water line over Europe: A record of increasing continentality[END_REF]. Secondly, there is a unique MAT-δ 18 Omw linear relationship for present-day (Figure 5A) for latitudes > 32°, even if it samples various continental areas from both hemispheres, which are submitted to various moisture-bearing prevailing winds.

The MAAT (in °C) vs δ 18 Ow linear relationship from equation ( 7), established for Europe [START_REF] Skrzypek | How cold was it for Neanderthals moving to Central Europe during warm phases of the last glaciation?[END_REF], is indistinguishable in terms of slope and intercept from equation ( 8) obtained at a global scale (Figure 5A). For a reliable comparison, equation ( 8) was computed by excluding data from the tropics as well as those corresponding to air temperatures lower than 0°C (Figure 5B).

MAAT (°C) = 1.41 δ 18 Ow-MAP + 21.63 (7) where MAT represents the mean annual temperature and δ 18 Ow-MAP, the weighted mean isotopic composition of annual precipitation.

MAAT (°C) = 1.40(±0.08) δ 18 Ow-MAP+ 21.86(±0.64) (8)

(n = 90; R 2 = 0.81) Equation ( 7) was thus used to calculate MAAT (reported in Table 1) during the LGM in southwestern Europe. There are three caveats when using this method:

1) Altitude effect: in Spain, the rate of air-temperature decrease with altitude is close to 5°C per 1000 m in mountain areas [START_REF] Navarro-Serrano | Estimation of near-surface air temperature lapse rates over continental Spain and its mountain areas[END_REF]. In the Pyrenees, altitudinal oxygen isotope gradients average -3‰ per 1000 m [START_REF] Hughes | Little Ice Age glaciers and climate in the Mediterranean mountains: a new analysis[END_REF], a value within the range observed in southern Europe (e.g., [START_REF] Tazioli | Estimating the Isotopic Altitude Gradient for Hydrogeological Studies in Mountainous Areas: Are the Low-Yield Springs Suitable? Insights from the Northern Apennines of Italy[END_REF]. Consequently, differences between present-day reference stations (cities are close to sea-level) and mammal ranges (between 0 and 1000 m) during the LGM may lead to overestimates of temperature differences by about 2°C in northern Spain and up to 4°C in south and southeast Spain.

2) Source effect: the source of drinking water may differ isotopically from local meteoric waters, for example in rivers since their sources are located at altitudes. A comparison between present-day river and precipitation δ 18 O maps [START_REF] Nan | Stable isotope composition of river waters across the world[END_REF]) reveals that differences between those two parameters are close to 2‰ over Europe on a yearly average. Once again, true air temperatures may be higher than those calculated following the above method. Similarly, animals that seasonally migrate may have a δ 18 Ow record that is biased relative to their place of death.

3) Reservoir effect: water stored in ponds, lakes or in snow cover may suffer evaporation or sublimation that tends to enrich the residual liquid-solid water in heavy isotopes. Such processes may induce oxygen isotope biases in the order of a few per mil.

However, a sporadic influence of the above boundary conditions on the isotopic approach should be detectable by altering the normal distribution of the studied population; which does not seem to be the case here.

Interpretation

δ 18 Ow-MAP have been calculated according either to ocean surface water having a δ 18 O of +1.2 or +0.8‰. Inferred temperatures using equation ( 7) are shifted by ≈ +0.6°C. In the case of a prevailing oceanic surface source of moisture having a δ 18 O of +1.2‰, calculated mean δ 18 Ow-MAP are -10.1±0.9‰ (n = 46) and -10.3±1.0‰ (n = 17) for northern Spain and Laugerie-Haute, respectively (Table 1). Pooled data follow a normal distribution (Figure 6) 2). Here also, this value is practically identical to LGM-SST values in the Western Mediterranean estimated by [START_REF] Hayes | Glacial Mediterranean sea surface temperatures based on planktonic foraminiferal assemblages[END_REF] or by [START_REF] Rodrigo-Gámiz | Sea surface temperature variations in the western Mediterranean Sea over the last 20 kar: A dual-organic proxy (UK′37and LDI) approach[END_REF].

according
The air temperatures reconstructed during the LGM are lower than the present-day ranges recorded in meteorological stations located either in the vicinity of the Atlantic sites (Santander: 14.8±0.4°C and Cestas-Pierroton: 13.1±0.6°C), or close to the Mediterranean sites (Almeria: 19.2±0.4°C; Valencia: 18.5±0.6°C). As shown in Figure 7, MAAT decreased during the LGM from South to North according to a trend roughly similar to the modern one.

As a whole, MAAT values in southern Spain were about 4-5°C below modern pre-industrial temperatures, and by 5-8°C in northern Spain and southwest France (Figure 7).

5.4.

LGM climate inferred from large mammal isotopes and biogeographic patterns.

Faunal associations described in the archaeological sites of southwestern Europe are generally compatible with cold and dry, open environments, such as tundra and dry steppe, but there are important regional differences. The sites sampled in this study are from three distinct biogeographic zones, following the classification suggested by [START_REF] Yravedra | Essai sur la biodiversité des associations de grands mammifères à la fin du Pléistocène dans le Sud-Ouest de l'Europe[END_REF] LGM, in zone IV (Cantabria/Asturias) diversity also decreases but reaches a minimum after the LGM, while carnivore diversity actually increases with time, whereas the opposite pattern occurs in zone IIIa (e.g., Valencia) where a gradual increase in herbivore diversity is observable over the course of the last Glacial [START_REF] Yravedra | Essai sur la biodiversité des associations de grands mammifères à la fin du Pléistocène dans le Sud-Ouest de l'Europe[END_REF].

In SW France (zone VI), horse dominates sites such as Laugerie Haute Est at the beginning of the LGM (e.g., in the Magdalenian 0 level; Table 1) but reindeer, a taxon associated with tundra environments, increases dramatically in importance throughout the sequence, eventually dominating during Magdalenian I and II (see [START_REF] Delpech | Biostratigraphie des niveaux solutréens de Laugerie-Haute (Les Eyzies, Dordogne, France). Implications archéologiques[END_REF]. In northwest Spain (zone IV), on the other hand red deer, a taxon generally considered more indicative of temperate, wooded contexts often dominates late Glacial faunal assemblages, while the proximity of sites to the coastal ranges is signaled by the presence of ibex and chamois [START_REF] Straus | The Pleistocene-Holocene transition in Cantabrian Spain: current reflections on culture change[END_REF]. Our results are consistent with the biogeographic patterns described above, indicating strong contrasts between Mediterranean sites (in Zone IIIa) and sites in zones IV and VI in terms of MAAT and MAP.

During the LGM, tundra extended across northern Europe while steppes dominated the Mediterranean areas, resulting in a geographic reduction and fragmentation of boreal evergreen forests (taiga) and temperate deciduous forests [START_REF] Tarasov | Last glacial maximum biomes reconstructed from pollen and plant macrofossil data from northern Eurasia[END_REF]. The mean value of 350±100 mm.yr -1 reconstructed for MAP from the δ 13 Ccarb of the dental apatite of herbivorous mammals (above) is consistent with this observation. However, a few samples of Cervus elaphus from zone IV, El Cierro and El Ruso, the closest to the Atlantic ocean depict a δ 13 Ccarb-value low enough to suggest that some years or periods in the year were characterized by MAP ≥ 500 mm.yr -1 . This could mean that the climate in northern Spain was rather unstable at inter-decadal or inter-centennial scales towards the end of the LGM, or simply that C. elaphus used a large seasonal altitudinal range, as observed in modern specimens in the Alps [START_REF] Luccarini | Red deer (Cervus elaphus) spatial use in the Italian Alps: home range patterns, seasonal migrations, and effects of snow and winter feeding[END_REF] or from the Rockies [START_REF] Morgantini | Nutritional Significance of Wapiti (Cervus Elaphus) Migrations to Alpine Ranges in Western Alberta, Canada[END_REF].

Geochemical proxies similar to those used in this study were applied to faunal remains (horse, deer and ibex) recovered from the archaeological site of Peña Capón [START_REF] Yravedra | Not so deserted… paleoecology and human subsistence in Central Iberia (Guadalajara, Spain) around the Last Glacial Maximum[END_REF], a rock shelter located in central Iberia which contains sedimentary deposits assigned to Marine Isotope Stage 2 (MIS 2). [START_REF] Yravedra | Not so deserted… paleoecology and human subsistence in Central Iberia (Guadalajara, Spain) around the Last Glacial Maximum[END_REF] estimated MAAT for archaeological levels with a 14 C cal. age of ~23.8 kyr, ranging from 9°C to 12°C, i.e. only slightly below the present-day MAAT of 12°C-13°C recorded in the locality of Muriel. These estimates for the central part of the Iberian Peninsula may reflect milder conditions prior to the onset of the

LGM. Nevertheless, it is interesting to note that they are bracketed by our results. It is also worthy to note that ibex recorded higher δ 13 C values and lower δ 18 O values than the other cohabitating mammals such as Cervidae and Equidae. These data point to a specific habitat for Caprinae as recently suggested by [START_REF] Jones | Adaptability, resilience and environmental buffering in European Refugia during the Late Pleistocene: Insights from La Riera Cave (Asturias, Cantabria, Spain)[END_REF] who analyzed both carbon and nitrogen isotope compositions of bone collagen. These authors concluded that the existence of refugia in south Europe allowed most herbivores to adapt despite harsh and fluctuating climatic conditions.

The LGM climate conditions have been reconstructed from palynological data (e.g., [START_REF] Peyron | Climatic reconstruction in Europe for 18,000 yr BP from pollen data[END_REF]. In the north of Spain (Lago Ajo) and in the south (Padul), MAAT values were estimated in the 10°C to 5°C range. Palynological estimates of MAP, range from 900-600 mm in the north of Spain, to close to 200 mm in the south. The temperatures from palynological records are congruent with the estimates of the present study. However, whereas paleoprecipitation estimates from palynology or from the present isotopic approach fit well in zone IIIa, on the Mediterranean side, the MAP estimates for Northern Spain from pollen proxies are somewhat higher than those based on δ 13 Ccarb values. This leads us to consider the possibility that isotopic reconstructions, such as used here, vs reconstructions from transfer function using palynological assemblages, do not necessarily encompass the same space-time dimensions. For example, differences in the seasonality of pollen production and patterns of altitudinal grazing may result in distinct signals.

As mentioned above, the study sites are close to either the Bay of Biscay (in zones IV SSTs on proximal inland MAAT. Moreover, these temperature estimates based on marine as well as on continental proxies, enforce the hypothesis of cold regional conditions as well as of a marked difference between the Atlantic and Mediterranean zones; the latter being less affected by global LGM cooling.

At odds with these results are the model outputs of paleo-glaciers [START_REF] Allen | Reconstructing glacier-based climates of LGM Europe and Russia-Part 2: A dataset of LGM precipitation/temperature relations derived from degree-day modelling of palaeo glaciers[END_REF], which provided both large amounts of precipitation (1000 to 1400 mm.yr -1 ) and temperature anomalies from -16°C to -11°C in northern and southern Spain, respectively. Such climatic patterns cannot be easily reconciled with the present knowledge of the LGM vegetation cover and faunal associations [START_REF] Sommer | Glacial refugia of mammals in Europe: evidence from fossil records[END_REF][START_REF] Benito Garzón | Predictive modelling of tree species distributions on the Iberian Peninsula during the Last Glacial Maximum and Mid-Holocene[END_REF] Morales-Molino and García-Antón, 2014; [START_REF] Binney | Vegetation of Eurasia from the last glacial maximum to present: Key biogeographic patterns[END_REF][START_REF] Quinzin | Glacial survival of trophically linked boreal species in northern Europe[END_REF].

Atmosphere-ocean global climate models (GCMs) computed global mean LGM-air temperatures 4°C to 7°C lower than Present, when adding feedback effects due to changes in the vegetation cover and the amount of atmospheric dust [START_REF] Crucifix | Impact of vegetation changes on the dynamics of the atmosphere at the Last Glacial Maximum[END_REF][START_REF] Jansen | Paleoclimate. Climate change 2007: the physical science basis; contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF][START_REF] Roche | Climate of the Last Glacial Maximum: sensitivity studies and model-data comparison with the LOVECLIM coupled model[END_REF][START_REF] Cao | Attribution of the Last Glacial Maximum climate formation[END_REF]. In southwest Europe, the air temperature anomaly was estimated between -8°C and -4°C based on a multi-model regression of proxy data and PMIP2 model output [START_REF] Annan | A new global reconstruction of temperature changes at the Last Glacial Maximum[END_REF]. Mean annual air temperatures anomalies were also estimated in the range -10°C to -5°C for ice-free areas of Europe, based on a regional climate model [START_REF] Strandberg | High-resolution regional simulation of last glacial maximum climate in Europe[END_REF]. [START_REF] Jost | High-resolution simulations of the last glacial maximum climate over Europe: a solution to discrepancies with continental palaeoclimatic reconstructions?[END_REF] emphasized the need of high-resolution spatial climate models to solve discrepancies with reconstructions based on proxies such as pollen assemblages. Models tend to underestimate the magnitude of cooling (especially in winter) with the consequence of overestimating the amount of precipitation.

When applied to Europe, most simulated climates, such as the those of the Paleoclimate Modelling Intercomparison Project (PMIP) or from refined General Climate Model (GCM), generate warmer and wetter LGM-conditions than those inferred from palynological or geochemical proxies [START_REF] Kageyama | The Last Glacial Maximum climate over Europe and western Siberia: a PMIP comparison between models and data[END_REF][START_REF] Latombe | Comparison of spatial downscaling methods of general circulation model results to study climate variability during the Last Glacial Maximum[END_REF]present study). [START_REF] Latombe | Comparison of spatial downscaling methods of general circulation model results to study climate variability during the Last Glacial Maximum[END_REF] study indicates temperatures and precipitation were close to 5-6°C and 200-250 mm y -1 in northern Spain, close to 7°C and 160 mm y -1 at Laugerie-Haute, and close to 11°C

and 225 mm y -1 at Nerja, in southern Spain. Within uncertainties, such air temperatures are coherent with those estimated from our isotopic proxy approach even though Latombe et al.

(2018) obtained slightly lower temperatures for Nerja. Our MAP estimates are generally slightly higher than those modeled by [START_REF] Latombe | Comparison of spatial downscaling methods of general circulation model results to study climate variability during the Last Glacial Maximum[END_REF], especially in the case of northern Spain.

Spatial resolution is critical in areas either close to the ocean (e.g., Atlantic Ocean) or along a mountain range (e.g., Pyrenees). [START_REF] Kageyama | Last Glacial Maximum temperatures over the North Atlantic, Europe and western Siberia: a comparison between PMIP models, MARGO sea-surface temperatures and pollen-based reconstructions[END_REF] performed a comparison of PMIP models, SSTs from MARGO (Multiproxy Approach for the Reconstruction of the Glacial Ocean surface) and pollen-based reconstructions over Europe. These authors highlighted the main data-model discrepancy observed for western Europe, with winter temperatures underestimated by all models that could be at least partly explained by a high interannual variability in the temperatures of the coldest months. However, more recent highresolution climatic simulations tend to reduce the discrepancy vs geochemical proxy-derived temperature estimates. This is for example the case of the central part of Iberia (~ 21 ka) with a computed MAAT close to 7±3°C [START_REF] Burke | Exploring the impact of climate variability during the Last Glacial Maximum on the pattern of human occupation of Iberia[END_REF] that are comparable or only slightly higher than those inferred from the oxygen isotope composition of large mammals [START_REF] Yravedra | Not so deserted… paleoecology and human subsistence in Central Iberia (Guadalajara, Spain) around the Last Glacial Maximum[END_REF].

The geochemical data presented in our study do not provide a temporal resolution at the seasonal scale, but large inter-annual variations in winter air temperatures should be imprinted on the half-year or full year time record provided by mammal tooth enamel. For the studied archaeological sites located on the Atlantic side, the range of δ 18 Ow-MAP is 4‰ (Figure 6), which corresponds to a temperature variation of ±5.5°C around the mean, according to equation ( 7). If interpreted as inter-decadal or inter-centennial variations in MAAT, such a magnitude is higher than that observed (2±1°C) during periods of relative climate steadiness (e.g., [START_REF] Corona | Millennium-long summer temperature variations in the European Alps as reconstructed from tree rings[END_REF][START_REF] Ljungqvist | Centennial-scale temperature change in last millennium simulations and proxy-based reconstructions[END_REF]. Whether this temperature variation of Climate conditions on the Iberian Peninsula are driven by large-scale processes that create strong regional contrasts in rainfall and temperature values. Our results demonstrate that these contrasts are reliably recorded in hard tissues and that isotope analyses of vertebrate remains can be used as proxies of mean annual temperature and precipitations on a regional scale. We confirm that, during the LGM, populations living in Franco-Cantabria and on the Mediterranean coast faced different environmental challenges. As the paleontological and archaeozoological records demonstrate it [START_REF] Yravedra | Essai sur la biodiversité des associations de grands mammifères à la fin du Pléistocène dans le Sud-Ouest de l'Europe[END_REF] this difference impacted the relative availability of different prey species, which would have required adjustments of human hunting strategies, potentially driving technological change. Indeed, the

LGM coincides with the appearance of the Solutrean technocomplex, which is geographically limited to the Iberian Peninsula and France. Three Solutrean "culture areas" (Andalucia, Portugal and Vasco-Cantabria) are described for the Iberian Peninsula [START_REF] Straus | A quarter-century of research on the Solutrean of Vasco-Cantabria, Iberia and beyond[END_REF] and these differences in the material culture used by human groups could reflect regional adaptations to contrasting environmental conditions such as those described here. These results are also consistent with archaeological research using an agent-based model that suggests that the population of Western Europe was spatially fragmented during the LGM as a result of climate-driven differences in habitat suitability (Wren and Burke, 2018). In this latter study, the SE region of the Iberian Peninsula was described as a "sink" during the LGM, with human population numbers declining periodically and being replenished from a larger "source" population, located in SE France and NE Spain.

Considering the uncertainties associated with our stable isotope proxy approach, including the altitudinal effect that is not easy to constrain, the present MAAT estimates are coherent with those inferred from other proxies (pollens, foraminifera) and some recent model outputs. In the case of MAP reconstructions, the isotopic approach suggests dry environments in agreement with our knowledge of LGM-flora and faunal associations. However, most modeling approaches still generate divergent patterns with some of them suggesting paleoenvironments as wet as those of today, or even more, along with similar MAAT.

Future research priorities may focus on the initial conditions required by the GCM models to reproduce the dry and cold conditions of the LGM in southwestern Europe. In addition, reconstructing the amplitude of precipitation and air temperature seasonality by high-resolution incremental sampling of tooth enamel may be of great interest. Indeed, mammals could have fed on resources (water, vegetation) mainly renewed during the winter and intermediate seasons, not in summer, hence biasing the isotopic record at least toward too negative mean air temperatures.

Conclusions

Our study documents the carbon and oxygen isotope composition of tooth enamel carbonate of Cervidae, Equidae and Caprinae remains dated to the LGM in southwest Europe.

The δ 13 C of tooth enamel carbonate was ultimately used to obtain semi-quantitative estimates of MAP. Precipitation in the range of ~ 250 to 550 mm.yr -1 were obtained for Atlantic sites and values close to 250 mm.yr -1 for the Mediterranean sites; both estimates are lower than local precipitation today (200-1000 mm.yr -1 ). Despite the large uncertainties associated with each step of the method, the reconstruction of a dry and cold LGM environment is consistent with the faunal and floral associations documented for this period. The limitation of the δ 13 C method for estimating paleo-precipitation is mainly related to the behavior of the studied vertebrate species. The observed δ 13 C differences between Caprinae and Cervidae or Equidae, for example, most likely reflects habitat use, with Caprinae occupying possibly drier habitats.

The δ 18 O of tooth enamel carbonate was used for estimating MAAT. In southern Spain, temperatures during the LGM were lower by about 4-5°C relative to modern times and lower by 5-8°C in northern Spain and southwest France. Several steps are necessary to calculate MAAT, and each one potentially propagates sizable uncertainties due to the combined use of empirical linear equations relating the isotopic compositions of the biomineral, the meteoric water and the air temperature. Moreover, we identified additional sources of uncertainties.

The first two are related to an "altitude effect" and a "source effect" that lead to an overestimate of the air temperature differences between the LGM and the present day on the order of a few °C. The third one, which we refer to as the "reservoir effect", relates to the 18 Oenrichment of meteoric waters stored in a natural reservoir (pond, lake, snow cover) open to the atmosphere and subsequently exposed to evaporation. In this case, the lower the residence time of water in the reservoir, the higher the increase in the δ 18 O of stored water, hence leading to a potential overestimation of air temperatures.

At first glance, our estimates of MAP and MAAT are compatible with the faunal associations described in the archaeological sites of southwestern Europe with generally cold

and dry, open environments such as tundra and dry steppe. Sea surface temperatures inferred from the oxygen isotope composition of planktonic foraminifera enforce the hypothesis of cold regional conditions as well as a marked difference between the Atlantic and Mediterranean sides; the latter being less impacted by the global cooling.

Close to montane areas, topography is highly variable with a juxtaposition of microclimates that makes a comparison of geochemical proxy results with model outputs risky and potentially fruitless, considering the unresolved differences in temporal and spatial resolution inherent in both scientific approaches. Nevertheless, there is room for improvement with, for example, the reconstitution of the seasonal amplitude of precipitation and air temperature performed by high-resolution incremental sampling of vertebrate tooth enamel. Precipitation and air temperature for the coldest and warmest months could be directly compared to model output with the aim to identify potential biases intrinsic to each method.

Québec Société et Culture, grant #SE-179537 (to AB). A NSERC-Canada Discovery Grant to chm and an infrastructure award from the Fonds de Recherche du Québec en Sciences Naturelles et Technologie to Geotop supported the analytical work. Three anonymous reviewers are thanked for their comments and suggestions that helped to improve our interpretations.

Table captions:

Table 1: Carbon and oxygen isotope compositions of tooth enamel from large herbivorous mammal remains of the LGM recovered from southwest France (Laugerie-Haute), north Spain (El Cierro, El Ruso, El Pendo, Hornos de la Peña), east Spain (Malladetes) and southeast Spain (Nerja). See the Figure 1 for the geographic locations of these archaeological sites. Studied taxa and 14 C ages are also provided along with the calculated oxygen isotope compositions of drinking and diet water (δ 18 Ow) and the carbon isotope composition of their diet (δ 13 Cdiet). See the text for more explanations.

Calendar ages were determined using the most recent calibration curve data for the Northern Hemisphere, IntCal20 [START_REF] Reimer | The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0-55 cal kBP)[END_REF]. Calibrated date ranges correspond to 92.9% probability (2σ).

2.

  Archaeological sites and sample collection 2.1. Description of the sampled archaeological sites Laugerie Haute Est (LHE) forms part of a much larger rock-shelter on the right bank of the Vézère river in the Dordogne region of France (Fig. 1). The shelter, which is about 30 m deep and 180 m long, is separated into two excavated sections (East and West) due to the presence of a house occupying the centre of the shelter. First discovered at the end of the 19th

  devoted to the reconstruction of MAP in Germany (Nussloch) during the Last Glacial, which was based on the δ 13 C of earthworm granule calcite crystals. In this study, the estimated MAP values range from ~ 270 to 460 mm with a mean value of ~ 370 mm for the Upper Pleniglacial (~LGM) and of ~ 300 mm, for the Middle Pleniglacial, covering a time interval between 17 ka and 55 ka. Uncertainties associated with the estimated MAP values are on the order of ± 100 to 200 mm. Based on δ 13 C values of organic matter preserved in the Nussloch loess deposits and using an inverse modelling approach (BIOME 4), Hatté and Guiot (2005) estimated MAP values ranging from 280 $/ ' 0 to 460 $ / '2 mm. Considering the large uncertainties associated with the various methods of MAP reconstructions, Prud'homme et al. (2018) concluded that the empirical method developed by Kohn (2010) may be used for MAP estimates in the absence of other proxies or modelling approaches.

  δ 13 C differences among C3-eating large mammals may result from differences in water stress and/or exposure to light (O'Leary, 1981). It means that δ 13 C values of Caprinae should be higher than those of Cervidae and Equidae if the Caprinae habitat was open (sparse vegetation cover in a rocky environment) and characterized by a vegetation submitted to water stress (e.g. dry sides of mountains).

  2 outliers were removed from the original database deduced from linear regression analysis -Wild Goat (Capra hircus and C. ibex) (Delgado-Huertas et al., 1995) δ 18 Ow = 1.09(±0.06)δ 18 Op -26.45(±0.98) (n = 21; R 2 = 0.95)

  The third step consists in correcting the δ 18 O of seawater from the volume of continental ice. Indeed, any change in the isotopic composition of oceanic sources will propagate throughout the entire surface water cycle. Knowing that the mean depth of the oceans is 3,700 m[START_REF] Charette | The volume of Earth's ocean[END_REF], a decrease in sea level by 120±10 m[START_REF] Lambeck | Sea level and global ice volumes from the Last Glacial Maximum to the Holocene[END_REF] results in a mean increase of the global ocean reservoir δ 18 O-value by +1.2±0.1‰ provided the mean δ 18 O of the LGM ice-sheets is close to -40‰[START_REF] Blunier | Timing of millennial-scale climate change in Antarctica and Greenland during the last glacial period[END_REF]. This isotopic offset is very close to that proposed by[START_REF] Schrag | The oxygen isotopic composition of seawater during the Last Glacial Maximum[END_REF] and[START_REF] Steele | Elements of physical oceanography: a derivative of the encyclopedia of ocean sciences[END_REF], but perhaps a little overestimated, based on isotopic composition for the Laurentide Ice Sheet proposed by[START_REF] Hillaire-Marcel | The late Pleistocene Laurentide glacier: ThU dating of its major fluctuations and δ 18 O range of the ice[END_REF] or[START_REF] Sima | Modeling the oxygen-isotopic composition of the North American Ice Sheet and its effect on the isotopic composition of the ocean during the last glacial cycle[END_REF].[START_REF] Duplessy | Constraints on the ocean oxygen isotopic enrichment between the Last Glacial Maximum and the Holocene: Paleoceanographic implications[END_REF] proposed an offset of 1.05±0.20‰ probably encompassing the effective overall change in the ocean isotopic composition. However, moisture sources relate to ocean surface water masses: their isotopic composition differs from that of the whole ocean, and may show large regional offsets, especially in the present case, between the Bay of Biscay and the western Mediterranean. The Bay of Biscay was likely the major moisture supplier for the Atlantic sites, whereas a mixture of Atlantic and western Mediterranean moisture supplies may be hypothesized for the Mediterranean sites. The two marine areas were of course characterized by distinct MAAT, mean salinity and isotopic composition.Based on[START_REF] Schrag | The oxygen isotopic composition of seawater during the Last Glacial Maximum[END_REF], one may assume an isotopic composition of surface waters from the northeast North Atlantic in the +0.7±0.1 to +0.8±0.1‰ range vs VSMOW, a range compatible with model estimates from[START_REF] Paul | Oxygen isotopic composition of the Mediterranean Sea since the Last Glacial Maximum: constraints from pore water analyses[END_REF]. According toLoncaric et al.

  Mediterranean and the Atlantic sites, would have been the northeastern North Atlantic Ocean, leading to adopt the δ 18 Osw ~ +0.8‰ value for the primary source water, and allowing for some18 O-enriched precipitation in summer, influenced by the high δ 18 Osw-values of the western Mediterranean waters (~ + 3‰). Perhaps worth of mention here is the fact that under the modern conditions, i.e., with a ∆δ 18 O offset of ~ +0.8‰ between surface waters of the western Mediterranean (δ 18 Osw ~ 1.3‰) and those of the northeastern North Atlantic (δ 18 Osw ~ 0.5‰;[START_REF] Benetti | Stable isotopes in surface waters of the Atlantic Ocean: Indicators of ocean-atmosphere water fluxes and oceanic mixing processes[END_REF], the long-term isotopic composition of precipitation in both areas is almost similar (~ -6‰;[START_REF] Araguas-Araguas | Isotope composition of precipitation and water vapour in the Iberian Peninsula: first results of the Spanish network of isotopes in precipitation[END_REF], thus pointing even today to a major contribution of Atlantic moisture sources to precipitation on the Iberian Peninsula.

  to the Shapiro-Wilk test (p = 0.26 and W = 0.976 in the 95% critical range [0.9620; 1.000]). The range of δ 18 Ow-MAP values comprised between ~ -12‰ and ~ -8‰ correspond to present-day compositions of precipitation in south Scandinavia and Central Europe (IAEA/WMO database and maps). Calculated mean δ 18 Ow-MAP are -5.4±0.8‰ (n = 8) and -5.8±0.8‰ (n = 6) for Malladetes and Nerja, respectively. Pooled data follow a normal distribution (Figure 6) according to the Shapiro-Wilk test (p = 0.93 and W = 0.975 in the 95% critical range [0.875; 1.000]). The range of δ 18 Omw values between ~ -7‰ and ~ -4‰ correspond to present-day compositions of precipitation along the Mediterranean rim (IAEA/WMO database and maps). The calculated δ 18 Ow-MAP values based on a δ 18 Osw of +0.8‰ range from ~ 6°C to ~ 10°C for Atlantic sites. This range may be explained by altitudinal gradients recorded by the faunal remains and/or by the large amplitude oscillations of LGM-SST and LGM-δ 18 Ocarb values in the Bay of Biscay reported by[START_REF] Loncaric | Late Quaternary sedimentation patterns on the Meriadzek Terrace, Bay of Biscay (ESSCAMP 02 core: 47°N 9°W)[END_REF]. It is interesting that the range of temperature inferred for the Atlantic sites in the present study is within the range proposed by these authors for surface waters of the Bay of Biscay during the LGM. On the Mediterranean side, the estimated LGM-MAAT values are close to 14-15°C for both sites (δ 18 Osw of +0.8‰) or close to 13.5-14.5°C (δ 18 Osw of +1.2‰) when accounting for some contribution of Mediterranean moisture to summer precipitations (Table

  and VI) or the western Mediterranean (zone IIIa).[START_REF] Hayes | Glacial Mediterranean sea surface temperatures based on planktonic foraminiferal assemblages[END_REF] and Rodrigo-Gámiz et al. (2014) estimated LGM-SSTs of ~ 13.5°C off the Mediterranean coastline of Spain, whereas Kucera et al. (2005) estimated LGM-SST of ~ 8°C for the Bay of Biscay on the basis of foraminifera δ 18 O values. Both values fit with the corresponding inland, but coastal temperatures reconstructed here, highlighting the strong influence of adjacent marine basin

  ±5.5°C indicates high-frequency climate instabilities in western Europe during the LGM, in agreement with Jones et al. (2020)'s or Loncaric et al. (1998), remains to be further investigated.

Figure 1 :

 1 Figure captions:

Figure 2 :.

 2 Figure 2: Frequency histogram of the calculated δ 13 C-value of mammal diets. Cervidae andEquidae are represented in blue; Caprinae, in red. The Gaussian curves were obtained by parametric estimation using PAST4.02 software[START_REF] Hammer | PAST: Paleontological Statistics Software Package for Education and Data Analysis[END_REF]. These curves may slightly differ from those obtained by least-squares approximations (performed using 10,000 iterations) according to the general equation of normal distribution as

Figure 3 :

 3 Figure 3: A) Modern and reconstructed MAP (mm yr -1 ) for the various French and Spanish archaeological sites. Modern cities (Cestas-Pierroton for Laugerie-Haute, Santander for northern Spain, Valencia for Malladates and Almeria for Nerja) in the vicinity of studied sites were ranked according to increasing MAP that roughly follow a S-N geographic trend. See Table2for the output of calculation. B) estimated MAP taking into account the lower atmospheric pCO2 during the LGM relatively to the Holocene (see section 5.2.1 for a detailed explanation).

Figure 4 :

 4 Figure 4: Present-day linear relationship between the dew point temperature in °C (Td) and the mean weighted δ 18 O-value (δ 18 Omw) of precipitations, extracted from the IAEA/WMO database. The database was filtered to only include meteorological sites above |32°| latitude. Td is the dew point (°C). The equation of the linear regression is Td = 1.30(±0.06)δ 18 Omw + 15.69(±0.67) with R 2 = 0.83.

Figure 5 :

 5 Figure 5: Present-day MAAT-δ 18 Omw linear relationships extracted from the IAEA/WMO database. A) Global data above 32°N with 15 outliers removed after linear regression analysis: MAAT = 1.42(±0.04)δ 18 Omw + 22.02(±0.44); B) Global data for sites above 32°N and with a MAAT > 0°C: MAAT = 1.40(±0.08)δ 18 Omw + 21.86(±0.64).

Figure 6 :

 6 Figure 6: Frequency histogram of δ 18 O of meteoric waters calculated from the δ 18 O of apatite carbonate from the sampled teeth of Cervidae, Equidae and Caprinae. Archaeological sites located from northern Spain and SW France, are represented in blue while the Mediterranean sites of Malladates and Nerja are represented in red. The Gaussian curves were obtained by parametric estimation using PAST4.02 software[START_REF] Hammer | PAST: Paleontological Statistics Software Package for Education and Data Analysis[END_REF]. These fitting curves may slightly differ from those obtained by least-squares approximations (performed by using 10,000 iterations) according to the general

Figure 7 :

 7 Figure 7: Modern and reconstructed MAAT (°C) of the sampled archaeological sites. Modern cities (Cestas-Pierroton for Laugerie-Haute, Santander for northern Spain, Valencia forMalladates and Almeria for Nerja) in the vicinity of studied sites were ranked according to a decreasing MAAT that roughly follow a S-N trend. See Table2for the output of calculation. Temperatures were calculated assuming an oxygen isotope composition of oceanic surface waters of +0.8‰ (VSMOW).

  for southwestern Europe. The zones sampled for this study are: southwest France (zone VI), Cantabria/Asturias (zone IV) and the Mediterranean coast (zone IIIa). A general reduction in faunal diversity is observable in Southwestern Europe during the LGM, signaling the severity of the climate downturn, but the relative diversity of herbivores and carnivores differs between zones. Diversity decreases in zone VI reaching minimum levels during the

Table 2 :

 2 Synthesis of calculated mean annual air temperatures (MAAT) and mean annual precipitation (MAP) during the LGM, estimated. from the δ 13 C-and δ 18 O-values of tooth enamel carbonate from large herbivorous mammal remains.

18 Oc δ 18 Op calc. d 18 Ow cor.ice (1.2) d 18 Ow cor.ice (0.8) T (d 18 Osw = 1.2) T (d 18 Osw = 0.8)

  d 13 Cc: carbon isotope composition of apatite carbonate d 13 Cdiet: calculated carbon isotope composition of the vertebrate diet ∆ 13 C leaf: d 13 Ccarb -d 13 Cdiet offset value P: calculated amount of precipitation (mm) d 18 Oc: oxygen isotope ratio of apatite carbonate d 18 Op calc.: calculated oxygen isotope composition of apatite phosphate d 18 Osw: oxygen isotope composition of Atlantic oceanic surface waters during the LGM d 18 Ow cor. Ice (+1.2): calculated oxygen isotope composition of vertebrate diet water corrected from contiental volume ice (corresponding to a d 18 Osw = +1.2‰) d 18 Ow cor. Ice (+0.8): calculated oxygen isotope composition of vertebrate diet water corrected from contiental volume ice (corresponding to a d 18 Osw = +0.8‰) T (d 18 Osw = +1.2): calculated mean air temperature assuming a d 18 O of +1.2‰ for the oceanic surface waters T (d 18 Osw = +0.8): calculated mean air temperature assuming a d 18 O of +0.8‰ for the oceanic surface waters

	Site	Lab#	Level	14 C dates BP	+/-	14 C cal BP δ IntCal20 (92.9%) Lithic scale Taxon Material δ 13 Cc δ 13 Cdiet ∆ 13 C leaf P δ 18 Oc ‰ VPDB ‰ VPDB (mm) ‰ VPDB ‰ VSMOW ‰ VSMOW	‰ VSMOW	‰ VSMOW	(°C)	(°C)
	El Cierro	EC-1B	G	15580	75	19011-18735	Lower Magdalenian	Cervus elaphus	M3 lower D	-11.21	-24.55	17.99	336	-6.6	24.1	15.1	-10.4	-10.0	6.9	7.5
	El Cierro	EC-2B	G	15580	75	19011-18735	Lower Magdalenian	Cervus elaphus	M3 lower D	-11.12	-24.46	17.89	323	-6.1	24.5	15.5	-10.0	-9.6	7.4	8.0
	El Cierro	EC-3B	G	15580	75	19011-18735	Lower Magdalenian	Cervus elaphus	M3 lower D	-10.94	-24.28	17.71	298	-6.6	24.0	15.0	-10.5	-10.1	6.8	7.4
	El Cierro	EC-4B	G	15580	75	19011-18735	Lower Magdalenian	Cervus elaphus	M3 lower D	-11.42	-24.76	18.22	367	-5.9	24.8	15.8	-9.8	-9.4	7.7	8.3
	El Cierro	EC-5B	G	15580	75	19011-18735	Lower Magdalenian	Cervus elaphus	M3 lower D	-10.92	-24.26	17.69	296	-6.7	24.0	15.0	-10.5	-10.1	6.8	7.4
	El Cierro	EC-17B	G	15580	75	19011-18735	Lower Magdalenian	Cervus elaphus	M2 lower D	-12.04	-25.38	18.86	463	-6.2	24.5	15.5	-10.1	-9.7	7.4	8.0
	El Cierro	EC-19B	G	15580	75	19011-18735	Lower Magdalenian	Cervus elaphus	M2 lower G	-12.28	-25.62	19.11	505	-5.5	25.1	16.1	-9.5	-9.1	8.2	8.8
	Malladetes	MA-2B	7 entry D	25120	240	> 25822-24366 < 29969-28917	Gravettian	Capra pyreneica	M3 sup G	-9.74	-22.58	15.94	111	-1.8	29.0	19.9	-6.0	-5.6	13.2	13.8
	Malladetes	MA-3B	7 E	25120	240	> 25822-24366 < 29969-28917	Gravettian	Capra pyreneica	M3 lower G	-10.84	-23.68	17.08	225	-2.2	28.6	19.5	-6.4	-6.0	12.5	13.1
	Malladetes	MA-4B	7 E	25120	240	> 25822-24366 < 29969-28917	Gravettian	Capra pyreneica	M3 sup D	-9.70	-22.54	15.90	108	-1.3	29.5	20.4	-5.4	-5.0	13.9	14.5
	Malladetes	MA-5B	8 E	25120	240	> 25822-24366 < 29969-28917	Gravettian	Capra pyreneica	M3 lower D	-9.93	-22.77	16.13	129	-0.7	30.1	21.0	-4.8	-4.4	14.9	15.5
	Malladetes	MA-6B	8 E	25120	240	> 25822-24366 < 29969-28917	Gravettian	Capra pyreneica	M3 lower G	-9.73	-22.57	15.93	111	-0.7	30.1	21.0	-4.8	-4.4	14.9	15.5
	Malladetes	MA-7B	9 entry D	25120	240	> 25822-24366 < 29969-28917	Gravettian	Capra pyreneica	M3 sup G	-9.84	-22.68	16.05	121	-2.0	28.8	19.7	-6.2	-5.8	12.9	13.5
	Malladetes	MA-8B	9 entry D	25120	240	> 25822-24366 < 29969-28917	Gravettian	Capra pyreneica	M3 lower G	-9.62	-22.46	15.81	101	-0.1	30.7	21.6	-4.1	-3.7	15.8	16.4
	Nerja	NE-1B	Vest VIIIc A-7	17940	200	> 22311-21181 < 23420-20987	Solutrean	Capra pyreneica	M3 lower G	-9.49	-22.33	15.68	90	-1.3	29.5	20.4	-5.5	-5.1	13.9	14.5
	Nerja	NE-2B	Vest VIIIc A-7	17940	N/A > 22311-21181 < 23420-20987	Solutrean	Capra pyreneica	M3 lower G	-9.25	-22.09	15.43	69	-0.7	30.2	21.1	-4.7	-4.3	14.9	15.5
	Nerja	NE-4B	Vest VIIIc A-7	17940	N/A > 22311-21181 < 23420-20987	Solutrean	Capra pyreneica	M3 lower D	-9.07	-21.91	15.25	55	-1.8	29.0	19.9	-6.0	-5.6	13.2	13.8
	Nerja	NE-6B	Vest VIIIc A-7	17940	N/A > 22311-21181 < 23420-20987	Solutrean	Capra pyreneica	M2 lower D	-9.21	-22.05	15.39	66	-1.6	29.2	20.1	-5.8	-5.4	13.4	14.0
	Nerja	NE-3B	N84 VIIIe A7-8 & B7-8	17940	N/A > 22311-21181 < 23420-20987	Solutrean	Capra pyreneica	M3 lower G	-9.52	-22.36	15.71	92	-2.8	28.0	18.9	-7.0	-6.6	11.7	12.3
	Nerja	NE-5B	N84 VIIIe A7-8 & B7-8	17940	N/A > 22311-21181 < 23420-20987	Solutrean	Capra pyreneica	M3 lower D	-9.52	-22.36	15.71	92	-1.7	29.1	20.0	-5.9	-5.5	13.3	13.9

Acknowledgements -We are grateful to A. Álvarez-Fernández and J. Bécares (Universidad de Salamanca), E. Aura Tortosa and E. Muñoz Fernández (Colectivo para la Ampliación de Estudios de Arqueología Prehistórica), C. Cacho (Museo Arqueológico Nacional), A.

Chauvin and R. Ontañón (Museo de Prehistoria y Arqueología de Cantabria), J-J. Cleyet-Merle and S. Madeleine (Musée national de Préhistoire), M. Pérez Ripoll (Universitat de València) and J. Yravedra (Complutense University) for their great and kind help during the first stages of this study, as well as for providing access to samples and allowing us to work in their different laboratories. We warmly thank J. Lepage for her help in the data collection, who was awarded an MEESR mobility grant. This study was possible thanks to an Early Career Research grant to MAJ. This research was supported by the Fonds de Recherche du