Aurélien Royer
email: aurelien.royer@u-bourgogne.fr

Blanca A García Yelo

Rémi Laffont

Manuel Hernández Fernández

Blanca A García

New bioclimatic models for the quaternary palaearctic based on insectivore

Keywords: Small mammalian communities, Paleoclimatology, Paleoecology, Pleistocene, Climate reconstructions

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

Understanding past climatic and environmental changes and their impact on continental ecosystems remains a major challenge, particularly in terms of human biological and social evolution. Currently, environmental reconstructions are built from mineralogical, geochemical, floral, and faunal data recovered from palaeontological or archaeological sites (e.g. [START_REF] Carrasco | Aplicaciones paleoclimáticas y paleoambientales de los estudios mineralógicos al yacimiento de vertebrados miocenos de Somosaguas[END_REF][START_REF] Frouin | Clay minerals of late Pleistocene sites (Jonzac and Les Cottés, SW France): Applications of X-ray diffraction analyses to local paleoclimatic and paleoenvironmental reconstructions[END_REF][START_REF] Domingo | Late Quaternary climatic and environmental conditions of northern Spain: An isotopic approach based on the mammalian record from La Paloma cave[END_REF][START_REF] Menéndez | Body-size structure of Central Iberian mammal fauna reveals semidesertic conditions during the middle Miocene Global Cooling Event[END_REF][START_REF] Villa | The Early Pleistocene herpetofauna of Rivoli Veronese (Northern Italy) as evidence for humid and forested glacial phases in the Gelasian of Southern Alps[END_REF][START_REF] Britton | Oxygen isotope analyses of Equus teeth evidences early Eemian and early Weichselian palaeotemperatures at the Middle Palaeolithic site of Neumark-Nord 2, Saxony-Anhalt, Germany[END_REF][START_REF] Girard | Fossil mega and micro-flora from Bernasso (Early Pleistocene, southern France): A multimethod comparative approach for paleoclimatic reconstruction[END_REF]. Vertebrates are particularly sensitive to changes in climate and habitat (e.g., Faunmap, 1996;[START_REF] Thomas | Extinction risk from climate change[END_REF], and it is generally assumed that these changes induce new biotic conditions, which, in turn, constrain vertebrate species distributions, biotic interactions and, ultimately, community organization [START_REF] Lyons | A Quantitative Assessment of the Range Shifts of Pleistocene Mammals[END_REF][START_REF] Moritz | Impact of a century of climate change on small-mammal communities in Yosemite National Park, USA[END_REF][START_REF] Hernández Fernández | Plio-Pleistocene climatic change had a major impact on the assembly and disassembly processes of Iberian rodent communities[END_REF][START_REF] Řičánková | Regional patterns of postglacial changes in the Palearctic mammalian diversity indicate retreat to Siberian steppes rather than extinction[END_REF][START_REF] Royer | Investigating the influence of climate changes on rodent communities at a regional-scale (MIS 1-3, Southwestern France)[END_REF][START_REF] Blanco | Differential responses of Miocene rodent metacommunities to global climatic changes were mediated by environmental context[END_REF]. Based on this assumption, numerous methods have been developed to reconstruct environmental conditions and quantify past climate parameters from mammalian fossil remains (for syntheses see [START_REF] Artemiou | Mammalian community palaeoecology: a review of recent methods with special reference to Miocene mammalian faunas of Europe[END_REF][START_REF] Andrews | Palaeoecology and hominoid palaeoenvironments[END_REF][START_REF] Lyman | Paleoenvironmental reconstruction from faunal remains: ecological basics and analytical assumptions[END_REF][START_REF] Nieto | Inferencia paleoecológica en mamíferos cenozoicos: limitaciones metodológicas[END_REF]. Several such models are based on modern species distributions and ecological niches (e.g., [START_REF] Atkinson | Seasonal temperatures in Britain during the past 22,000 years, reconstructed using beetle remains[END_REF][START_REF] García | Evolución de la diversidad taxonómica de los micromamíferos en la península ibérica y cambios paleoambientales durante el Pleistoceno Superior[END_REF][START_REF] Jeannet | L'écologie quantifiée. Essai de description de l'environnement continental à l'aide des microvertébrés[END_REF][START_REF] Fagoaga | Improving the accuracy of small vertebrate-based palaeoclimatic reconstructions derived from the Mutual Ecogeographic Range. A case study using geographic information systems and UDA-ODA discrimination methodology[END_REF], while others employ species richness (Montuire et al., 1997;[START_REF] Legendre | Rodents and climate: a new model for estimating past temperatures[END_REF][START_REF] Araújo | Quaternary climate changes explain diversity among reptiles and amphibians[END_REF]. Relatively few, however, utilize whole mammal communities [START_REF] Valverde | Remarques sur la structure et l'évolution des communautés de Vertébrés terrestres I. structure d'une communauté II. Rapports entre prédateurs et proies[END_REF][START_REF] Legendre | Analysis of mammalian communities from the late Eocene and Oligocene of southern France[END_REF][START_REF] Legendre | Les communautés de mammifères du Paléogène (Eocène supérieur et Oligocène) d'Europe occidentale: structures, milieux et évolution[END_REF][START_REF] Kay | Mammals and rainfall: paleoecology of the middle Miocene at La Venta (Colombia, South America)[END_REF][START_REF] García Yelo | Palaeoenvironmental analysis of the Aragonian (middle Miocene) mammalian faunas from the Madrid Basin based on body-size structure[END_REF][START_REF] Escarguel | Unearthing deep-time biodiversity changes: the Palaeogene mammalian metacommunity of the Quercy and Limagne area (Massif Central, France)[END_REF] or large mammal species associations [START_REF] De Bonis | Multivariate study of late Cenozoic mammalian faunal compositions and paleoecology[END_REF][START_REF] Griggo | Etablissement de courbes climatiques quantifiées à partir des communautés animales pléistocènes suivi d'une application aux gisements de l'abri Suard (Charente) et de la grotte du Bois-Ragot (Vienne)[END_REF][START_REF] Hernández Fernández | Plio-Pleistocene climatic change in the Turkana Basin (East Africa): evidence from large mammal faunas[END_REF]. Most environmental reconstructions are, in fact, based on small mammals or uniquely rodents [START_REF] Hokr | A method of the quantitative determination of the climate in the Quaternary period by means of mammal associations[END_REF][START_REF] Van De Weerd | Quantitative composition of rodent faunas in the Spanish Neogene and paleoecological implications (I and II)[END_REF][START_REF] Daams | Paleoenvironmental and paleoclimatic interpretation of micromammal faunal successions in the Upper Oligocene and Miocene of North Central Spain[END_REF][START_REF] Chaline | Les Rongeurs: leurs significations paléoécologiques et paléoclimatiques in La Baume de Gigny (Jura)[END_REF][START_REF] Sesé | Interpretación paleoclimática delas faunas de micromamíferos del Mioceno, Plioceno y Pleistoceno de la Cuenca de Guadix-Baza (Granada, España)[END_REF][START_REF] Van Der Meulen | Evolution of Early-Middle Miocene rodent faunas in relation to long-term palaeoenvironmental changes[END_REF][START_REF] Van Dam | Geographic and temporal patterns in the late Neogene (12-3 Ma) aridification of Europe: the use of small mammals as paleoprecipitation proxies[END_REF], as these species present several unique particularities. First, they are frequently preserved in archaeological and palaeontological deposits in caves and rock shelters, sometimes in association with human remains and artifacts and, unlike larger mammals, are less susceptible to human-induced biases (e.g., [START_REF] Chaline | Les rongeurs au pilori? Analyse critique des reconstitutions paléoécologiques en microthériodontologie[END_REF][START_REF] Discamps | Reconstructing palaeoenvironmental conditions faced by Mousterian hunters during MIS 5 to 3 in southwestern France: A multi-scale approach using data from large and small mammal communities[END_REF]. Second, small mammals are key prey for numerous predators and inhabit a wide variety of ecological habitats (e.g., [START_REF] Quéré | Les rongeurs de France: Faunistique et biologie. 3e édition revue et augmentée[END_REF]Wilson et al., 2017Wilson et al., , 2018)). Finally, due to their small size and large taxonomic diversity, rodent communities reflect local environments and are highly sensitive to environmental parameters and climatic changes, leading to their being considered one of the most reliable proxies for inferring past environmental and climatic conditions (e.g., [START_REF] Hinton | Monograph of the voles and lemmings[END_REF][START_REF] Tchernov | Rodent faunas and environmental changes in the Pleistocene of Israel[END_REF][START_REF] Avery | Late pleistocene coastal environment of the Southern Cape province of South Africa: Micromammals from Klasies river Mouth[END_REF][START_REF] Montuire | Palaeoecological analysis of mammalian faunas and environmental evolution in the South of France during the Pleistocene[END_REF][START_REF] Hernández Fernández | Bioclimatic analysis of rodent palaeofaunas reveals severe climatic changes in Southwestern Europe during the Plio-Pleistocene[END_REF][START_REF] García-Alix | Late Miocene-Early Pliocene climatic evolution of the Granada Basin (southern Spain) deduced from the paleoecology of the micromammal associations[END_REF][START_REF] Cuenca-Bescós | The reconstruction of past environments through small mammals: from the Mousterian to the Bronze Age in El Mirón Cave (Cantabria, Spain)[END_REF][START_REF] López-García | Palaeoenvironmental and palaeoclimatic reconstruction of the Latest Pleistocene of El Portalón Site, Sierra de Atapuerca, northwestern Spain[END_REF][START_REF] Cuenca-Bescós | The Early-Middle Pleistocene environmental and climatic change and the human expansion in Western Europe: a case study with small vertebrates (Gran Dolina, Atapuerca, Spain)[END_REF][START_REF] Royer | Late Pleistocene (MIS 3-4) climate inferred from micromammal communities and δ 18 O of rodents from Les Pradelles, France[END_REF][START_REF] Royer | Summer air temperature, reconstructions from the last glacial stage based on rodents from the site Taillis-des-Coteaux (Vienne), Western France[END_REF][START_REF] Rofes | Combining small-vertebrate, marine and stable-isotope data to reconstruct past environments[END_REF][START_REF] Royer | Investigating the influence of climate changes on rodent communities at a regional-scale (MIS 1-3, Southwestern France)[END_REF][START_REF] Fernández-García | Palaeoecological implications of rodents as proxies for the Late Pleistocene-Holocene environmental and climatic changes in northeastern Iberia[END_REF][START_REF] Laplana | Los micromamíferos y la interpretación ambiental de los yacimientos de la Cuenca de Villarroya[END_REF][START_REF] Piñero | Paleoenvironmental reconstruction of the Early Pleistocene site of Quibas (SE Spain) using a rodent assemblage reconstruction[END_REF][START_REF] Berto | Paleoenvironmental and paleoclimatic context during the Upper Palaeolithic (late Upper Pleistocene) in the Italian Peninsula. The small mammal record from Grotta Paglicci (Rignano Garganico, Foggia, Southern Italy)[END_REF][START_REF] López-García | Middle to Late Pleistocene environmental and climatic reconstruction of the human occurrence at Grotta Maggiore di San Bernardino (Vicenza, Italy) through the small-mammal assemblage[END_REF][START_REF] Berto | Palaeoenvironment and palaeoclimate in the western Liguria region (northwestern Italy) during the Last Glacial. The small mammal sequence of Riparo Mochi (Balzi Rossi, Ventimiglia)[END_REF][START_REF] Stoetzel | Quaternary terrestrial microvertebrates from mediterranean northwestern Africa: State-of-the-art focused on recent multidisciplinary studies[END_REF][START_REF] Fernández-García | Combined palaeoecological methods using small-mammal assemblages to decipher environmental context of a long-term Neanderthal settlement in northeastern Iberia[END_REF].

Among the many methods for reconstructing palaeonvironments, Bioclimatic Analysis [START_REF] Hernández Fernández | Bioclimatic discriminant capacity of terrestrial mammal faunas[END_REF] produces both qualitative and quantitative climatic reconstructions built from the structure of mammal communities, including multiple taxonomic orders, and their global distribution along broad climatic gradients and biozones. One of its main advantages is that entire mammal communities are taken into consideration, giving equal value to all species as climatic indicators related to their ecological tolerance. Here we further refine this approach by using Eulipotyphla (moles, shrews and hedgehogs, commonly referred to simply as insectivores) combined with Rodentia (rodents), or uniquely Rodentia recovered from Pleistocene to Late Holocene Palaearctic faunal assemblages. Although recent research has begun to combine rodents and insectivores in analyzes, these small mammals had almost always been previously considered separately by specialists. Fossil material from palaeontological and archaeological sites commonly result from the same accumulation agent and are recovered using similar excavation methods.

Moreover, Eulipotyphla are regularly preserved in European fossil sites and are highly sensitive to several climatic parameters, such as rainfall, (e.g., [START_REF] Reumer | The evolution of shrews (Mammalia, Soricidae): a key role for humidity[END_REF][START_REF] Hernández Fernández | Bioclimatic discriminant capacity of terrestrial mammal faunas[END_REF][START_REF] Héran | Mammifères et climat : reconstitutions paléoclimatiques du Cénozoïque d'Europe occidentale sur la base des faunes fossiles[END_REF][START_REF] Van Dam | Geographic and temporal patterns in the late Neogene (12-3 Ma) aridification of Europe: the use of small mammals as paleoprecipitation proxies[END_REF]Furio et al., 2011), making them considerably more robust proxies for inferring past climates. Furthermore, integrating these taxa is equally essential for generating a broader, more holistic vision of the past ecology of small mammal communities. Quaternary sites in the Palaearctic realm frequently yield both human and mammalian fossil remains (e.g., [START_REF] Puzachenko | Evolution of mammal species composition and species richness during the Late Pleistocene-Holocene transition in Europe: A general view at the regional scale[END_REF] and provide essential evidence for reconstructing the influence of climate changes on human biological and cultural evolution. Developing new palaeoclimatic modelling techniques focused on this realms's rich fossil record potentially offers a means to increase both the precision and spatiotemporal resolution of palaeoclimatic inferences. To test these new models, we mobilized data from Balma de l'Abeurador (Languedoc-Roussillon, France) and El Mirón (Cantabria, Spain), two archaeological sequences that cover the period from the Last Glacial Maximum to the beginning of the Holocene.

Material and methods

The Bioclimatic Analysis technique originally described by Hernández [START_REF] Hernández Fernández | Bioclimatic discriminant capacity of terrestrial mammal faunas[END_REF] and Hernández Fernández andPeláez-Campomanes (2003, 2005) consists of two sets of analyses.

The first uses linear discriminant functions deduced from the bioclimatic structure of modern mammalian communities in a specific climate zone (defined following Walter's typology, 1970), which are subsequently used to classify additional observations (extinct communities in our case) in each climatic zones. The second is built from transfer functions by means of multiple linear regression analyses of climatic parameters and modern bioclimatic mammalian community structures. These models are ultimately used to infer climatic variables for additional data (i.e. extinct communities). These two sets of analyses use a single global dataset of localities for which the composition of mammalian communities and climatic variables are known.

-Localities, faunal lists, and climate variables

While multiple global bioclimatic classification systems (e.g. Rivas-Martinez, 2008) are available, we employed the climate typology defined by [START_REF] Walter | Vegetationszonen und Klima[END_REF], which is based in the interaction between monthly precipitation and average temperatures. This typology was originally used in the original qualitative Bioclimatic Analysis [START_REF] Hernández Fernández | Bioclimatic discriminant capacity of terrestrial mammal faunas[END_REF], due to its simplicity and correlations with major terrestrial biomes. We collated faunal data from 49 modern localities from across the Palearctic (Table 1, Fig. 1), seven for each of the climate zones present in this biogeographic realm (see Hernández Fernández, 2001 for more details). Localities were selected to be representative of the average climatic conditions of each climate zone and to be as widely distributed as possible (see Table 1 for details). Our sample includes 18 localities from Hernández Fernández (2001), 13 from García Yelo et al. (2014), to which we added an additional 18 new localities (Table 1), for a total of 49 data points. Localities in climate zones I (equatorial climate), II (tropical climate with summer rains), and II/III (transitional tropical to semi-arid climate) are not present in the Palaearctic realm and, therefore, contrary to the methodology used by Hernandez Fernandez (2001), were excluded in the development of the new models.

In order to compare results from the two approaches, we built two faunal lists from published data [START_REF] Rutilevskiy | The zoogeography of the Yamal-Gydan area, Western Siberia[END_REF][START_REF] Hernández Fernández | Bioclimatic discriminant capacity of terrestrial mammal faunas[END_REF][START_REF] Héran | Mammifères et climat : reconstitutions paléoclimatiques du Cénozoïque d'Europe occidentale sur la base des faunes fossiles[END_REF][START_REF] García Yelo | Palaeoenvironmental analysis of the Aragonian (middle Miocene) mammalian faunas from the Madrid Basin based on body-size structure[END_REF] for each locality: one comprising only rodent species, the other including both rodents and insectivores. Species taxonomies were revised for each list, integrating the most recent species distribution data (Wilson et al., 2017;Wilson and Mittermeier, 2018;IUCN red list, 2019).

Nine climatic variables were collected for each locality (Table 1): mean annual temperature (MAT); annual positive temperature (TP), calculated as the sum of monthly mean temperature for months with average temperatures above 0 °C; mean temperature of the warmest month (Tmax); mean temperature of the coldest month (Tmin); mean annual thermal amplitude (MTA), calculated as the difference between Tmax and Tmin; thermicity index (It), which reflects the intensity of winter cold (eq .1); It = 10*(MAT +2xTmin) (eq.1)

compensated thermicity index (Itc), which was designed to balance the harsh winter in the continental climates and the extremely mild winter in the oceanic regions in order to make this index comparable worldwide (see Rivas [START_REF] Rivas-Martínez | Clasificación bioclimática de la Tierra[END_REF]; annual total precipitation (P); drought length (D), which is the period when P values are lower than 2*MAT.

2.2-Calculation of bioclimatic spectra

A climate zone vs species presence matrix was constructed for each locality. When a species is absent from a given climate zone, it is encoded as 0. When the species is present in that climate zone, a value of 1/n (referred to as the Climatic Restriction Index, CRI; Hernández [START_REF] Hernández Fernández | Bioclimatic discriminant capacity of terrestrial mammal faunas[END_REF], where n is the number of climates in which the species is present (the division by n ensures that the sum of CRI values for a given species equals 1). The more climatically restricted a species is, the higher its CRI values it has. As species with higher CRI imply more specific environmental requirements, they offer more reliable indications of climatic conditions.

CRI values for Rodentia and Eulipotyphla species were reexamined (Supplementary material 1) and modified following updated taxonomic and geographical range information (Wilson et al., 2017;Wilson and Mittermeier, 2018). Although localities from the three tropical climatic zones were not included, the large range distributions of some small mammal species extend into these zones. Consequently, all ten global climate zones were used to calculate the different Bioclimatic Components (BC) per locality (see below).

For a given locality, the Bioclimatic Component (BC) is the representation of each of the 10 existing climates. All 10 BC values for each locality were calculated according to the formula (eq. 2):

𝐵𝐵𝐵𝐵 𝑖𝑖 = �100 ∑ 𝐵𝐵𝐶𝐶𝐶𝐶 𝑖𝑖𝑖𝑖 𝑘𝑘 𝑖𝑖=1 � 𝑆𝑆 ⁄ (eq. 2)
where i stands for the i th climate zone, j for the j th species present in the locality, with S being the number of species in the locality. The 10 BC values for a given locality constitute its bioclimatic spectrum. The bioclimatic spectrum of each locality was calculated using both datasets: one comprised uniquely of rodents (Supplementary material 2), the other including both Rodentia and Eulipotyphla species (Supplementary material 3).

2.3-Palaeoclimatic inferences

All statistical analyses were performed using the R software package (R core team, 2019).

Linear Discriminant Analysis (LDA - [START_REF] Sneath | Numerical Taxonomy -Freeman and Company[END_REF][START_REF] Rencher | Methods of Multivariate Analysis -2 nd edition[END_REF][START_REF] Legendre | Numerical Ecology -3 rd edition[END_REF] with the LDA function from the R package MASS [START_REF] Venables | Modern Applied Statistics with S-PLUS[END_REF] was used to classify localities into different climate zones. We subsequently used a Leave-One-Out Cross Validation (LOOCV - [START_REF] Efron | The Jackknife, the Bootstrap and Other Resampling Plans -Society for Industrial and Applied Mathematics[END_REF][START_REF] Chernick | Bootstrap Methods: A Practioner's Guide[END_REF] to evaluate classification quality.

Predictive equations for each of climatic variable were obtained from multiple linear regressions [START_REF] Sokal | Biometry -3 rd edition[END_REF][START_REF] Zar | Biostatistical Analysis -4 th edition[END_REF][START_REF] Rencher | Methods of Multivariate Analysis -2 nd edition[END_REF][START_REF] Legendre | Numerical Ecology -3 rd edition[END_REF] using the lm() function, with the ten bioclimatic components as independent variables and the climatic variable as the dependent ones. We generated a script to easily calculate the bioclimatic spectra and inferred climatic zone, which is available as supplementary material (Supplementary material 4 and available at https://github.com/AurelienRoyer/PalBER).

2.4-Testing sampling biases on palaeoclimatic inference through randomization

Small mammal associations from fossil sites potentially present biased representations of past community structure due the nature of the accumulation and subsequent taphonomic alterations [START_REF] Clark | Oligocene sedimentation, stratigraphy, paleoecology and paleoclimatology: in the Big Badlands of South Dakota[END_REF][START_REF] Denys | Nouveaux critères de reconnaissance des concentrations de microvertébrés d'après l'étude des pelotes de chouettes du Botswana (Afrique australe)[END_REF][START_REF] Andrews | Owls, Caves and Fossils[END_REF][START_REF] Kowalski | Some problems of the taphonomy of small mammals[END_REF][START_REF] Stahl | The recovery and interpretation of microvertebrate bone assemblages from archaeological contexts[END_REF][START_REF] Reed | Serengeti micromammals and their implications for Olduvai paleoenvironments[END_REF], which may, in turn, impact climatic inferences drawn from Bioclimatic Analysis. For instance, many fossil assemblages of small mammals are accumulated by predators, especially nocturnal raptors, who are relatively selective in their choice of prey (e.g., [START_REF] Andrews | Owls, Caves and Fossils[END_REF]. Bioclimatic Analysis has the advantage of dampening this potential bias by taking into account the general representativeness of the taxonomic composition of micromammal death assemblages (i.e. the presence or absence of particular taxa) rather than their relative abundance.

This method does not differentially weight taxa, and the occurrence of a single fossil of a unique species could theoretically alter a locality's bioclimatic spectrum and, therefore, its climate reconstruction. However, although high levels of ecological redundancy in mammalian communities [START_REF] Cano | Cenogramas, análisis bioclimático y muestreo en faunas de mamíferos: implicaciones para la aplicación de métodos de análisis paleoecológico[END_REF] somewhat minimize this problem, it can be better overcome with adequate sample size. Although the fragile nature of small mammal fossils means they are unlikely to be redistributed [START_REF] Sierra | Taphonomic and sedimentary factors in the fossil record of mammals. European neogene mammal chronology[END_REF][START_REF] Fernández-López | Temas de Tafonomía[END_REF], potential species loss associated with different processes of accumulation, fossilization, and sampling bias cannot be ruled out [START_REF] Andrews | Owls, Caves and Fossils[END_REF][START_REF] Behrensmeyer | Terrestrial vertebrate accumulations. Taphonomy: releasing the data locked in the fossil record[END_REF][START_REF] Lyman | Vertebrate taphonomy[END_REF][START_REF] Andrews | Taphonomic effects of faunal impoverishment and faunal mixing[END_REF][START_REF] Cano | Cenogramas, análisis bioclimático y muestreo en faunas de mamíferos: implicaciones para la aplicación de métodos de análisis paleoecológico[END_REF][START_REF] Turvey | The past is another country: is evidence for prehistoric, historical and present-day extinction really comparable[END_REF][START_REF] Fernández-Jalvo | Taphonomy in palaeoecological interpretations[END_REF], raising questions as to how these potential biases affect Bioclimatic Analysis (Gómez [START_REF] Cano | Cenogramas, análisis bioclimático y muestreo en faunas de mamíferos: implicaciones para la aplicación de métodos de análisis paleoecológico[END_REF]).

To investigate the potential impact of species loss, we randomly removed species in the modern localities in order to replicate impoverished faunal associations. Multiple LDA and multiple linear regressions were performed for each locality by randomly removing k taxa (k varying from 1 to n-2, with n being number of species in the locality). For a given locality and a

given k, 100 impoverished samples were generated. A total of 69,000 randomized samples were produced, for which simulated climatic inferences were calculated using the new rodents and rodents plus insectivores models. Three locations in the arctic climate zone (IX) were left unchanged, as they contain less than three species. We used logistic regressions [START_REF] Sokal | Biometry -3 rd edition[END_REF][START_REF] Legendre | Numerical Ecology -3 rd edition[END_REF] to explore the coherence of simulated climatic inferences for each climate zone. Quantitative inferences of climatic variables were evaluated using standard deviations of predicted values calculated for a given locality and a given k taxa removed.

ANCOVAs [START_REF] Zar | Biostatistical Analysis -4 th edition[END_REF][START_REF] Sokal | Biometry -3 rd edition[END_REF] were performed for each climate parameter using log-transformed standard deviations of predicted values, with the percentage of taxon removed and climate zones used as independent variables.

2.5-Application of the new models to palaeoclimatic reconstruction

The two new approaches for Bioclimatic Analysis allow us to analyze small mammal associations from the Palaearctic Quaternary and reconstruct climatic conditions. We used small mammal associations from La Balma de l'Abeurador (Languedoc-Roussillon, Hérault, France)

and El Mirón (Cantabria, Spain), two sites located in typical temperate climate areas and that yielded sequences spanning from the end of the Late Pleistocene to the beginning of the Holocene (Fig. 2). All radiocarbon dates were calibrated at 2σ (95%) using Oxcal v4.3 (Bronk Ramsey, 2017) and the IntCal13 calibration curve [START_REF] Reimer | IntCal13 and Marine13 radiocarbon age calibration curves 0-50,000 Years cal BP[END_REF].

The site of Balma de l'Abeurador is located in southern France and currently has a mean annual temperature of 13.6°C and annual precipitation around 700 mm. The central sector of the cave produced six multi-level deposits [START_REF] Heinz | La grotte de l'Abeurador (Félines-Minervois, Hérault): paléoécologie d'après l'anthracologie et la carpologie[END_REF][START_REF] Vaquer | Mesolithic plant exploitation at the Balma Abeurador (France)[END_REF][START_REF] Vaquer | La grotte de l'Abeurador Félines-Minervois (Hérault) : occupations humaines et environnement du Tardiglaciaire à l'Holocène. De Méditerranée et d'ailleurs…[END_REF], whose chronology was determined by the nature of the lithic assemblage, palynology, and radiocarbon dates [START_REF] Vaquer | La grotte de l'Abeurador Félines-Minervois (Hérault) : occupations humaines et environnement du Tardiglaciaire à l'Holocène. De Méditerranée et d'ailleurs…[END_REF]. Small mammals from Balma de l'Abeurador were first studied by [START_REF] Marquet | Les modifications de l'environnement Postglaciaire en France méditerranéenne d'après les rongeurs de l'Abeurador et de Font-Juvenal. Premières communautés paysannes en méditerranée occidentale[END_REF] and then by [START_REF] Mistrot | Contribution des micromammifères de la Balma de l'Abeurador à la connaissance de l'évolution des paysages tardiglaciaires et holocènes en Languedoc-Roussillon[END_REF], who identified at least 19 taxa in an assemblage comprising at least 2000 individuals (Table 2). We used [START_REF] Mistrot | Contribution des micromammifères de la Balma de l'Abeurador à la connaissance de l'évolution des paysages tardiglaciaires et holocènes en Languedoc-Roussillon[END_REF] data, as it benefits from a taphonomical analysis demonstrating no bias in the representativity of the small mammal assemblage. This detailed taphonomical approach, together with the high number of remains (more than 100 first and second upper and lower rodent molars for each level) and an adequate minimum fossil sample size of fossils, allows for reliable palaeoclimatic inferences (van de Weerd and [START_REF] Van De Weerd | Quantitative composition of rodent faunas in the Spanish Neogene and paleoecological implications (I and II)[END_REF][START_REF] Daams | Stratigraphy and sedimentology of the Aragonian (Early to Middle Miocene) in its type area (North-Central Spain)[END_REF]. Several differences do exist between [START_REF] Mistrot | Contribution des micromammifères de la Balma de l'Abeurador à la connaissance de l'évolution des paysages tardiglaciaires et holocènes en Languedoc-Roussillon[END_REF] and [START_REF] Marquet | Les modifications de l'environnement Postglaciaire en France méditerranéenne d'après les rongeurs de l'Abeurador et de Font-Juvenal. Premières communautés paysannes en méditerranée occidentale[END_REF] analyses in terms of species identification, with the latter employing extreme caution when identifying fossils to species. For Apodemus, Arvicola, and Microtus (Terricola), [START_REF] Mistrot | Contribution des micromammifères de la Balma de l'Abeurador à la connaissance de l'évolution des paysages tardiglaciaires et holocènes en Languedoc-Roussillon[END_REF] made no specific identifications, whereas [START_REF] Marquet | Les modifications de l'environnement Postglaciaire en France méditerranéenne d'après les rongeurs de l'Abeurador et de Font-Juvenal. Premières communautés paysannes en méditerranée occidentale[END_REF] Crosses indicate species present, circles represent absence species. Data are from [START_REF] Mistrot | Contribution des micromammifères de la Balma de l'Abeurador à la connaissance de l'évolution des paysages tardiglaciaires et holocènes en Languedoc-Roussillon[END_REF].

The site of El Mirón is located in northern Spain, and enjoys relatively similar climate conditions to Balma de l'Abeurador, but with slightly higher annual precipitation (reaching 1,000 mm) and smaller seasonal differences in rainfall (Fig. 2). The site yielded a long stratigraphic sequence (Strauss and Morales, 2019) and numerous radiocarbon dates (e.g., Strauss and Morales, 2003[START_REF] Reed | Serengeti micromammals and their implications for Olduvai paleoenvironments[END_REF], 2012[START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF]. Twenty-one insectivore and rodent taxa were identified from throughout the 44 studied levels of the stratigraphy [START_REF] Cuenca-Bescós | Paleoclima y paisaje del final del Cuaternario en Cantabria: los pequeños mamíferos de la Cueva del Mirón (Ramales de la victoria)[END_REF], 2009). Although over 3,000 individuals were identified, they are highly unevenly distributed among the different levels. To ensure the reliability of palaeoclimatic inferences, we selected only deposits with more than 50 individuals, or grouped together several successive levels (i.e. the Solutrean levels 120 to 126 from the Last Glacial Maximum). For Apodemus and Sorex, no distinctions were made between Ap.

sylvaticus and Ap. flavicollis or between S. araneus and S. coronatus.

When reliable criteria for confidently distinguishing between species are unavailable, two approaches can be taken: either building bioclimatic spectra from all the potentially present species and calculating a palaeoclimate inferences for each; or using an average species, sometimes referred to as "quimeric" species. When a taxon is only identified to genus, the biome occupation of all the potential species can be used. For instance, if one species inhabits biomes VI and VII, and the other inhabits IV and VI; the genus can be considered to inhabit IV, VI, and VII. Even if relatively little uncertainty is introduced, "quimeric" species have a lower discriminatory power than any of the potential original species while still providing useful information for community structure. Consequently, both approaches were used and their results compared.

Extinct species potentially present an additional limitation, as no past climate data is available for their geographic distribution. Following Hernández Fernández et al. ( 2007), we assigned extinct taxa to different BCs according to the distribution of their extant relatives. As "quimeric" species, they provide reduced discriminatory power. The El Mirón sequence has the particularity of yielding remains of the now-extinct rodent Pliomys coronensis (=Pliomys lenki),

for which we used the bioclimate characterization proposed by Hernández Fernández et al. ( 2007).

3-Results and discussion

3.1-Qualitative climatic classification and quantitative climatic inferences

The number of species in modern localities from the different climate zones ranges on average between 17 and 24, except for bioclimate zones III and IX (Supplementary material 6).

Localities in these two climate zones yielded a maximum of 15 species, with some containing less than 5 species. Eulipotyphla produced around 25 % of all faunal lists on average, but faunal lists from bioclimate zones IV, VI and VIII included 30 % or more of Eulipotyphla.

The two linear discriminant analyses based on the bioclimatic spectra correctly classified nearly all localities, inducing a very low rate of error. When only Rodentia was used, incorrectly classified localities rose to 6 % ( with a posterior probability > 0.9) for each approach (Rodentia, Rodentia plus Eulipotyphla) for LDA. The prediction error is based on LOOCV.

The Leave-One-Out Cross Validation (LOOCV) indicated very low prediction error; the rodent-only dataset produced approximately 12 % misclassifications, while the dataset including both sets of mammals had both a lower level of error (Table 3) and predictions associated with a high posterior probability (Table 3). This would indicate climate inferences based on the combination of rodents and insectivores to be slightly more reliable. The high predictive power of both models reveals a close relationship between small mammal communities and climate, which makes them suitable to infer past climate conditions [START_REF] Hernández Fernández | Bioclimatic discriminant capacity of terrestrial mammal faunas[END_REF].

Most predictive equations generated by the multiple linear regressions for each climatic factor, as a function of the bioclimatic components of either Rodentia or Rodentia and Eulipotyphla, produced highly significant determination coefficients (Table 4). The different transfer functions for each climatic variable had similar coefficients of determination for both datasets, with minor trends towards slightly higher values for the combined Rodentia and Eulipotyphla dataset. Values were also relatively similar to those obtained by Hernández

Fernández and Peláez-Campomanes (2005). Only MTA showed very low coefficients of determination (<0.6). Precipitation shows a slightly better coefficient of determination than the one obtained by Hernández Fernández and Peláez-Campomanes ( 2005), although it is still lower than 0.8 (Table 4), which resulted in calculated prediction with large intervals.

3.2-Loss of species and bias in climate inferences

Overall, few incorrect qualitative inferences were generated when species were removed from the faunal lists (usually below 15 % per climate zone -fig. 3, supplementary material 7). The S-shaped curves of logistic regressions are skewed upwards (Fig. 3; Supplementary material 7), meaning there is no "turning point" when climatic inferences become erroneous and that reliable classifications can nevertheless be obtained using relatively few species, although they may have a high probability of error (Fig. 4). The two approaches behaved differently according to climate zone when species are randomly removed; inferences drawn from rodent communities are more robust for climate zones III, V, and IX, while the second approach (Rodentia plus Eulipotyphla) is more robust for climate zones IV, VI, and VIII. Incorrect inferences for communities in the polar climate zone (IX) are low with the rodent-only dataset (4.2%) and much higher with the Rodentia plus Eulipotyphla, with incorrect inferences mainly classified into the boreal climate zone (VIII).

The percentage of incorrect classifications increased for all climate zones when more taxa are removed (Figs. 3 and4, supplementary material 7). Although classification error generally increased after the removal of 30-40 % of taxa from the localities (Fig. 4), most reached the 25 % incorrect classifications threshold when more than 50 % of taxa were removed. Nevertheless, the warm temperate (climate zone V; except for the Batumi locality) and arid temperate climate (VII) zones appeared to be much more resistant to species loss independent of the dataset used (Fig. 3, supplementary material 7). Substantial variability is evident between localities in the same climate zone, demonstrating faunal communities from the same climate zone to be highly variable. This variability reflects varying responses of different taxa inhabiting the same climate zone.

Quantification of climatic parameters follows similar trends. Differences between theoretical and calculated values in impoverished communities for each climatic parameter are highly dependent on the percentage of remaining taxa (Supplementary material 8). These differences vary as a function of the climate zone considered and the dataset used. For instance, in the case of precipitation estimates from climate zones III and VIII (subtropical arid and boreal), less impoverished communities produced larger differences when the rodent-only dataset was used. Figures 5 and6 show standard deviations of predicted values obtained for impoverished faunal associations computed for each climatic parameters. These standard deviations follow exponential variations, which increase with the percentage of taxa removed. Overall, standard deviations of predicted values become important when half of the faunal list is removed (Fig. 5 and6). ANCOVAs for each climate parameter show no significant difference between climate zones (F = (6, 676), p > 0.7). The approach based on rodent communities generated steeper regression curves than those based on Rodentia plus Eulipotyphla, which suggests the former to be more sensitive to species loss. This is probably related to the fact that faunal lists comprised exclusively of rodents contain less species than those combining both orders of small mammals.

When confronted with species loss, the combined dataset therefore appears better adapted to producing robust climate quantifications.

As expected, our results indicate that palaeoclimatic inferences obtained by Bioclimatic Analysis depend significantly on the integrity of the past faunal communities used. Satisfactory levels of confidence can be obtained only when species loss is below 50 %, which is similar to the initial estimates produced by Gómez [START_REF] Cano | Cenogramas, análisis bioclimático y muestreo en faunas de mamíferos: implicaciones para la aplicación de métodos de análisis paleoecológico[END_REF]. When quantifying climatic parameters, the approach combining Rodentia and Eulipotyphla proved less sensitive to faunal impoverishment, and should be preferred when species loss is suspected.

4-Palaeoclimatic inferences from past small mammal communities

Climate zone classifications and mean annual temperature predictions generated by the two approaches were applied to multiple faunal lists from El Mirón and Abeurador (Figure 7). All possibilities related to uncertainties in species identifications (see above) were tested. Climate prediction remained the same for both datasets, except for the case of Abeurador level 5b, for which two different climate zones (IV and VI) were predicted depending on the Eulipotyphla species used. When uncertainties concern only two species, variations in mean annual temperature values remain minimal (Fig. 7, Supplementary material 5). For instance, temperature variation was less than 0.5°C between communities in which Apodemus specimens were identified as A.

sylvaticus instead of A. flavicollis, and vice versa. In the case of N. anomalus and N. fodiens, variation reached almost 2°C. Observed variations can reach up to 4°C (e.g. Abeurador level 8) when all possibilities are run using either of the two approaches. The highest degree of variation primarily concerns uncertain species inhabiting only a few climate zone that are different between them, such as S. araneus and S. coronatus. As a consequence, faunal lists employed in bioclimatic modelling should only include reliably identified species. Finally, small variations could also be observed according to the approach used, the one based on rodents tending to estimate mean annual temperature two degrees lower.

Figure 8 depicts the mean annual temperature obtained from small mammal lists of Abeurador and El Mirón. For levels with uncertainties in species classification, we used "quimeric" species for climatic inferences. Overall, the two approaches show similar trends in both sequences; an attribution to climate zone IV during the Holocene and a progressive warming with the onset of the Bølling interstadial. While the same trend can be observed for predictions based on Rodentia and Eulipotyphla and only Rodentia, the latter tended to produce lower mean annual temperature estimates. Mean annual temperatures obtained for the Abeurador sequence using the rodent dataset ranges from 3.8°C to 9.2°C, whereas the approach based on Rodentia and Eulipotyphla ranges from 6.3 to 11.6°C. Similarly, mean annual temperature varied throughout the El Mirón sequence, averaging around 1.2°C for both approaches, which is included in the confidence intervals (and so in the prediction intervals). The lowest mean annual estimations were for the LGM levels 112-114. Temperatures rose around 3.0°C between Heinrich stadial 1 and the beginning of the Holocene at El Mirón, which is slightly lower than was calculated for Abeurador.

No major differences are observed between the two approaches for precipitation and thermal index conditions (Fig. 9, Supplementary material 5). Throughout the Abeurador sequence, the Eulipotyphla approach showed an increase of precipitation between Heinrich Stadial 1 and the Boreal interstadial, suggesting wetter climate conditions to have prevailed during the Holocene at

Abeurador. Nevertheless, no differences were inferred with the rodent approach, suggesting that if changes in precipitation did occur, they were not important enough to have a significant impact on rodent communities. A progressive increase in the Thermicity index can be observed throughout the sequence, confirming milder conditions associated with the Holocene deposits.

Apart from a slight decrease at the end of the sequence, no abrupt changes in precipitation were recorded anywhere in the El Mirón record, and no differences are evident between predictions for both approaches. Several variations in the Thermicity index are evident across the sequence, with several notable decreases in levels 112-114 and 107-106, whereas a slight increasing trends emerges at the end of the Late Glacial. These variations were recorded for a level where low mean annual temperature were predicted, suggesting these changes were associated with modifications in seasonality and the intensity of winter colds. Climate parameters from levels 110 and 111 of El Mirón were previously explored using the Mutual Climatic Range of the rodents [START_REF] Bañuls-Cardona | The end of the Last Glacial Maximum in the Iberian Peninsula characterized by the small-mammal assemblages[END_REF]. Mean annual temperatures were estimated to be 3 to 4 °C higher, averaging around 10.6°C, whereas the mean temperature of the coldest months diverge by 6°C. The Bioclimatic approach suggests colder climatic conditions for the LGM levels of El Mirón.

Our analysis of the Balma de l'Abeurador and El Mirón sequences primarily focused on three climate variables (mean annual temperature, precipitation, and thermicity index) considered susceptible to changes during this time period. These two terrestrial sequences are in relative agreement with evidence from regional climate evidence recorded in speleothems as well as lake and marine core records [START_REF] Jalut | La végétation des Pyrénées françaises lors du dernier épisode glaciaire et durant la transition Glaciaire-Interglaciaire (Last Termination). Els Pirineus i les àrees circumdants durant el Tardiglaciar. Mutacions i filiacions tecnoculturals, evolució paleoambiental (16 000[END_REF][START_REF] Moreno | Northern Iberian abrupt climate change dynamics during the last glacial cycle: a view from lacustrine sediments[END_REF]. The LGM appears as a cold, relatively dry period (but not as dry as Heinrich Stadial 1), inducing a hiatus in speleothem and lake records. Bioclimatic models built from LGM deposits tend to estimate lower mean annual temperature and a lower Thermicity index, as is the case in El Mirón level 112.114. With that said, our models revealed no major change in precipitation coincident with the LGM, with all variations falling within the large confidence intervals. Speleothems from northern Spain and southern France (Chauvet Cave) have suggested an increase in temperature and humidity just before the Bølling interstadial (15.4 kyr), which reached their maximum during the Allerød [START_REF] Genty | Timing and dynamics of the last deglaciation from European and North African δ 13 C stalagmite profiles-comparison with Chinese and South Hemisphere stalagmites[END_REF][START_REF] Moreno | A speleothem record of glacial (25-11.6 kyr BP) rapid climatic changes from northern Iberian Peninsula[END_REF][START_REF] Martínez Pillado | Upper Pleistocene and Holocene palaeoenvironmental records in Cueva Mayor karst (Atapuerca, Spain) from different proxies: speleothem crystal fabrics, palynology, and archaeology[END_REF]. Mean annual temperature based on small mammals closely mirror the mean annual temperature model built from speleothems. The Younger-Dryas at El Mirón is represented by two levels with radiocarbon dates compatible with this period. Overall, speleothem and pollen records suggest the Younger-Dryas to be characterized by an abrupt cold event associated with a reduction in arboreal communities [START_REF] Fletcher | Orbital-and sub-orbital-scale climate impacts on vegetation of the western Mediterranean basin over the last 48,000 yr[END_REF][START_REF] Morellón | Early Holocene humidity patterns in the Iberian Peninsula reconstructed from lake, pollen and speleothem records[END_REF]. There is currently little evidence supporting this event to comprise two phases; an initial dry phase followed by a progressive increase in humidity and warmer conditions (e.g. [START_REF] Muñoz-Sobrino | New data on the Lateglacial period of SW Europe: a high resolution multiproxy record from Laguna de la Roya (NW Iberia)[END_REF][START_REF] Bartolomé | Hydrological change in Southern Europe responding to increasing North Atlantic overturning during Greenland Stadial 1[END_REF][START_REF] Morellón | Early Holocene humidity patterns in the Iberian Peninsula reconstructed from lake, pollen and speleothem records[END_REF].

Two hypotheses could explain the pattern observed for the Younger-Dryas at El Mirón; low archaeological resolution or the relatively insignificant impact of this event on small mammal communities.

Despite being located on different sides of the Pyrenees, these two sequences are relatively close geographically and today lie in typical temperate climates. Our results suggest that warming during the Holocene differed slightly of either side of this natural barrier, being less marked in northern Spain compared to southern France.

Conclusions

By Table 2: Faunal lists of identified Rodentia and Eulipotyphla per level at Balma de l'Abeurador.

Crosses stand for species presence, circles for species absence. Data are from [START_REF] Mistrot | Contribution des micromammifères de la Balma de l'Abeurador à la connaissance de l'évolution des paysages tardiglaciaires et holocènes en Languedoc-Roussillon[END_REF].

Table 3: Apparent error rate, prediction error, and prediction quality (proportion of observations with a posterior probability > 0.9) for each approach (Rodentia, Rodentia plus Eulipotyphla) for LDA. The prediction error is based on LOOCV.

Figure 5 :

 5 Figure 5: Fitted curves (black) and associated confidence intervals (blue) for the exponential

Figure 6 :

 6 Figure 6: Fitted curves (black) for each climatic parameter and associated confidence intervals (in

Figure 7 :

 7 Figure 7: Influence of the species identification uncertainty on Mean Annual Temperature

Figure 8 :

 8 Figure 8: Evolution of the estimated Mean Annual Temperature along the El Mirón (bottom) and

Figure 9 :

 9 Figure 9: Evolution of estimated precipitation and thermicity index (TI) along the El Mirón

  localities added for this study are indicated with asterisks. Abbreviations: MAT = mean annual temperature; TP = annual positive temperature; Tmax = mean temperature of the warmest month; Tmin = mean temperature of the coldest month; MTA = mean annual thermal amplitude; It = thermicity index; Itc = compensated thermicity index; P = annual total precipitation; D= drought length.

Table 2 :

 2 duodecimcostatus throughout the sequence. For the insectivore genus Sorex, Neomys, and Crocidura, no distinction was made between Sorex araneus and Sorex coronatus, between Neomys anomalus and Neomys fodiens, or between Crocidura leucodon and Crocidura russula. Faunal lists of identified Rodentia and Eulipotyphla per level at Balma de l'Abeurador.

[START_REF] Brunet-Lecomte | Les Campagnols souterrains, Terricola, Arvicolidae, Rodentia, actuels et fossiles d'Europe occidentale[END_REF]

identified almost exclusively Ap. sylvaticus, Ar. Amphibius, and M.

(T.)

  Table 3), a similar figure to the one calculated by Hernández[START_REF] Hernández Fernández | Bioclimatic discriminant capacity of terrestrial mammal faunas[END_REF]. When the combined Rodentia and Eulipotyphla dataset is used, only Blagoveshchensk was incorrectly classified in climate zone VIII instead of VI. However, the cooltemperate and boreal forests throughout the region of Blagoveshchensk[START_REF] Kolbek | Forest Vegetation of Northeast Asia[END_REF] are home to numerous species typical of colder climates, including Lemmus amurensis,

Clethrionomys rufocanus, Myopus schisticolor, Pteromys volans, or Sorex tundrensis), which likely goes some way in explaining the incorrect climate zone attribution. Table 3: Apparent error rate, prediction error, and prediction quality (proportion of observations

  focusing on Palaearctic small mammal communities, our work refines the Bioclimatic Analysis first developed by Hernández Fernández two decades ago and extends its use to rich fossil record of the Eurasian Quaternary. Although rodent communities alone generate informative qualitative and quantitative climatic reconstructions, combining rodents and insectivores increases precision for several climatic parameters, particularly precipitation. As previously indicated by[START_REF] Hernández Fernández | Quantitative palaeoclimatic inference based on terrestrial mammal faunas[END_REF], due to the statistical uncertainty present in every climate transference function, the inferences drawn should not be taken as absolutes but as indicators of general trends. Although it is robust against sampling bias, this method nevertheless requires data from secure contexts with proper sampling efforts and adequate taphonomical control in order to ensure accurate taxonomic identification and the representativity of faunal assemblages. Although imprecisions in taxonomic identification do not impact general trends, they could engender variations in climatic estimations, which highlights the importance of confidence intervals. The R script developed for this analysis (and available as supplementary material 4) provides a simple means for using Bioclimatic Analysis for past climate reconstructions.The application of these new models to the Quaternary sequences of Balma de l'Abeurador and El Mirón to reconstruct paleoclimate revealed relatively similar changes with those built from other proxies, such as speleothems. The bioclimatic inferences from these two archaeological sequences suggest that, depending on the side of Pyrenees, mean annual temperatures warmed from 3 to 5°C between the Late Glacial and the Holocene.Figure1: Maps of the 10 climate zones and 49 selected localities. See Table1for the list of localities. Fitted curves (red) for logistic regressions of the correctness of the LDA classification for each climatic zone (vertical axes: above, correct predictions, below, incorrect predictions) against the percentage of remaining taxon (horizontal axes). This percentage was calculated from the randomization procedure of 49 rodent and Eulipotyphla faunal associations (seven for each climate zone), which were impoverished by randomly removing species. Histograms illustrate frequencies of classification as a function of the percentage of species remaining in the association (number of impoverished localities is shown on the right axis). Boxplots illustrate median, quartile, and range of incorrect inferences for all the randomly impoverished communities in each climate zone. A detailed indication of the percentage for each climate zone estimated by incorrect classifications is provided in the middle of each graph (numbers as in Table1).Figure4: Evolution (colored curves) of the percentage of incorrectly classified localities (vertical axes) depending on the percentages of remaining taxa (horizontal axes) obtained with randomly impoverished communities for each locality (numbers as in Table1) and per climate zone.Numbers of localities are reported in table 1. Blagoveshchensk was excluded from climate zone IV. Three localities(44, 45 and 46) were not randomly impoverished for climate zone IX due to very low species numbers.

	Figure 2: A) location of Abeurador and El Mirón, B) Walter and Lieth climate diagrams
	(calculated between 1985 and 2015) and, C) available radiocarbon dates for levels from both
	sequence with associated cultural attributions. Radiocarbon dates calibrated using OxCal v4.3 and
	IntCal13.
	Figure 3:

Table captions :

 captions Table 1: Coordinates, climate zone, and climatic parameters values for the 49 localities. New

Table 4 :

 4 Coefficient of determination values obtained for multiple linear regressions for each climatic variable according to the two approaches (Rodentia, Rodentia plus Eulipotyphla).

Acknowledgments

This research was conducted by the PMMV Team (Palaeoecology, Macroecology and Macroevolution of Vertebrates (http://pmmv.com.es), as part of Evolution of Cenozoic Mammals and Continental Palaeoenvironments research group (UCM-910607). This research was supported by projects of the Spanish Ministries of Education, Science and Innovation (PGC2018-094122-B-I00; PGC2018-094955-A-I00) and the Biogéosciences Laboratory of the University of Burgundy.

Supplementary captions:

Supplementary material 1: Climatic restriction index of Rodentia and Eulipotyphla. The last column provides details for each locality (numbers as in Table 1), from which these taxa were found.

Supplementary material 2: Bioclimatic spectra of Rodentia for each locality. Supplementary material 3: Bioclimatic spectra of Rodentia and Eulipotyphla for each locality. Supplementary material 4: Palber.R Script applied to bioclimatic method and bioclimatic spectra files (.csv) necessary to run it. Supplementary material 5: Predictions calculated for each level of Abeurador and El Mirón, obtained from the two approaches. Supplementary material 6: Number of rodents and Eulipotyphla per modern locality. Supplementary material 7: For each climatic zone, S-shape curves (red) for logistic regressions of the correctness of the LDA classification for each climate zone (vertical axes: above, correct predictions, below, wrong predictions) against the percentage of the remaining taxon (horizontal axes). This percentage was calculated from the randomization procedure of 49 rodent faunal associations (seven for each climate zone), which were impoverished by randomly removing species. Histograms illustrate frequencies of classification as a function of the percentage of species remaining in the association (number of impoverished localities is shown on the right axes). Boxplots illustrate median, quartile, and range of incorrect inferences for all the randomly impoverished communities in each climate zone. A detailed of incorrect classifications (%) are provided in the middle of each graph (numbers as in Table 1), with the average details of the climate zones into which they are classified. Supplementary material 8: Evolution of differences between estimated climate values by means of multiple regression for each climate variable and both approaches (Rodentia or Rodentia + Eulipotyphla) according to the percentage of remaining taxa and values calculated when the faunal list is complete. Each figure corresponds to a combination of a given climate variable and a given approach.