Paleoenvironmental and paleoclimatic reconstruction in the Western Mediterranean during the Late Early Pleistocene

The major climate change recorded during the Early-Middle Pleistocene Transition (1.4–0.4 Ma) is characterized by an increase in the length of climatic cycles from 41 to 100 ka. During this period, the Mediterranean climate underwent aridification associated with a drop in temperatures. This period also corresponds to the arrival of the first hominins in Western Europe. This study aims to establish the climatic framework at the end of the Early Pleistocene between MIS 37-31 (∼ 1.25 to 1.06 Ma). To this end, a multiproxy approach was applied to assemblages of planktonic foraminifera, pollen and clay mineralogy, constituting a multi-method approach to climate reconstruction. Comparisons with other sequences from the Central and Western Mediterranean show major differences between the Iberian Peninsula and southern Italy. Indeed, the MIS 31 appears cooler and dryer in the south of the Iberian Peninsula than in central Mediterranean. Clay mineralogy analysis highlights a sudden change in oceanic and atmospheric circulation in the Alboran Basin around 1140 ka. Climate reconstructions follow climatic cycles, with temperatures that appear to be underestimated compared to temperature reconstructions based on fauna from southern Spain. Precipitation reconstructions, on the other hand, are more consistent, suggesting that the climate and vegetation of this region are more likely to be influenced by variations in precipitation. These new data enhance our understanding of the climate of the South-western Mediterranean at the beginning of the Early-Middle Pleistocene Transition.

1 Introduction

The boundaries of the Early Middle Pleistocene Transition (EMPT), also known as the Mid-Pleistocene Transition, have been set at 1.25 to 0.75 Ma by Clark et al. (2006). This chronological window has been widely discussed and extended to 1.4 to 0.4 Ma by Head and Gibbard (2015). At that time, the Earth's climatic cycles lengthened from 41 to 100 kyr (Pisias and Moore, 1981; Lisiecki and Raymo, 2005), and cycles became asymmetrical, with gradual onsets of glaciation and abrupt glacial terminations (Lisiecki and Raymo, 2005). This difference in structure is linked to a combination of changes in the influence of orbital parameters on the Earth's climate system, and internal changes involving feedback loops such as the long-term cooling and a decrease in CO₂, the removal of regolith and ice dynamics (Chalk et al., 2017; Head and Gibbard, 2015; Berends et al., 2021; Herbert, 2023). Before the EMPT, climate cycles were controlled by obliquity, whereas after the EMPT, they are often presented as being controlled by eccentricity. However, this cyclicity is much more complex, involving the influence of obliquity (41 ka), precession (23 ka) and an eccentricity cycle (100 ka) (Maslin and Brierley, 2015; Herbert, 2023). During the EMPT, the cooling trend that began in the Late Pliocene continued and accelerated in parallel with a decrease in the concentration of atmospheric CO₂ (Martínez-Botí et al., 2015). Glacial periods became more intense and longer (Herbert, 2023). The increase in cycle duration associated with the 100 ka cyclicity is recorded towards the end of the Early Pleistocene, between Marine Isotope Stages (MIS) 37 to 31 (1.24 to 1.06 Ma), comprising the particularly long MIS 36-35 cycle (∼ 75 ka; Barker et al., 2022; Hodell et al., 2023). It is also marked by the MIS 31 “super interglacial” (DeConto et al., 2012; Melles et al., 2012) and the glacial MIS 34, which has recently been characterized as extreme cooling (Margari et al., 2023).

The start of the EMPT coincides with the arrival of the first hominins in Western Europe, as shown by evidence of occupation dating to around 1.4–1.2 Ma in Italy and on the Iberian Peninsula (Peretto, 2006; Carbonell et al., 2008; Moyano et al., 2011; Duval et al., 2024; Cucart-Mora et al., 2026). The Mediterranean is a key region for studying the spread of hominids, because the presence of environmental studies associated with archaeological evidence (Leroy et al., 2011; Messager et al., 2011), and the pathways they took to reach Europe are particularly concentrated in this area (Abbate and Sagri, 2012). In addition, the region is highly sensitive to climate change and is thus appropriate for cross-checking climate results from the high- and mid-latitude systems of the Northern Hemisphere and the African monsoon regimes. The impact of global climate change during the EMPT on the Mediterranean is still poorly understood due to the scarcity of reliable sedimentary archives for this region. Some previous studies have shown that aridification in the Mediterranean was a major parameter throughout this period, leading to changes in vegetation, including forest regression and expansion and the progressive disappearance of relict taxa with affinities for warm and humid environments (Joannin et al., 2011; Pross et al., 2015; Magri et al., 2017; Donders et al., 2021). Foraminifera assemblages reveal that changes in Mediterranean circulation, and particularly variations in Mediterranean outflow water (MOW), resulted in a modulation of Atlantic oceanic circulation, heat transport (Bahr et al., 2018, 2022) and ventilation (Guo et al., 2020). To date, changes in the atmospheric system based on the mineralogy of clays in the Mediterranean at the end of the Early Pleistocene are only continuously documented at ODP Site 964 (Ionian Sea) and MD 964 (Levantine Basin), in the Central and Eastern Mediterranean, respectively. Both records highlight changes in sedimentary sources, varying between North Africa and the southern Apennine Peninsula (Zhao et al., 2016), and changes in precipitation in the Nile source areas linked to monsoon activity (Zhao et al., 2012). Climate reconstruction efforts for this period are still in their infancy. Based on pollen assemblages, only three sequences have been used for climate reconstruction at the end of the Early Pleistocene: Lake Acigöl in Türkiye (Robles et al., 2026), which covers the entire EMPT period, previous data from ODP 976 (Joannin et al., 2011) and the Palominas site (Altolaguirre et al., 2020), which covers the beginning of the transition. Other climate reconstructions based on faunal assemblages from Spanish archaeological and palaeontological sites also provide information on specific time windows, depending on site chronology: Gran Dolina (Blain et al., 2009), Fuente Nueva 3, Barranco Leon (Agustí et al., 2009; Blain et al., 2011, 2016; Saarinen et al., 2021; Sánchez-Bandera et al., 2023), and Almena-Casablanca 3 (Blain et al., 2007; Marquina-Blasco et al., 2025). Comparing different types of climate reconstruction (based on pollen and fauna) can yield a more complete picture of the climate. It is therefore essential to provide new long-term climate data in order to understand global and regional differences in the Mediterranean climate during the EMPT.

Despite the widespread use of marine sequences, the possibilities for studying climate change throughout this period with a multiproxy and integrated approach are limited by the restricted number of long sequences covering the entirety of the EMPT. Thus, our knowledge of the EMPT climate is fragmentary, and questions persist as to the impact of the EMPT on hydroclimatic dynamics, terrigenous sources and vegetation changes. How did the different basins of the Mediterranean respond to climatic changes during the EMPT? How did this transition begin in the Mediterranean? In this paper, we aim to better constrain the land–sea climate dynamics across different Mediterranean basins, with a particular focus on the Alboran Sea, during the onset of the EMPT.

This study provides the first multiproxy and multi-method analysis covering the period ranging from MIS 37 to 31 (1.24 to 1.06 Ma) in the Western Mediterranean using the ODP 976 marine sequence (Alboran Sea). It connects isotope records, pollen data from previous studies (Joannin et al., 2011; Catrain et al., 2025), with new reconstructed Mean Annual Air Temperatures (MAAT) and Mean Annual Precipitation (MAP) based on this pollen data, planktonic foraminifera assemblages, Sea Surface Temperatures (SSTs) based on foraminifera data and clay mineralogy. For this time period, comparisons with the sites of Montalbano Jonico in south Italy (Joannin et al., 2008; Girone et al., 2013), Monte San Giorgio in Sicily (Catrain et al., 2025) and Palominas in southern Spain (Altolaguirre et al., 2020) provide the first land-sea comparison between the central and western Mediterranean showing that differences in vegetation already existed (Catrain et al., 2025). This approach aims to enhance our understanding of climate system variability through: (i) examining links between sea and continent, (ii) a regional focus and comparison of climate response in the Central and Western Mediterranean, and (iii) a methodological discussion of pollen-based climate reconstruction methods.

Figure 1 ODP site 976 climate and location maps. Climate data from ERA5 reanalysis between 1960–2022 with a lat-lon grid of 0.25°. (A) Mean annual precipitation with sea surface circulation. AJ: Atlantic Jet, WAG: Western Alboran Gyre, CCG: Central Cyclonic Gyre and EAG: Eastern Alboran Gyre. (B) Mean annual temperature with sea deep circulation. MOW: Mediterranean Outflow Water, LIW: Levantine Intermediate Water and WMDW: Western Mediterranean Deep Water. (C) Atmospheric circulation. NSAL: Northern Saharan Air Layer.

2 Environmental setting

According to the Köppen classification, the Mediterranean climate is a temperate climate characterized by dry summers (Cs) with precipitation mainly concentrated in winter. Some areas are characterized by more arid steppe (BS) and desert (BW) climates (Köppen, 1936; Daget, 1977; Quézel and Médail, 2003; Peel et al., 2007; Beck et al., 2023). Current annual rainfall around the Alboran Sea varies between 200 and 1100 mm and annual temperatures between 12 and 19 °C, depending on areas and altitudes (Fig. 1). The surface circulation of the Alboran Sea is linked to the entry of the Atlantic jet (AJ) and two main surface gyres: the Western Alboran Gyre (WAG) and the Eastern Alboran Gyre (EAG), both with an anticyclonic circulation. The WAG and EAG are connected to each other by a smaller gyre; the Central Cyclonic Gyre (CCG) (Fig. 1A) (Sánchez-Garrido and Nadal, 2022). The surface waters from the Atlantic Ocean entering through the Strait of Gibraltar move towards the Almeria Oran Front (AOF) where they mix with saltier waters. From the east, the Levantine Intermediate Water (LIW) and the colder and more saline Western Mediterranean Deep Water (WMDW) move west before leaving through the strait. The mixing of these two water masses upon arrival in the Atlantic Ocean is called Mediterranean Outflowing Water (MOW) (Fig. 1B) (Fabres et al., 2002; Macías et al., 2008; Moreno et al., 2005; Pérez-Folgado et al., 2003; Tintore et al., 1988).

In the Western Mediterranean region, atmospheric circulation is influenced by the westerlies and variations in the positioning of the intertropical convergence zone and therefore the monsoon system. The Western Mediterranean is subject to winds from Africa, such as Northern Saharan pluri-directional winds, which originate in the Sahara. On the Atlantic side, the Saharan Air Layer (SAL) originating from the monsoon system is deflected by the Azores high pressure system, thus creating the northern branch of the SAL (NSAL), which goes up to the Western Mediterranean (Bout-Roumazeilles et al., 2007; Leblanc, 2022). The Alboran basin is subject to westerly and south-westerly winds, the Poniente and the Vendaval, and warm easterly winds such as the Levanter.

Around the Alboran Sea, Mediterranean vegetation is organized altitudinally according to temperature and precipitation distribution. The Thermo-Mediterranean belt (0 to 300–700 m) is mainly composed of Mediterranean species, such as olive trees (Olea), pistachio trees (Pistacia) and carob trees (Ceratonia). The Meso-Mediterranean (300–700 to 1000–1200 m) and Supra-Mediterranean (1000–1200 and 1500–1700 m) belts consist of evergreen mixed forest and deciduous oak (Quercus) and pine (Pinus) forests respectively. Above these altitudes, forests are composed of conifers such as pine (Pinus), fir (Abies) and cedar (Cedrus) in Morocco (Benabid, 1985; Ozenda, 1994; Quézel and Médail, 2003; Médail, 2022).

3 Site location and chronology

ODP Leg 161 site 976 (36°12.30^′ N, 4°18.8^′ W, 1108 m b.s.l. (below sea level); Fig. 1) is located in the northwest of the Alboran Sea. The sediment core records a continuous lithostratigraphy from the Upper Miocene to the Holocene (Comas et al., 1996). Our study focuses on the depth interval of the master core from 281.23 to 256.53 mcd (metre composite depth), corresponding to the MIS 37-31 interval, with an age model based on nannofossil events, sapropel and the correlation between ODP site 976 (Catrain et al., 2025) and global stack oxygen isotope curves (Lisiecki and Raymo, 2005). The sedimentary sequence consists of homogeneous clays rich in nannofossils and silty nannofossils. The lithogenic particles in the Alboran Sea come from aeolian and fluvial sources from Spain and Morocco.

4 Materials and Methods

This paper uses the ODP Site 976 palynological data obtained by Catrain et al. (2025) which followed a standard protocol of chemical treatment. For every sample a minimum of 150 grains excluding Pinus (overrepresented) and a minimum of 15 taxa were recorded to a have a reliable representation of the vegetation.

Additional proxies and statistical analyses were used to cross environmental information and to better constraint the land-sea climate dynamics.

4.1 Planktonic foraminifera analysis and isotope analysis

Planktonic foraminifera analysis was carried out on one hundred and sixteen samples with a spacing of 10 to 80 cm, with an average of every 23 cm. The sediments were sieved to 150 µm. Foraminifera were counted using a binocular magnifying glass with a minimum of 450–500 tests per sample to be representative. Identification is based on publications (Orbigny, 1839) and a book (Schiebel and Hemleben, 2017). No evidence of diagenesis was identified.

Between eight and 12 individuals of the foraminifera Globigerina bulloides were selected from 56 of these same samples to conduct δ¹⁸O and δ¹³C analyses. The samples were cleaned with methanol before being inserted into the dual-entry Isoprime isotope ratio mass spectrometer. Analyses were performed at the Laboratory for Sciences of Climate and Environment (LSCE, France) on the analytical platform PANOPLY.

4.2 Clay mineral analysis

Clay mineral associations were studied using X-ray diffraction following the protocol of Bout-Roumazeilles et al. (1999). In this study, “clay minerals” refer to the main phyllosilicate minerals in the clay-size fraction (generally less than 2 µm particles). All samples were first decalcified with 0.2 M HCl, before clay defloculation by successive washing with distilled water. The clay-size fraction was separated by settling according to Stokes's law, concentrated by centrifugation, and oriented by wet smearing on glass slides. X-ray diagrams were obtained using a Bruker D4 Endeavor diffractometer with CuKα radiation and a Ni filter, coupled with a Lynxeye detector. A tube voltage of 30 kV and a tube current of 35 mA were utilised. Three X-ray diagrams were performed: air-dried sample (normal run), ethylene-glycol vapour saturation for 12 h (glycol run) and heating at 490 °C during 2 h (heating run). The goniometer scanned from 2.49 to 32.49° 2θ for normal and glycol run and from 2.5 to 14.5° 2θ heating run. Each clay mineral is characterized by its basal layer plus interlayer interval (d), as revealed by XRD analysis (Brown and Brindley, 1980). Smectite (S) is characterized by a peak at 14 Å on the normal run, which expands to 17 Å after saturation by ethylene-glycol and retracts to 10 Å after heating. Illite (I) presents a basal peak at 10 Å on the three runs. Palygorskite (P) is characterized by the presence of a peak at 10.34Å on the normal run. Chlorite (C) is characterized by peaks at 14 Å (001), 7 Å (002), 4.75 Å (003) and 3.53 Å (004) on the three runs. Kaolinite (K) is characterized by peaks at 7 Å (001) and 3.58 Å (002) on the normal and glycol runs. Both peaks disappear or are strongly reduced after heating. To distinguish kaolinite from chlorite, the portion of the spectrum containing the basal peaks of kaolinite (002) and chlorite (004) around 3.55 Å is step-scanned in a high-resolution mode following standard procedures described in detail by Petschick et al. (1996). Semi-quantitative estimation of clay mineral abundance is based on peak areas and totalled to 100 % (S + I + K + C + P =100 %). Peak area measurements were taken in the glycol runs using the Macintosh MacDiff^® 4.2.5 software (Petschick, 2000). The error on measurement reproducibility is estimated to be ±5 % for each clay mineral.

4.3 Sea surface temperature reconstruction from planktonic

foraminifera

Sea surface temperatures (SST) were reconstructed based on planktonic foraminifera assemblages, using PaleoAnalogs software (version 3.0). This software uses the Modern Analogue Technique (MAT), which compares a database of modern assemblages with fossil assemblages to find the closest modern samples (Theron et al., 2003). The degree of dissimilarity between samples is calculated using the Squared-Chord-Distance. A minimum of ten temperature analogues were selected and then averaged to estimate those of the fossil samples. Reconstructions were based on a North Atlantic Mediterranean database (Siccha and Kucera, 2017).

4.4 Pollen-inferred climate reconstructions

For climate reconstructions, pollen counts were used to reconstruct four climatic parameters: Mean Annual Air Temperature (MAAT), Mean Temperature of the Coldest month (MTCO), Mean Temperature of the Warmest month (MTWA), and Mean Annual Precipitation (PANN). A multi-method approach was adopted for greater reliability (Peyron et al., 2011, 2013; Salonen et al., 2019; Chevalier et al., 2020), and taxa that have now disappeared from the Mediterranean (relict taxa) were included in certain methods (Fauquette et al., 1998). Five methods with different mathematical bases were selected: the modern analogue technique (MAT; Guiot, 1990), weighted averaging partial least-squares regression (WAPLS; Ter Braak et al., 1993), random forest (RF; Breiman, 2001; Prasad et al., 2006), boosted regression trees (BRT; De'ath, 2007; Elith et al., 2008) and the climatic amplitude method (CAM; Fauquette et al., 1998). To discuss the climate of the Western Mediterranean Basin with a multi-site approach, additional pollen-based climate reconstructions were carried out using data from Monte San Giorgio and Palominas, following the same multi-method approach (CAM, MAT, WAPLS, RF and BRT).

The CAM is based on the overlap of the climatic requirements of taxa and takes relict taxa into account. The climatic requirements of ∼ 135 taxa were defined on the basis of 8000 modern spectra and the literature. The past climate is estimated by transposing the climatic requirements of the maximum number of modern taxa to the fossil data. The most probable climate corresponds to the suitable climatic interval for the highest number of taxa and to a “Most Likely Value” corresponding to a weighted mean (Fauquette et al., 1998). Some taxa represented by a single grain in the record (Lygeum, Croton, Rhus, Myrthus, Calligonum, Taxodium, and Sapotaceae) and Cichorioideae were not taken into account in the climate reconstruction as they generate noise. For the CAM, Poaceae were not taken into account for the reconstruction of the Palominas site for the same reason as for Cichorioideae in ODP site 976.

RF and BRT methods are based on machine learning and have recently been used for palaeoclimatology in Northern Europe (Salonen et al., 2019) and the Mediterranean Basin (d'Oliveira et al., 2023, 2025; Robles et al., 2023, 2025; Sassoon et al., 2025; Charton et al., 2025). For RF, each tree is estimated from a random set of different subsets defined separately from the modern data assembly. This is called bootstrapping and each sample has the same probability of being selected (Breiman, 2001; Prasad et al., 2006). For BRT, the random creation of sets from modern assemblages means that samples that were not sufficiently taken into account in the previous tree have a higher probability of being selected in the next training sessions. This is known as boosting (De'ath, 2007; Elith et al., 2008). WAPLS is a non-linear regression technique that creates relationships between climate parameters and pollen taxa within the modern database, then applies them to fossil assemblages (Ter Braak et al., 1993; Chevalier et al., 2020). MAT method works in a similar way to the one use for reconstruct SST based on foraminifera. The MAT and WAPLS methods were applied with the R rioja package (Juggins, 2024), RF with randomForest package (Breiman, 2001) and BRT with dismo package (Hijmans et al., 2024).

The Eurasian and Mediterranean Current Surface Sample Database (EAPDB; n=3373 sites), compiled by Peyron et al. (1998, 2011, 2013), and completed by Dugerdil et al. (2021) and Robles et al. (2023), was used for the first four methods. In this dataset, each sample is associated with a biome (Tarasov et al., 1998; Elenga et al., 2000). Some biomes were selected and others were not taken into account in order to better constrain the database. Here the biomes selected: cold deciduous forest (CLDE), cold desert (CODE), cool mixed forest (COMX), cold steppe (COST), hot desert (HODE), pioneer (PION), temperate deciduous forest (TEDE), warm mixed forest (WAMX), warm steppe (WAST) and xerophilous (XERO). Anthropic (ANTH), cold mixed forest (CLMX), cold conifer forest (COCO), taiga (TAIG) and tundra (TUND) biomes were not taking into account because they are too far from the study area. One hundred and seven sites from the Eurasian modern database (EMPD) (Davis et al., 2020) were added to the EAPDB for a more reliable representation of certain taxa (Ericaceae). Particular attention was paid to these new samples to avoid selecting anthropized sites. Pinus, kept in the modern database, was not included in the fossil sequence due to its overrepresentation, except for the Palominas site which corresponds to a continental basin.

It is difficult to rank climate reconstruction methods, as each has its own advantages and disadvantages. The order in which the methods appear in this section reflects our subjective ranking. CAM offers a key advantage for this time period because, unlike other methods, it accounts for relict taxa. Furthermore, it considers taxa individually, allowing them to be distinguished within pollen assemblages. Machine learning methods such as BRT and RF do not rely on ecological theories but are based on statistical methods, which helps limit ecological biases in the database. MAT appears to be the least suitable method for this study period, as there is no known true analogue of past vegetation in the Mediterranean region, hence the use of a database covering Eurasia to identify the closest possible assemblages. However, the dominance of Quercus in the arboreal component of the pollen signal makes it difficult to constrain the reconstructions effectively, regardless of the method used, since trees exert a stronger constraining influence on climate reconstructions than grasses due to their higher taxonomic resolution.

4.5 PERMANOVA

The use of PERMANOVA in this study aims to compare different methods of climate reconstruction with one another, taking into account different time periods and climate parameters. PERMANOVA is a non-parametric test based on the F-test. The p-value of PERMANOVA is estimated by permutation, which makes it particularly well adapted to data that do not follow multivariate normality or homogeneity of variances, as is the case with climate reconstruction data. These two criteria are important for the application of traditional analyses of variance methods (ANOVA, MANOVA), as confirmed using a Mardia test. Whilst ANOVA compares group means to verify whether or not they come from the same population, PERMANOVA projects data points into multivariate space, and therefore compares group centroids instead of means. The null hypothesis (H0) of the PERMANOVA test is that the reconstructed values belong to the same centroids and that their dispersion is equivalent across all groups. In other words, reconstructed values are considered as permutable or “equivalent” under the null hypothesis. When the p-value of PERMANOVA is less than or equal to 0.05, the test is significant and rejects H0. The centroids of reconstructed values and their spread is different for each group, which are considered not permutable in multivariate space and therefore, “different”. This test was applied to reconstructed temperature and precipitation data, which were centred and reduced to avoid the influence of scale and give equal weight to each variable and to compare trends.

5 Results

As a reminder, the pollen data presented by Catrain et al. (2025) are composed of 31 families and 72 genera distributed over the different belts of Mediterranean vegetation. The temperate taxa curves, mainly represented by the Quercus deciduous type, Isoetes and elements of steppe environments (Artemisia, Amaranthaceae, Ephedra) highlight alternating warm/humid periods and cold/dry periods between MIS 37 and 31. The high percentages of high-altitude forest elements, such as Cedrus and Abies, and of Ericaceae, point to the onset of colder but still humid climatic conditions at the end of the interglacial. The constant values of Asteroideae, Cichorioideae and Poaceae show a continuous presence of an open vegetation environment (Fig. S1 in the Supplement).

Figure 2 Planktonic foraminifera assemblages from ODP site 976 and SST reconstructed with the Modern Analog Technique (MAT) associated with MIS chronology. The red curves correspond to foraminifera with warm affinities and the blue curves to cold affinities. A map of the analog distribution (Fig. S2) and a table of analog's distances (Table S1) are available in the Supplement.

5.1 Late Early Pleistocene planktonic foraminifera assemblages and SST reconstruction

Altogether, ninety-seven samples were analysed. However, fourteen samples could not be used due to the low number of foraminifera. Thirty-nine species were identified, dominated by Neogloboquadrina incompta (on average 30 %), Globigerina bulloides, Globorotalia inflata, Neogloboquadrina pachyderma, and Globigerinoides ruber white (Fig. 2). Globigernella calida (around 1 %), Globigerinita glutinata (4 %), and Globigerina falconensis (3 %) are consistently present. Taxa with affinity for warm waters, such as Orbulina, Globorotalia crassaformis, and G.calida, are poorly represented in the assemblage. Warm taxa, Globigerinoides ruber white type, and taxa with an affinity for temperate waters, such as Globigerina bulloides and Globorotalia inflata, are found in interglacial periods. Two peaks of G Globigerinoides ruber white are recorded at 1229 (39 %) and 1084 ka (49 %). Turborotalita quinqueloba develops significantly (up to 30 %) between 1215 and 1200 ka. The high Neogloboquadrina pachyderma percentages, with peaks between 40 % and 70 %, are synchronized with glacial periods.

Mean annual sea surface temperature values reconstructed from planktonic foraminifera assemblages vary between 5 and 19.6 °C. The values show marked variation with a trend of cooling during glacial periods and warming during interglacial periods (Fig. 2).

Figure 3 Clay mineralogy from ODP site 976 associated with MIS chronology. Solid lines represent the mean and dotted lines the standard deviation.

5.2 Late Early Pleistocene clay analysis in the Alboran Sea

Over the studied interval, the clay assemblage is mainly composed of illite (35 %) and smectite (30 %) associated to 15 % chlorite and 10 % kaolinite, palygorskite being less abundant (6 %). This composition is very similar to the composition of recent sediments deposited in the Alboran Sea (MIS1-MIS3), consisting mainly of illite and smectite, which account on average for 35 % and 30 % of the clay assemblage respectively (Bout-Roumazeilles et al., 2007; Hassouni and El Hmaidi, 2007), with kaolinite and chlorite at around 16 %, whilst palygorskite accounts for less than 5 % (Bout-Roumazeilles et al., 2007; Hassouni and El Hmaidi, 2007). The clay mineral signals can be divided into four main phases of variation between stages 37 to 31 (Fig. 3):

• Phase I, corresponding to MIS 37 is slightly enriched in illite (40 %) and depleted in smectite (20 %), with a I $/$ K ratio around 2.4. The end of stage 37 is marked by an increase in kaolinite (∼ 20 %–21 %) and chlorite (∼ 19 %) and a decrease in smectite (∼ 12 %) and the I $/$ K ratio (∼ 2).

• Phase II is the longest phase, covering MIS 36-35, and records fairly constant values for each mineral, very similar to those for phase I.

• Phase III, which begins with the transition from MIS 35 to 34, shows the greatest change in clay composition. A general decreasing trend is observed for illite (∼ 31 %), kaolinite (∼ 10 %), chlorite (∼ 11 %) and palygorskite (∼ 6 %), although the latter peaks at around 1114, 1103, 1001 and 1088 ka. On the other hand, during this phase, smectite increases significantly, reaching average values of around 42 %. The I $/$ K ratio fluctuates more regularly, with values rising slightly to 3.69.

• Phase IV corresponding to MIS 31 presents a return to the values of phases II and I.

Figure 4 ODP site 976 pollen climate reconstruction. In light purple, RF (Random Forest), in purple BRT (Boosted Regression Trees), in yellow WAPLS (Weighted Averaging Partial Least Squares regression), in light brown MAT (Modern Analogue Technique) and in green CAM (Climatic Amplitude Method). Four climatic parameters: MAAT (Mean Annual air Temperature), MTWA (Mean Temperature of the WArmest month), MTCO (Mean Temperature of the COldest month) and PANN (mean annual precipitation). Grey bands correspond to interglacial periods and the thick coloured lines correspond to loess-smoothed curves and shaded areas to the 95 % confidence interval.

Table 1 R² (Pearson correlation coefficient) and RMSE (Root Mean Square Error) values for climatic parameters. The lower the RMSE and the higher the R², the closer the method is to the target value from the modern database. R² and RMSE are based on bootstrapping for MAT and WAPLS and on leave-one-out validation for BRT and RF. The calibration range associated is given in Table S2 and the WAPLS has 3 components. Figure S3 presenting the observation vs predicted plots. For CAM, R² and RMSE are obtained by applying the method to an independent validation dataset and comparing the predictions with the observed data.

5.3 Climate reconstruction during the Late Early Pleistocene

Reconstructions from pollen data at ODP site 976 show similar trends between all the different methods and for each parameter (Fig. 4). BRT and MAT have the best R² and RMSE, while WAPLS yield the poorest for all climatic parameters (Table 1). MAT and WAPLS show greater variation amplitude than BRT, RF and CAM, where differences between warm, wet phases and cold, dry phases are less pronounced. The RF shows the same variations as the CAM for the MAAT and PANN; for the MTWA, the CAM exhibits greater amplitudes, whereas for the MTCO the RF has greater amplitudes. Parameter variations follow those of interglacial and glacial periods. Interglacial phases are characterized by higher temperatures and precipitation than for glacial periods. RF appears colder than other methods, with a general mean value of around 10 °C in MAAT and 0 °C for MTCO. WAPLS and MAT present more humid conditions.

All methods except BRT show that MIS 37 (1245–1215 ka) was the warmest period of the study, based on MAAT, MTWA and MTCO. Precipitation levels were high, but were nevertheless similar to those reconstructed for MIS 31 (1081–1062 ka), which had lower temperatures. MIS 35 (1190–1141 ka), meanwhile, shows temperature values close to those of MIS 31 but with slightly lower precipitation. For this stage, the BRT is the only method that shows a decrease in annual temperatures and in the temperature of the warmest month. The RF, MAT and WAPLS methods indicate that MIS 33 (1114–1104 ka) was the coolest and driest interglacial of the period.

According to all methods, MIS 36 (1215–1190 ka) was characterised by cold and dry conditions, although these were less severe than during the subsequent glacial phases; the only exception is the BRT method, which shows no decrease in temperature or precipitation relative to the flanking interglacial periods. According to the MAT and WAPLS methods, the MIS 34 glacial phase (1141–1114 ka) is the coldest and driest of the cold phases of the period for the MAAT and MTCO. The other methods show that the end of MIS 32 (1104–1081 ka) was the most extreme period. However, all methods unanimously indicate that MIS 34 is the stage with the driest conditions.

Figure 5 PERMANOVA p-values between climate reconstruction methods for each isotopic stage and for Mean Annual Air Temperature (MAAT) and Mean of Annual Precipitation (PANN). A p-value less than or equal to 0.05 is significant and indicates that the two methods are different from each other.

5.4 Statistical analysis of climate reconstruction

The test comparing the methods for each MIS (Fig. 5) shows that, for most of them, the values reconstructed by the methods belong to the same centroid. Indeed, for PANN, H0 is not rejected for any methods from MIS 37 to MIS 33 (Fig. 5). For MAAT, H0 cannot be rejected for any of the methods (Fig. 5). Therefore, the various tests carried out with PERMANOVA on the centred and reduced climate reconstruction data suggest that the methods can be considered as consistent between each other. H0 is only rejected for certain MIS or for PANNs. For MIS 32, only CAM has a p-value <0.05 and is significantly different to MAT (p-value =0.004), WAPLS (p-value =0.004) and BRT (p-value =0.005). RF differs from MAT (p-value =0.039) and BRT (p-value =0.049) during MIS 31. In those cases, reconstructed values are considered to be significantly different between each other (Fig. 5).

PERMANOVA teststhat compare the glacial and interglacial periods with each other for each method, however, all reject H0 (p-value <0.05; Fig. 6). All interglacial and glacial phases can therefore be considered to be different from the others, regardless of the methods.

Figure 6 PERMANOVA p-values between the same types of periods (glacial-glacial and interglacial-interglacial) for each climate reconstruction method on centred-reduced data. A p-value less than or equal to 0.05 is significant and indicates that the two periods are different from each other.

6 Discussion

6.1 Land-sea responses from the Western to Central Mediterranean

Planktonic foraminifera data from the same marine sequence as the pollen data highlight simultaneous changes in oceanic characteristics and in the floristic composition of terrestrial ecosystems. The Mediterranean Basin, and in particular the Western Mediterranean, only comprises rare sites combining these two proxies. This foraminifera/pollen data approach is recorded in the Montalbano Jonico site over a large EMPT time window (Joannin et al., 2008; Girone et al., 2013).

Figure 7 Land-sea variability comparison at ODP site 976 (left), compared to the site of Montalbano Jonico (right) in southern Italy. Isotope data for the ODP site 976 are from Catrain et al. (2025). For the Montalbano Jonico site, pollen and oxygen isotope data come from Joannin et al. (2008), and planktonic foraminifera data from Girone et al. (2013).

Between MIS 37 and MIS 31 in ODP site 976 (Western Mediterranean basin), all the proxies highlight the typical interglacial and glacial variations recorded during this period, except clay mineralogy, which shows longer cycles. The SST estimates should be interpreted with caution, as the high abundance of Neogloboquadrina pachyderma is likely to introduce a cold bias, as already observed in a previous study (Mega et al., 2025). However, the SST variations are of great interest and may be compared to the other proxies. The comparison between pollen data and planktonic foraminifera assemblages highlights a synchronization of steppe elements with the polar affinity foraminifera Neogloboquadrina pachyderma. A similar pattern is observed during MIS 20-19 (Toti et al., 2020). This correlation emphasizes that the steppes from the Iberian Peninsula and/or North Africa are not only linked to aridity but also to low temperatures (Fig. 7a). The constant dominance of the subpolar affinity Neogloboquadrina incompta (∼ 30 %) throughout the period also attests to relatively cool surface waters (average SST of 12 °C). Air and ocean temperatures (MAAT and SST) reconstructed from pollen and planktonic foraminifera also show these glacial-interglacial variations (Fig. 7a). The SST signal is highly influenced by the quantity of Globigerinoides ruber of subtropical affinity during interglacial phases, and Neogloboquadrina pachyderma of polar affinity, during glacial phases. For both air and ocean temperatures, these are lower than modern values. The current SST is 17.7 °C between 1981–2010 (36° N, 03° W; Guijarro et al., 2015; Sánchez-Laulhé et al., 2021) and the current mean air annual temperature in a perimeter of 400 km around ODP 976 is 16.8 °C between 1960–2022 (calculated with ERA5 reanalysis). Furthermore, multiproxy analyses from ODP site 976 show that the warmest phase is MIS 37 (1245–1215 ka). During this phase, temperate forests reached their maximum development, averaging 35 %, while steppes were reduced to 10 %. This stage also has the highest percentages of Globigerinoides ruber (∼ 20 %), a warm-water foraminifera, and SSTs around 16.5 °C, as well as low δ¹⁸O values (∼ 0.56 ‰; Catrain et al., 2025) (Fig. 7).

This observation counters current knowledge, which presents MIS 31 (1081–1062 ka) as the warmest interglacial during this period. It is generally referred to as the “Super Interglacial”, associated with maxima in obliquity, eccentricity and insolation (DeConto et al., 2012; Melles et al., 2012). However, this result is consistent with interpretations made at sites off the Iberian Peninsula. Indeed cores U1385 and U1387 show that MIS 31 is not associated with particularly high temperature conditions compared to those recorded at high latitudes (Voelker et al., 2015; Oliveira et al., 2017). Finally, in ODP site 976, MIS 31 is less marked than MIS 37 in terms of isotope values, pollen and foraminifera assemblages, and is closer to MIS 33 (1114–1104 ka) and MIS 35 (1190–1141 ka).

Concerning glacial phases, MIS 34 (1141–1114 ka) seems to be the most pronounced, with significant steppe expansion (∼ 28 %) and forest regression (10 %). Neogloboquadrina pachyderma reached 40 % and reconstructions estimate an annual air temperature of ∼ 9 °C, which is corroborated by high δ¹⁸O values of up to 1.62 ‰. However, these values remain very close to those observed during MIS 36 (1215–1190 ka). Reconstructed SSTs show colder values for MIS 36 (9.5 °C) than MIS 34 (9.8 °C), even though they remain within the same range of values. This slight difference can be explained by the combined presence of Neogloboquadrina pachyderma and the subpolar taxa Turborotalita quinqueloba in MIS 36, whereas the latter is almost absent in MIS 34.

6.2 Provenance of clay minerals in the Western Mediterranean

Variations in sediment mineralogy – with smectite ranging from 20 % to 60 %, illite from 20 % to 45 %, and kaolinite present in lower proportions, up to 20 % – suggest changes in the provenance and/or transport patterns of clay particles during the studied interval.

However, the origin and interpretation of smectite in the Alboran basin is still debated (Holeman, 1968; Milliman and Syvitski, 1992; Stanley et al., 1992; Krom et al., 1999; Weldeab et al., 2002; Garzanti et al., 2006; Ehrmann et al., 2007; Hassouni and El Hmaidi, 2007; Hamann et al., 2009; Zhao et al., 2012; Bout-Roumazeilles et al., 2013). Smectite is ubiquitous in various types of environments, as it is easily transported over long distances due to its buoyancy. There are several hypotheses regarding the potential sources of the smectite collected in the Alboran Basin, as well as for illite and kaolinite : (1) fluvial origin (the Nile, major European rivers, coastal rivers of the Iberian Peninsula, Moroccan wadis, etc.), (2) eolian origin, and (3) advection/redistribution by oceanic water masses (turbidites, Suspended Particulate Matter (SPM)) from the central and eastern Mediterranean, or from the Gulf of Cádiz and the North Atlantic via the Gibraltar Strait.

The Nile, according to its average discharge (from 120 to 230.10⁶ t yr⁻¹) is one of the primary fluvial sources of smectite in the eastern Mediterranean basin (Venkatarathnam and Ryan, 1971; Foucault and Mélières, 2000; Bout-Roumazeilles et al., 2007; Hamann et al., 2009; Zhao et al., 2012) but previous studies have suggested that clays originating from the Nile dispersed only negligibly westward beyond the Mediterranean Ridge, (Krom et al., 1999; Weldeab et al., 2002).

The river systems in southern Spain are unlikely to influence clay sedimentation in the Alboran Sea. Indeed, the Guadalquivir flows into the Gulf of Cádiz, on the other side of Gibraltar, and the Segura, Júcar and Ebro river systems have no outlet into the Alboran Sea, as they flow further north.

The Po River, with an annual discharge of 17.10⁶ t yr⁻¹, as well as most SE European rivers, is mainly supplying illite toward central Mediterranean flows via the Adriatic Sea and the Aegean Sea (Venkatarathnam and Ryan, 1971; Dominik and Stoffers, 1979; Chamley, 1989; Alonso and Maldonado, 1990; Tomadin, 2000). The Rhone River may supply illite and chlorite to the Gulf of Lion, but its outlet is distant from the studied area, reducing its direct influence on the Alboran Sea sedimentation.

Riverine contribution to the central and western Mediterranean from areas bordering the southern shores is generally enriched in illite and is severely reduced due to the narrow drainage basin and sparse rainfall (Martin and Milliman, 1997). Nevertheless, the Medjerda and Miliane rivers, discharging into the Gulf of Tunis, may influence the sedimentation in the study area. Smectite is a dominant component of suspended particulate material in the Alboran Sea (Pierce and Stanley, 1975; Grousset et al., 1988; Baringer and Price, 1999), especially within the Western Alboran Gyre (Bulian, 2022), resulting from the interaction between the smectite-rich entering surface Atlantic water and smectite-rich Levantine Intermediate Water and Western Mediterranean Deep Water (Ambar and Howe, 1979; Bryden and Stommel, 1984; Millot, 2014).

Given the physiographic context of the Alboran Sea, part of the terrigenous fraction may be transported from the South Iberian Plateau to ODP Site 976 via submarine canyon/delta systems such as the Baños, Calahonda or Fuengirola systems, which allow sediments to be transported to depths of over 1000 m (Ercilla et al., 1992; Alonso et al., 1999; Vázquez et al., 2015). Nevertheless, the sedimentological description of ODP 976 (Comas et al., 1996) indicates that sedimentation over the studied interval is predominantly hemipelagic with no major turbiditic intervals. This contrasts with the older Pliocene section, characterised by numerous gravitational episodes, reflecting reduced sediment transport through these canyons during the interval studied. These observations suggest that gravitational transport processes did not play a major role in sedimentation. Sea-level changes may also affect the clay composition through modifications of the source-to-sink transfers. During glacial period, low sea-level would promote the transfer from South Spain riverine systems enriched in illite and kaolinite and depleted in smectite, which is however not supported by the clay mineral assemblages at site ODP 976. Similarly, during deglaciations, sea level rise would favour the rapid remobilization of sediments from the shelf. By contrast, deglaciations are not associated with major modifications of the clay mineral fraction, given low confidence in this hypothesis.

The advection of smectite-rich sediments from the Gulf of Cádiz via surface water masses entering the Mediterranean is a potential candidate for smectite supply. A recent study evidenced that sediment drifts from the Gulf of Cádiz likely formed under the influence of the MOW are enriched in smectite (Moal-Darrigade et al., 2022). As a consequence, we should not rule out the possibility that sediments passing through Gibraltar are returned to the Mediterranean by ocean circulation, thereby allowing accumulation in the Alboran Basin via the WAG gyre.

Finally, the generally low direct fluvial input and potential eolian contribution, as well as the overall presence of palygorskite suggest that a significant proportion of the clay sediments may be of eolian origin. Previous studies evidenced the dominance of illite and kaolinite, associated with palygorskite in eolian dust from northern Africa, and highlighted the scarcity of smectite in Saharan dust (Guerzoni et al., 1999; Bout-Roumazeilles et al., 2007, 2013).

6.3 Atmospheric and oceanographic implications

To highlight the potential eolian contribution to marine clay sedimentation, the comparison between palygorskite content and I $/$ K ratio may be used to refine the identification of the various eolian sources. Physical erosion and chemical weathering are responsible for the diversity of clay assemblages, but previous studies (Thiry, 2000), on the paleoclimatic interpretation of clay minerals in marine deposits, showed that direct interpretation of variations in the composition of the clay fraction may be tricky. In particular, processes like chemical weathering and hydrolysis, which are responsible for the transformation of clays in soils, are slow-acting mechanisms, which raises issues of synchronisation and time lag between the moment when hydrolysis conditions change, the moment when the nature of the minerals formed changes, and the moment when they are finally transported and then deposited in the oceanic environment. Based on these assumptions, any direct link between modifications of weathering conditions on the continent and contemporaneous variations in clay composition within the Alboran sea is difficult to assess and is not supported by the whole dataset.

Based on these findings, the clay mineralogy analysis likely records changes in atmospheric and oceanographic conditions in the Western Mediterranean with reduced impact of gravity transfer processes. Considering the almost permanent presence of palygorskite, and based upon previous studies in the central and western Mediterranean including the Alboran sea (Bout-Roumazeilles et al., 2007, 2013), the illite/kaolinite (I $/$ K) ratio is used in order to tentatively retrace variations in dust sources from North Africa. Currently, illite decreases from 60 % of clay content in Northern Algeria to 30 % in the Sahelian area (Paquet et al., 1984), whereas kaolinite, scarce in the northern Sahara, increases progressively southward and eastward and becomes dominant in the Sahel (Caquineau et al., 1998; Guerzoni et al., 1999; Bout-Roumazeilles et al., 2007; Stumpf et al., 2011). A high I $/$ K ratio therefore suggests inputs preferentially from the Sahara, while lower values indicate a provenance associated with the Sahel. In ODP 976, the I $/$ K ratio strongly increases after the transition between MIS 35 and 34 (around 1140 ka), reflecting a dominant Saharan origin for aeolian particles (Fig. 7). This suggests the development of arid conditions with an enhanced supply of Saharan vs Sahelian origin at that time (Bout-Roumazeilles et al., 2007; Skonieczny et al., 2011).

While variations in the I $/$ K ratio are commonly used to retrace aeolian supplies, the interpretation of individual clay minerals like smectite is more complex. Given the various scenarios outlined above, the variations observed in smectite from ODP 976 can be interpreted as reflecting the oceanic transport of terrestrial material into the Alboran Sea, whereas the I $/$ K ratio may partly retrace aeolian supply. As a consequence, the shift in the smectite signal is observed at ∼ 1140 ka, concomitant to the observed shift in the I $/$ K ratio, suggesting a change in terrigenous advection and thus in the oceanic configuration in the Alboran Sea.

In the Central Mediterranean, the vegetation signal from Montalbano Jonico shows extensive steppe developments during glacial periods and different interglacial periods compared to the western basin, more pronounced in MIS 31 and 35 than MIS 37 (Fig. 8b). This shows the influence of the Atlantic mitigated forest development during the MIS 31 super interglacial period, compared to the Eastern Mediterranean. In marine environments, the signal also differs from ODP site 976 with the dominance of Turborotalita quinqueloba and higher values for Globigerinoides ruber and lesser values for Neogloboquadrina pachyderma and N. incompta. T. quinqueloba is a cold-water taxon which represents very productive waters and reflects the distribution of water masses and the stratification of the water column (Capotondi et al., 2004; Girone et al., 2013; Schiebel and Hemleben, 2017; Bazzicalupo et al., 2018). The central Mediterranean waters therefore appear warmer than in the western basin, but still relatively cool as attested by the presence of T. quinqueloba. This is consistent with a lesser influence of the input of cold Atlantic waters through the Gibraltar Strait compared to the Alboran Sea. Throughout this period, Turborotalita quinqueloba variations appear to be partly linked to interglacial-glacial variations. At the end of MIS 35, it experienced a drastic decline from 80 % to 20 %. This major decline coincides with changes in the percentage of smectite in the Alboran basin (Fig. 7). However, this event is not found in the clay mineralogy from ODP Site 964 in the Ionian Sea (Zhao et al., 2016), nor further east at Site MD90-964 in the Levantine Basin (Zhao et al., 2011). This suggests that the change in marine circulation in the Alboran Sea is likely associated with Atlantic circulation and attenuates eastwards in the Mediterranean.

6.4 Late Early Pleistocene climate in the Southwest Mediterranean

Among the small number of continuous palynological sequences from the Early-Middle Pleistocene period in the Mediterranean, only two provide climate reconstruction data; Palominas site in southern Spain (Altolaguirre et al., 2020) with the Coexistence approach (CA) and Acigöl site in Türkiye (Robles et al., 2026), with a multi-method approach (CAM, MAT, WAPLS, RF and BRT).

Figure 8 Mean Annual Air Temperature (MAAT) and Annual Precipitation (PANN) reconstructions of three sites in the Western Mediterranean basin between 36 and 38° N: ODP site 976, Palominas and Monte San Giorgio. Five methods were used: Weighted Averaging Partials Least Square (WAPLS), Modern Analog Technique (MAT), Random Forest (RF), Boosted Regression Tree (BRT) and Climatic Amplitude Method (CAM). For Palominas, the Coexistence Approach (CA) method was performed by Altolaguirre et al (2020). m b.s.l.: metres below sea level; m a.s.l.: metres above sea level. The dotted lines correspond to modern values: for the ODP site 976 the mean value for a 400 km radius is taken from the ERA5 reanalysis dataset (1960–2022), for Palominas and Monte Sangiorgio the values are taken from http://www.climate-date.org (last access: last access 3 July 2025) respectively for the cities of Baza and Caltagirone.

The three Western Mediterranean sites present consistent temperature variations (MAAT), which mirror the alternation of glacial and interglacial phases (Fig. 8). The MAT, WAPLS, BRT, and CAM methods show average annual temperatures of around 11 to 12 °C for both marine sites (ODP site 976 and Monte San Giorgio) over the entire period, while the RF drops to an average of 10 °C (Fig. 8a, c). In the continental sequence of Palominas, average temperature varies depending on the method used: 15 °C for CA, 12 °C for CAM, 10 °C for MAT and BRT, 9.5 °C for RF and 8 °C for WAPLS (Fig. 8b). The lower values appear to be consistent with the altitude of the site, around 780 m above sea level. However, these values seem very cold for annual temperatures in this Southern Mediterranean region, where current values are around 16 to 19 °C (Brun et al., 2022). The MAT method shows the greatest fluctuations over time across all sites. For precipitation (PANN), the signal from ODP site 976 mirrors the interglacial-glacial variations with large amplitudes for MAT and WAPLS (Fig. 8a). The other two sites show wetter conditions for certain glacial phases: MIS 34 at Monte San Giorgio, except for CAM, and MIS 36 at Palominas (Fig. 8b, c). In addition, MIS 35 appears drier at Palominas than at ODP 976 and Monte San Giorgio. MIS 34 appears quite marked at ODP site 976 and Monte San Giorgio with temperature decreases down to 9 °C depending on the method. However, precipitation is reversed at Monte San Giorgio with more rainfall during this cold period. Even though this phase seems to be the most intense, it remains close to the conditions obtained for MIS 36, which are dry at Monte San Giorgio and ODP site 976, but wet at Palominas. This observation is consistent with current knowledge presenting MIS 34 as a very cold phase (Hodell et al., 2023; Margari et al., 2023), but we cannot exclude the possibility that conditions were less extreme in other Mediterranean areas. CAM reconstruction shows higher precipitation and temperature fluctuations for Monte San Giorgio and Palominas than for ODP 976. Two hypotheses based on pollen data may explain this difference: (i) forest diversity at the ODP 976 site is lower than at the other sites as the signal is almost exclusively composed of deciduous Quercus; (ii) the slightly higher presence of relict taxa, among them Tsuga, Zelkova or Cathaya. The deciduous Quercus pollen type encompasses a large number of species that develop within a very wide temperature range (MAAT: 5 and 23 °C; Fauquette et al., 1998), which does not allow for effective constraints on climate reconstruction. The other temperate deciduous forest taxa, which are present within narrower temperature ranges, are not present in sufficient quantities to exceed the presence and abundance thresholds, and therefore to constrain the climate estimate. For the other two sites, Quercus also dominates the assemblage, but other taxa are sufficiently represented to constrain the model, such as Ulmus or Ostrya and Celtis for the Spanish site, and Carpinus for the Italian site. Monte San Giorgio also presents sufficient quantities of Arecaceae to better constrain the reconstruction.

The methods applied to the Palominas dataset (CAM, MAT, BRT, RF, and WAPLS) are consistently colder than the results provided by the Coexistence Approach (CA) (Altolaguirre et al., 2020). They rely on different modern databases: PALEOFLORA for the CA (https://www.neclime.de/external-data.html, last access: 3 July 2025), the database of modern pollen spectra used by Fauquette et al. (1998) for the CAM and the Eurasian and Mediterranean Current Surface Sample Database (EAPDB; Peyron et al., 1998, 2011, 2013) for the other methods. The difference between the CAM and CA is all the more striking since they are based on similar concept and parameterization choices (i.e., the removal of taxa present only once, Altolaguirre et al., 2020). However, unlike the CA, CAM takes into account not only the presence/absence criteria but also the pollen percentages of the taxa, that certainly refines the climatic estimates.

Figure 9 Climate reconstruction based on faunal assemblages from archaeological and palaeontological sites in south-eastern Spain. Palaeontological sites are represented by the image of a horse, and archaeological sites are represented by a stone tool. MAAT, MTCO and MTWA are in degrees and PANN in mm yr⁻¹.

6.5 Comparison with climate reconstruction from archaeological/palaeontological sites in South-western Spain

In southern Spain, archaeological faunal data is an additional proxy for climate reconstructions and comparisons with pollen-based climate reconstructions from Palominas and ODP 976. The faunal assemblages come from archaeological sites (Fuente Nueva 3, Barranco Leon and Cueva Victoria) and palaeontological sites (Quibas and El Chaparral) (Fig. 9). Faunal-based climate reconstructions follow three different methods in terms of concepts and types of fauna. All these sites have higher temperature estimates (Fig. 9) than the temperature reconstructions based on pollen assemblages (Fig. 8). Only the results of the Coexistence Approach from Palominas pollen data (Fig. 8) are in the same range of values (15.5 °C). On the other hand, pollen-based reconstructions of precipitation are more consistent with results from faunal data (Fig. 9). The WAPLS values in ODP 976 and Palominas (∼ 750 and ∼ 700 mm yr⁻¹) are close to those reconstructed at Fuente Nueva (738–753 mm yr⁻¹; Agustí et al., 2009; Blain et al., 2011, 2016), Barranco Léon (732–741 mm yr⁻¹; Agustí et al., 2009; Blain et al., 2011), Quibas (723 mm yr⁻¹; Agustí et al., 2009) and Cueva Victoria (716 mm yr⁻¹; Agustí et al., 2009). BRT values (∼ 475 mm yr⁻¹) in ODP 976 are close to the lowest reconstructed precipitation values at Quibas (417 mm yr⁻¹; Hernández Fernández et al., 2007). This probably suggests that climate and vegetation in the southern Mediterranean region were more strongly influenced by, and sensitive to, variations in rainfall and drought than temperatures. Underestimations of temperatures using pollen-based reconstruction methods raise questions about the methodology, the modern databases used, the functioning of each method, and the application of these methods to marine pollen records and to the ancient Quaternary period (except for CAM which was developed for Neogene and Early Quaternary climate reconstruction). In the case of the two marine sites ODP 976 and Monte San Giorgio, temperatures appear to be underestimated. However, this same trend is observed at the continental site of Palominas. Unlike faunal data, ODP 976 records plants in the Alboran Sea from a large area stretching from the coasts to the mountains of both Morocco and southern Spain. The reconstructed climate is therefore an average climate over a large geographical area and not a local one.

7 Conclusions

This high-resolution study of the ODP 976 site provides multiproxy reconstructions of climate change in the South-western Mediterranean at the end of the Early Pleistocene. The sequence continuously documents marine and continental environments during the EMPT, when climate cycles began to lengthen from 41 to 100 kyr. This multiproxy approach highlights the synchronicity between the continental and marine response to glacial-interglacial cycles in the Western Mediterranean between MIS 37 and MIS 31 (∼ 1235–1062 ka). Comparisons with other sequences from the Central and Western Mediterranean show major differences between the Iberian Peninsula and southern Italy. Results highlight a marked climatic gradient between the Western and Central Mediterranean during MIS 31, constraining forest development under the influence of the Atlantic climate system, compared to southern Italy, which was less affected by this dynamic. Clay mineralogy shows an abrupt change in smectite and I $/$ K signals, reflecting atmospheric and marine conditions in the Alboran Sea around 1140 ka, i.e., at the transition between MIS 35 and 34. Multi-method climate reconstructions based on foraminifera (MAT) and pollen assemblages (CAM, MAT, WAPLS, FR and BRT) show variations in temperature and precipitation mirroring climate cycles. However, some glacial periods appear wetter than some interglacial periods, depending on the method used. MIS 34 reconstructions appear more pronounced in terms of temperature and precipitation, compared to other glacial phases. Pollen-based reconstructions systematically show colder temperatures than those based on faunal assemblages from palaeontological and archaeological sites in southern Spain, which could be explained by local vs regional climate reconstruction. On the other hand, precipitation reconstruction seems more consistent for the two types of archives, suggesting that vegetation changes are driven more by precipitation than temperature fluctuations in this part of the Southern Mediterranean. These results enhance our understanding of the environmental changes that occurred in the Western Mediterranean during the period between MIS 37-31, providing land-sea information. This type of research must now be extended to the entire EMPT in order to better understand the climate changes that took place at a time of major modification of the Earth's climate system, and to comprehend the environmental context in which the first hominins settled in Western Europe.