Anatolia (Türkiye) during the late Eocene and the Eocene–Oligocene Transition: successive warming and cooling, aridification, and implications for the westward dispersal of Asian terrestrial mammals

The Eocene–Oligocene Transition (EOT dated at ∼ 34 Ma) represents one of the most significant climatic shifts of the Cenozoic, marking the transition from the last warmhouse state to a coolhouse state. This global cooling had major consequences for terrestrial ecosystems and was synchronous with the dispersal of numerous Asian mammalian clades towards western Europe. However, the terrestrial expression of the EOT exhibits strong regional heterogeneity questioning its role in establishing dispersal corridors associated with floral and faunal turnovers.

Here, we describe, date, and document the paleoenvironments of a continental sedimentary section from Balkanatolia, a biogeographic province corresponding to the present-day Balkans (NE Mediterranean region) and Anatolia (Türkiye). This region most likely functioned as a critical stepping stone for the dispersal of Asian mammals toward western Europe. Our sedimentary record represents a fluvio-lacustrine system exposed over a ∼ 200 m section located in Büyükteflek (Çiçekdağı/Kırşehir area, Türkiye), dated by magnetostratigraphy to the Priabonian and the lower Rupelian, including the Oi-1 glaciation (∼ 33.65 Ma). Clumped isotope analyses on pedogenic carbonates across our record show evidence for a Late Eocene Warming starting during the middle Priabonian (ca. 37 Ma), followed by a marked cooling event at the Eocene–Oligocene Glacial Maximum (EOGM). Stable isotopic data and sedimentary facies further indicate that this complete interval is associated with a long-term aridification trend, starting during the Late Eocene warming and culminating at the EOT. Our results provide the first quantitative record of late Eocene warming on land, and our temperature estimates for the earliest Oligocene cooling are consistent with other Eurasian clumped-isotope records. These temperature shifts and associated aridification steps may have acted as contributing drivers of the late Eocene decline of Balkanatolian endemic taxa and likely facilitated the westward expansion of Asia-derived mammals ultimately resulting in the colonization of western Europe.

1 Introduction

The Eocene–Oligocene Transition (EOT) is one of the main climatic events of the Cenozoic. It marks the transition between the last warmhouse state of the Eocene to a Coolhouse state in the Oligocene (Westerhold et al., 2020). The EOT is commonly interpreted as a phase of accelerated climatic and biotic change, linked to an abrupt global drop in temperatures (Coxall and Pearson, 2007). This drop is primarily explained by a decline in atmospheric CO₂ concentrations and changes in oceanic gateways, which, in combination with orbital forcing, acted as triggers for the first major glaciation of Antarctica (DeConto and Pollard, 2003). Within the deep-sea record, the EOT is defined by a homogeneous and marked positive excursion of δ¹⁸O from ∼ 34.4 to 33.6 Ma (Katz et al., 2008; Hutchinson et al., 2021). The EOT is followed by the Early Oligocene Glacial Maximum (EOGM), dated from 33.65 to ∼ 33.16 Ma, corresponding to the peak of glacial expansion of the Oligocene associated with the major growth of the Antarctic ice sheet (Miller et al., 1991; Zachos et al., 1996; Hutchinson et al., 2021)

While the EOT and the EOGM are well documented in marine records, their impact on land remains poorly understood. This lack of understanding is mainly due to the heterogeneous nature of terrestrial responses: studies conducted so far show contrasting results regarding changes in seasonality, temperatures, and their timing, or the absence of these changes (e.g., Pound and Salzmann, 2017). Pollen compilations have shown a heterogeneous response of vegetation, especially in Eurasia where ecosystems show no, little, or a gradual shift to flora adapted to colder and more arid conditions (Pound and Salzmann, 2017). Stable and clumped isotope analysis of continental carbonates (pedogenic and lacustrine carbonates) provide complementary insights into surface temperature and hydroclimate evolution, further highlighting this heterogeneity. In North and South America, stable and clumped isotope records from pedogenic carbonates indicate a cooling of temperatures and changes in seasonality, although with varying intensity and chronology (Kohn et al., 2015; Fan et al., 2018; Antoine et al., 2021; Meijer et al., 2025). In Eurasia, the response to the EOT is much more heterogeneous. In Western Europe, geochemical proxies for weathering intensity in paleosols and stable isotopes in pedogenic carbonates in the Ebro Basin (Spain) and on the Isle of Wight (UK) indicate no significant changes in Mean Annual Temperature (MAT) or in seasonality (Sheldon, 2009; Sheldon et al., 2012). Alternative investigations on lacustrine carbonates, however, suggest a progressive shift toward more arid conditions, as observed in the Vistrenque Basin (France; Semmani et al., 2024), with marked aridity peaks predating the EOT, such as in the Alès Basin (France; Lettéron et al., 2017). In addition, clumped isotope analyses of gastropod shells from the Hampshire Basin document a decrease in MAT by ∼ 4–6 °C across the EOT (Hren et al., 2013). Further east, in the central Asian Tajik Basin (Tajikistan), Wang et al. (2020) found evidence for increasing aridification beginning prior to the EOT based on carbon isotopes in pedogenic carbonates. Studies conducted in the Qaidam Basin (north-central Tibet) provide further evidence of dramatic climatic changes across the late Eocene and the early Oligocene. Kent-Corson et al. (2009) documented a decrease in oxygen isotope values prior to the EOT, interpreted as reflecting a decline in temperatures. By contrast, Sun et al. (2020) reported a positive shift in oxygen isotopes from carbonate mudstones spanning the late Eocene to the middle Oligocene, suggesting a long-term trend toward enhanced aridity. In the nearby Xining Basin (northeastern Tibet), Page et al. (2019) documented a temperature decrease of ∼ 20 °C occurring both prior to the EOT and during the EOGM. This decline is interpreted as resulting from both a shift in the seasonality of carbonate formation and a decrease in surface temperatures, estimated at ∼ 9 °C. In the Lanzhou Basin (northeasternmost Tibet), Li et al. (2016) reported a shift toward more positive oxygen isotope values in pedogenic carbonates after 33 Ma, suggesting the onset of aridification following the EOT.

This heterogeneous Eurasian response to global climate change at the EOT and EOGM has been traditionally linked to regional changes in atmospheric circulation impacting monsoonal rainfall in eastern Asia (Page et al., 2019), and to Paratethys sea retreat impacting specific areas of central Asia and the Middle East (Toumoulin et al., 2022). It has also been linked to growing topography, impacting atmospheric moisture sources and pathways, especially in Europe around the nascent Alps where moisture sources significantly shifted between 35 and 31 Ma (Kocsis et al., 2014).

The variable and regional climatic and vegetation responses in Eurasia contrasts with the dramatic faunal changes witnessed at the continental scale. In Mongolia, a major faunal turnover, known as the Mongolian Remodelling, is interpreted as a response to the climatic shift toward more arid conditions (Meng and McKenna, 1998). But the most important faunal turnover spanning the late Priabonian and the earliest Oligocene is known as the “Grande Coupure” and corresponds to the extinction of 55 % of European placental mammals (Weppe et al., 2023), coinciding with the arrival of Asia-derived mammals. Recent diversification modeling highlights the role of climate stress on western European taxa as a main driver for their decline, rather than competition with the arriving taxa (Weppe et al., 2023). Moreover, climate modeling suggests that the two proposed dispersal routes (Fig. 1), via high latitudes in north-central Asia (Becker, 2009; Tissier et al., 2018; Mennecart et al., 2021) and via mid latitudes through Balkanatolia (Grandi and Bona, 2017; Licht et al., 2022a) were both influenced by pCO₂ changes and orbital forcing during the latest Eocene. These factors affected biome stability along each route, limiting dispersal opportunities (Tardif et al., 2021). Collectively, these changes might have favored the dispersal of Asia-derived taxa towards western Europe during the EOT, though the ecosystem and climatic evolution is yet poorly documented in these two areas, where no terrestrial EOT record currently exists.

Figure 1 Paleogeographic map of Balkanatolia and the Neotethysian realm at 34–33.5 Ma, modified from Montheil et al. (2026). Arrows indicate the proposed dispersal routes of the “Grande Coupure” event. The studied location is shown with a red star.

This study provides the first magnetostratigraphically dated continental section from Balkanatolia that records a hydroclimatic response to the EOT. Our analyses combine sedimentological facies description with stable and clumped isotope analyses of pedogenic carbonates. We compare our hydroclimatic records with other Eurasian records to better characterize the regional climatic evolution around the EOT and assess how it could have impacted faunal dispersal occurring in the region.

2 Geological context

2.1 Balkanatolia

Balkanatolia is a semi-continuous strip of land extending from the Alpine region to the Lesser Caucasus (Licht et al., 2022a) (Fig. 2). It comprised various terranes of Laurasian (Tisza, Dacia, Rhodope, Strandja Massif and the Pontides) and Gondwanan (Alcapa, Greater Adria, Anatolide-Tauride, South Armenia) affinities that progressively accreted throughout the Cretaceous and the early Paleogene due to the northeastward motion of Afro-Arabia and Greater Adria and consequent closure of the Anatolian Neotethys and the Sava Ocean (Van Hinsbergen et al., 2020; Mueller et al., 2022; see Fig. 2B).

Figure 2 (A) Current structural overview of the Balkanatolian terranes. Terrane abbreviations: ST (Standja Massif), An-TD (Anatolide-Tauride), D (Dacia, SA South Armenia), LK (Luda Kamchiya Trough), MP (Moesian Platform), Pd (Pontides), R (Rhodope), ACP (Alcapa), T (Tisza), GA (Greater Adria). (B) Middle Eocene structural overview of Balkanatolia (based on Licht et al., 2022a, and updated by Montheil et al., 2025). The star indicates the studied locality of Büyükteflek. (C) Simplified geological map of the Çiçekdağı Basin, with the synthetic stratigraphic log and the different basin units. Maps in panels (A) and (B) are modified from Montheil et al. (2025); geological map in panel (C) modified from Gülyüz et al. (2013).

During the Paleocene and Eocene, Balkanatolia consisted of large islands isolated from Afro-Arabia, Western Europe and Eastern Asia by the Paratethys to the north and east, the Neotethys to the south, and the Dobrujan, Valaisan and Carpathian seaways to the northwest (Barrier et al., 2018; Palcu and Krijgsman, 2023; Montheil et al., 2026). This geographic isolation led to faunal endemism that prevailed until the Late Paleogene (Métais et al., 2018).

Tectonic activity and sea-level fluctuations created a dynamic landscape, often depicted as a mosaic of islands during the Paleocene and Early Eocene (Popov et al., 2004). These fluctuations were caused by alternating phases of extension, thrusting, and shortening, resulting in submergence and emergence of the terranes (Van Hinsbergen et al., 2020). Between the Lutetian and Bartonian, the retreat of shallow seaways separating Balkanatolian islands enhanced land connectivity (Licht et al., 2022a; Montheil et al., 2026). As a result of renewed crustal shortening and sea level fall, terrestrial continuity with Laurasia was progressively re-established from the late Middle Eocene to Early Oligocene: first in the east via Anatolia and the Cimmerian terranes (Barrier et al., 2018), then in the west through continental sediment infill of the Luda Kamchiya Trough, connecting to the Moesian Platform (Doglioni et al., 1996). The first arrival of Asia-derived taxa in Balkanatolia is dated to the Bartonian (40–38 Ma) and associated with the decline of Balkanatolian endemic taxa (Licht et al., 2022a). This first wave of dispersal from Asia did not reach western Europe and Balkanatolia acted as a cul de sac (Métais et al., 2023). The study of rodent assemblages of Balkanatolia shows that the Asian newcomers remained largely confined to Balkanatolia until the early Oligocene (Maridet et al., 2025). Only a few Asia-derived mammals colonized western Europe by passing through Balkanatolia during the late Priabonian (mammalian biohorizon MP18-MP20; ca. 37–33.9 Ma). These mammals include the anthracotheriid artiodactyl Elomeryx, possibly the gelocid artiodactyl Phaneromeryx, and the amphicyonid carnivoran Cynodictis (Métais et al., 2023). Faunal exchanges through Balkanatolia appear to have declined during the EOT and the EOGM in favor of the higher latitude routes, and resumed later during the Rupelian (see Fig. 1; Mennecart et al., 2021; Métais et al., 2023, Maridet et al., 2025).

2.2 The Çankırı Basin

Anatolia, as we know it today, is divided into the Pontides and the Anatolide-Taurides terranes, located in its northern and southern parts, respectively. The separation of these two terranes is marked by the Izmir-Ankara-Erzincan Suture (IAES). In central Anatolia, the Kırşehir Massif (also known as the Central Anatolian Crystalline Complex, CACC) is sandwiched between the Pontide and Anatolide-Taurides terranes. The Kırşehir Massif has been interpreted either as a former volcanic arc located between the Pontides and the Anatolide-Taurides or as the subducted and exhumed northern tip of the Taurides (Lefebvre et al., 2013; Van Hinsbergen et al., 2016). The exhumation of the CACC was followed by renewed contraction related to a Paleogene collision of the massif with the Pontides (Gülyüz et al., 2013).

The Çankırı Basin, at the northern edge of the CACC, is a sedimentary basin that can be subdivided into three main structural and stratigraphic units cropping out along its southern border, each reflecting distinct tectonic settings and depositional histories (Akgün et al., 2002; Kaymakcı et al., 2003) (Fig. 2C). The (1) Çiçekdağı Belt, represents the basin basement and belongs to the crystalline core of the CACC. This unit corresponds to Cretaceous volcanic arc units (Akgün et al., 2002) or is considered part of the Central Taurides (Lefebvre et al., 2013; Van Hinsbergen et al., 2020), which was accreted onto the Pontides during its collision with the Taurides. It is overlain by the (2) Çankırı Basin Fill, a Paleogene sequence filled with shallow marine to continental deposits (Akgün et al., 2002). This unit records alternating episodes of marine and lacustrine transgression and regression events throughout the Paleogene. The lower part consists of a Paleocene-Eocene regressive sequence, transitioning from flysch to molasse, and later overlained by nummulitic limestones. From the Middle Eocene to the Oligocene, the sea progressively retreated, giving way to predominantly more continental red clastics, intercalated with and eventually overlain by Oligocene evaporites (Kaymakcı et al., 2003). The third unit, referred to as the (3) cover series, overlies the Cretaceous–Paleogene deposits and began to accumulate during the Early Miocene with fluvio-lacustrine sediments alternating with evaporitic layers, continuing up to the Plio-Quaternary (Kaymakcı et al., 2003; Akgün et al., 2002).

During most of the Eocene, the basin was bounded by the sea: the northern part of Çankırı corresponded to a deep-sea environment, which was tectonically active due to its proximity to the IAES zone, while the southern part corresponded to a partly shallow marine to continental deltaic environment, which was relatively tectonically quiescent (Erdogan et al., 1996). During the late Eocene (Priabonian), marine deposition had ceased, and the basin had lost its connection to the sea, evolving into a large, closed fluvio-lacustrine system (Campbell et al., 2026).

Our study site is located in the Çiçekdağı area (Fig. 2C). In the area, the Çiçekdağı Belt crops out along the Çiçekdağı anticline, which separates the Çiçekdağı Syncline and the Yerköy Syncline, located to the south and north of the Çiçekdağı area, respectively, where Paleogene deposits of the Çankırı Basin fill are exposed (Gülyüz et al., 2013; Licht et al., 2022a).

The sedimentary sequence begins with the Yoncalı Formation, which contains conglomerates with volcanic clasts and lignites (Licht et al., 2022a). It is overlain by the Late Lutetian–Early Priabonian marine limestones of the Kocaçay Formation, which are rich in nummulites (Akgün et al., 2002; Gülyüz et al., 2013; Licht et al., 2022a). Above these beds, lie the continental deposits of the Incik Formation, where fossils of Brontotheriidae and Hyracodontidae (Perissodactyla) of Asian affinity have recently been discovered (Licht et al., 2022a). These are then overlain by lacustrine deposits including marls, gypsum beds, and tuffs, known as the Sekili Member of the Incik Formation (Fig. 2C).

2.3 The Büyükteflek section

Our study section is located within the Yerköy Syncline, near the village of Büyükteflek (see Figs. 2C, and S1 in the Supplement), and spans the Incik Formation and the Sekili Member. In this area, the Incik Formation is represented by red beds, sandstones, conglomerates, and pedogenic carbonates. In the upper part of the section, the base of the Sekili Member is marked by thick conglomeratic units and sandstone bodies containing reworked carbonate clasts and invertebrate fossils derived from the upper Kocaçay Formation, where caliche horizons are commonly developed (Licht et al., 2022a). This section has been previously studied by Gülyüz et al. (2013) and Licht et al. (2022a). The latter work provided improved constraints on the basal part of the section, assigning a Priabonian age based on marine microfossils (ostracod and planktonic foraminifera) assemblage within the uppermost beds of the Kocaçay Formation. The upper part of the section, covering the Sekili Member, was dated using U–Pb analyses of interbedded tuff layers, yielding ages of 32.1 ± 1.1 and 32.6 ± 1.1 Ma and was correlated with the reversed polarity chron C12r (early Rupelian) based on magnetostratigraphic data. The lower part of the section, extending from the Incik Formation to the base of the Sekili Member, was insufficiently sampled, preventing a precise assignment by magnetostratigraphy within the Priabonian and the early Rupelian. Given the age constraints for the lower (Priabonian) and the upper (Rupelian) parts, the section should cover the EOT.

3 Methods

3.1 Sedimentary logging

We measured and logged a 197 m section spanning the Incik Formation and the Sekili Member at ∼ 50 cm resolution, located at the Büyükteflek locality of Licht et al. (2022a) (see Fig. S1 for the overall picture of the Büyükteflek section). Our work consists of a detailed stratigraphic log starting at the level 148 m of Licht et al. (2022a), ending at the level 345 m. The aim was to improve the resolution of the log to characterize the sedimentary facies of this section. These facies were classified into 6 distinct lithofacies, adapted according to the classification of Miall (2013). Photos of the sedimentary facies are provided in Fig. 3, microphotographs and complementary photos are available in the Fig. S2.

Figure 3 Examples of sedimentary facies described in Table 1. (A) Well-developed caliche (facies Smp); hiatus 2, marked by red dotted lines, corresponds to the stratigraphic level at 292 m; hammer for scale. (B) Well-developed caliche (facies Fmp), eroded at the top by channelized conglomerates (facies Gmc); this surface corresponds to hiatus 1 at 263 m; hammer for scale. (C) Massive clayish red beds (facies Fmp); hammer for scale. (D) Massive greenish sandstones containing gravel-sized clasts up to 2 cm (facies Sm), overlying red beds (facies Fmp); knife for scale. (E) Clast-supported conglomerates overlying caliche developed in sandstones (facies Smp); hammer for scale. (F) Massive sandstone with diffuse carbonate development (facies Smp); knife for scale. (G) Alternation of matrix-supported (facies Gmm) clast-supported conglomerates (facies Gmc) and the channelized sandstones (facies St); Leny Montheil (180 cm) for scale. (H) Trough cross-bedded sandstones (facies St); hammer for scale.

3.2 X-Ray diffraction (XRD)

24 carbonate samples were collected all along the section and consist of nodules or blocks from continuous massive, carbonated layers. All the samples were reduced to powder with a Dremel, and filtered with a 160 µm sieve. The mineralogy of the carbonate samples was characterized using an X'Pert Pro MPD X-ray diffractometer (Malvern Panalytical), equipped with a cobalt X-ray tube (wavelength λ = 1.79 Å), operating at 40 kV and 40 mA. Detailed procedures can be found in Couturier et al. (2025). Sample powder was deposited onto low-background Silicium (Si) plates. The powders were homogenized using drops of ethanol and allowed to dry. During analysis, samples were rotated at 15 rpm to improve measurement quality.

XRD scans were conducted over a 2θ range from 34 to 38°, with a step size of 0.026°. Each scan lasted 10 min, allowing for the identification of carbonate mineral phases based on their primary diffraction peaks (see XRD patterns in Supplementary data 1).

3.3 Magnetostratigraphy

For the magnetostratigraphic study, 110 rock samples were collected using a portable electric drill (core diameter of 1 in., i.e. 2.54 cm) and oriented with a magnetic compass. Samples were retrieved from red clastic layers ranging from clay to silt grain size and from pedogenic carbonates (from nodules to well-developed caliches) developed on silty sandstones or red clastics. We also collected nine oriented blocks in clay rocks when they were too soft to be drilled.

Natural Remanent Magnetization (NRM) and either its stepwise alternating field (AF) demagnetization with field increments of 5 mT from 0 to 70 mT or its stepwise thermal (TH) demagnetization with temperature increments of 20, 30, 40, or 50 °C from room temperature (∼ 20 °C) up to 660 °C were measured with the Superconducting Rock Magnetometer SRM760R (2G enterprises) of CEREGE (Aix-Marseille University). Remanent signal from some clay samples was resistant to alternating field (AF) demagnetization and consequently required thermal demagnetization.

89 samples were demagnetized using AF demagnetization, and 16 samples were subjected to thermal demagnetization of which 14 samples were analyzed using both methods. Data were plotted in orthogonal vector diagrams (Zijderveld, 1967) and stereographic projections and interpreted using the PuffinPlot software (Lurcock and Wilson, 2012). The Characteristic Remanent Magnetization (ChRM) directions were then calculated using the principal component analysis of Kirschvink (1980). The dataset of this study has been supplemented with the paleomagnetic data from Licht et al. (2022a; 60 samples). The sample levels were adjusted to match our stratigraphic log using field-related descriptions.

Samples exhibiting low magnetic intensity or lacking a coherent directional trend, which likely reflects local noise rather than a primary signal, were excluded from the dataset. The ChRM for each sample was isolated using progressive thermal or alternating field (AF) demagnetization, performed in at least four steps covering the ranges 120–660 °C or 5–70 mT (see demagnetization steps, e.g. samples BT_82, BT_48_2 and BT28_2 Fig. 5A, E, F; see Table S1 in the Supplement and Supplementary data 2 and 3), respectively, to ensure stable endpoint determination. Reliability was assessed based on the stability of demagnetization trajectories (ensuring the absence of erratic behaviour or secondary components) and a maximum Maximum Angular Deviation (MAD) of 15° (Fisher, 1953; McFadden and McElhinny, 1990). Polarity consistency was further validated through reversal tests, reversal angle calculations, and Fisher direction statistics, performed using custom Python scripts and the https://paleomagnetism.org/ (last access: 20 May 2026) platform (Koymans et al., 2016; Fisher, 1953; Arason and Levi, 2010; Deenen et al., 2011; Heslop et al., 2023; McFadden and McElhinny, 1990; King, 1955; Tauxe, 2010; Tauxe et al., 2010; Tauxe and Watson, 1994; Tauxe et al., 2008; Zijderveld, 1967). These tests were applied to both our data and those from Licht et al. (2022a), and 28 samples were subsequently classified as outliers. Given the intrinsic uncertainties in paleomagnetic measurements, which can exceed ±30° in Virtual Geomagnetic Pole (VGP) latitude for equatorial sites due to shallowing and low inclination resolution (Tauxe et al., 2010), VGPs with latitudes between −30 and 30° were classified as uncertain, and no polarity interpretation was attempted for these cases.

3.4 Stable and clumped Isotopes

23 carbonate sample powders identified as pure calcite from XRD, and one composed of approximately equal proportions of calcite and dolomite, were selected for stable isotope analysis. The samples were analysed on a Thermo 253+ mass spectrometer attached to a Thermo Kiel IV device, run under the Qtegra™ software and in Dual Inlet mode at CEREGE. Sample powder corresponding to around 80–120 µg of pure carbonate was loaded into glass vials and placed on a carousel along with international standards (IAEA-603, NBS18) and laboratory standards (BDH and VIA-1, which are made of pure calcite) in the Kiel oven held at 70 °C. After each vial is in turn put under vacuum, the powder reacts with 2 to 3 drops of 105 % phosphoric acid at 70° to produce CO₂ for analysis. The released CO₂ is then trapped and cleaned via a series of cryo-traps, and finally sent to the Thermo 253+ mass spectrometer for measurements of the oxygen and carbon isotope values. Isotopic values are reported relative to IAEA-603 international standards in the V-PDB (Vienna Peedee Belemnite) reference materials.

14 of these samples were selected for clumped isotope analysis based on their carbonate content and their distribution along our stratigraphic log. They were treated overnight with 3 %–8 % NaOCl, rinsed and centrifuged 4 times and dried with a freeze-drier, to remove nitrate contamination (Fiebig et al., 2024). The samples were then analysed at CEREGE on the same previously described system. Details about the clumped isotope procedure at CEREGE are given in Supplement Sect. S1. Briefly, primary standards (ETH-1 to -4), secondary standards and unknown sample powders were loaded into the glass vials and placed on the carousel of the Kiel oven still held at 70 °C. The procedure is similar to the one described above for stable isotope measurement with the exception of an additional step of cryotrapping involving a trap filled with a Porapaq molecular sieve, held at −20 °C by a Peltier system. The resulting purified CO₂ gas is transferred to the 253+, where extra Faraday cups will allow the measurement of clumped isotopes ($m/z$ 47, 48, 49) on top of carbon and oxygen isotopes ($m/z$ 44, 45, 46) and pressure baseline corrections ($m/z$ 47.5 and 48.5).

Baseline-corrected δ data are screened with an in-house Matlab script for statistical outliers and are then translated into δ¹⁸O, δ¹³C and Δ₄₇ projected to the Intercarb-Carbon Dioxide Equilibrium Scale (I-CDES90; Bernasconi et al., 2021) using D47crunch (Daëron, 2021), pooling over data from all sessions of the dataset. In total, we run 10 to 20 replicates per unknown sample spread over 28 sessions, with an inter-session repeatability of 30 ppm in Δ₄₇ for all replicates (unknown and standards). Paleotemperatures are calculated using the unified Δ₄₇-T calibration of Anderson et al. (2021). Uncertainties around temperature estimates are provided with and without propagating intersession-related uncertainty (difference of ±1 to 2 °C at 2 s).

Note that to test the reproducibility of our results and the impact of the NaOCl pre-treatment, three samples were divided into two aliquots: one with pre-treatment, and one without. Two out of three samples show no statistical difference with or without pre-treatment (difference < 2 s); for a third sample, pre-treated and non-treated aliquots lie within 3 s of each other (see Table S2).

Data presented in the manuscript are only from pre-treated samples. Raw data for unknown samples and standards are provided together with interpreted data, to allow possible future reprocessing of Δ₄₇ values and paleo-temperature estimates, in Table S2.

4 Results

4.1 Sedimentology and XRD analysis of the Büyükteflek section

The lower part of the studied section (from ∼ 100 to 263 m) is characterized by three distinct sedimentary facies (Table 1). The first facies (Fmp; see Fig. 3B, C, D) consists of clayey to silty red beds, generally massive in structure, and occasionally exhibiting discontinuous planar laminations. These deposits may contain root traces and mottling due to reducing conditions. Above level 240 m, they display well developed pedogenic features, including small carbonate nodules to well-developed caliche horizons. X-ray diffraction (XRD) analysis indicates that all carbonate samples from these beds are composed exclusively of calcite. The second facies (St; see Fig. 3G, H) is composed of channelized, trough cross-bedded sandstones, with grain-size ranging from very fine to medium sand. Some beds contain thin layers of gravel or pebbles. The fossiliferous horizon described by Licht et al. (2022a), containing Brontotheriidae (Embolotherium aff. andrewsi) and Hyracodontidae (Prohyracodon sp.), occurs at approximately the 215 m level within these facies. The third facies (Sm; see Fig. 3D) consists of structureless greenish medium to coarse sandstones, containing isolated clasts ranging in size from gravel to approximately 2 cm. The three facies Fmp-St-Sm are grouped under Facies Association FA1 (Table 2).

Table 1 Sedimentary facies of the Büyükteflek section.

Table 2 Facies associations of the Büyükteflek section.

The level at 263 m is marked by a field-identifiable surface, corresponding to a channelized erosional surface developed above a well-developed caliche. This hiatus was interpreted by Licht et al. (2022a) as the boundary between the Incik Formation and the Sekili Member (hiatus 1 on Fig. 4). From 263 to 292 m, the facies consist of channelized sandstones (St) and pedogenized clay to silty red beds (Fmp). This interval marks the appearance of three new facies. The first one (Smp; see Fig. 3A, E, F) consists of reddish to pinkish, massive, very fine to fine grained sands. These beds are strongly pedogenized and feature carbonate nodules, well-developed caliche horizons or diffuse carbonate. The second facies (Gmm; see Fig. 3G), which is coarser, is composed of massive, matrix-supported conglomerates with a sandy matrix, containing clasts ranging from 1 mm to 5–10 cm and exhibiting weak normal grading. Some of these conglomerates include reworked Nummulites, probably derived from underlying marine deposits. The last facies (Gmc; see Fig. 3B, E, G) consists of clast-supported conglomerates with a sandy matrix, containing imbricated clasts ranging from 1 mm to 5–20 cm, and may also contain reworked Nummulites. These three new facies, together with the channelized sandstone Smp-St-Gmm-Gmc are grouped under Facies Association FA2 (Table 2). X-ray diffraction analysis indicates that the nodules found at these depths contain only calcite.

Figure 4 Stratigraphic log of the Büyükteflek section (see Sect. S2 in the supplement for the exact location of hiatuses and tuffs). Colour coding illustrates the different facies and facies associations. The VGP latitude diagram shows the paleomagnetic data and the corresponding magnetozones (R1, N1, N2, R2, N3, N4).

At 292 m, a second unconformity is observed, marked by the development of a 2–4 m thick caliche (hiatus 2; see Fig. 4). This paleosol is the thickest of the section and forms a topographic feature that can be tracked over hundreds of meters (Fig. 3A). Five facies are present from this unconformity up to the top of the logged section in this study (level 355 m). These include the same channelized sandstones (St), along with pedogenized red beds assigned to the Fmp facies. The Smp facies is also identified, which is particularly abundant in this part of the section. X-ray diffraction analysis shows that carbonate in these beds may contain either calcite (from 292 to 325 m) or dolomite (from 302 to 355 m). In addition, the two conglomerate facies, Gmm and Gmc, are present within this part of the section. This interval of the log includes both facies associations FA1 and FA2.

4.2 Magnetostratigraphy of the Büyükteflek section

The magnetostratigraphic data presented here complement the work of Licht et al. (2022a) on the same section. Natural Remanent Magnetization (NRM) measurements reveal that clay samples generally exhibited weaker intensities, ranging from 6.9 × 10⁻⁴ to 1.5 × 10⁻¹ A m⁻¹, while coarser samples, from silt to fine sand, displayed higher intensities, ranging from 1.5 × 10⁻² to 8.9 × 10⁻¹ A m⁻¹. Some samples show a low-coercity secondary of a component in both fine and coarse grain sizes (e.g., BT82, Fig. 5A). For the silt to fine sand fraction, demagnetization typically shows a stable decay toward the origin (e.g., BT23, Fig. 5B). Similarly, several clay samples also display a consistent decay toward the origin, though often accompanied by a higher degree of noise (e.g., BT71, BT80, Fig. 5C, D).

Figure 5 Representative intensity plots and orthogonal demagnetization diagrams obtained after Alternating Field (AF) and Thermal demagnetization.

Thermal demagnetization results suggest the presence of both magnetite and hematite, with unblocking temperature ranges (performed from 120 to 660 °C) supporting their coexistence. Zijderveld diagrams for these samples typically show a single-component decay, indicating the presence of a stable Characteristic Remanent Magnetization (ChRM) (see representative example in Fig. 5). More samples recorded a normal geomagnetic polarity than a reverse geomagnetic polarity, consistent with the investigated time interval, which is characterized by a predominance of long-lasting normal polarity chrons. To complete the geomagnetic polarity record, we added the magnetostratigraphic dataset from Licht et al. (2022a), which covers the upper part of the section and includes a higher proportion of reversed magnetized samples attributed to Chron C12r.

Regarding ChRM directional statistics (see Fig. S5), the inclination derived from geographic coordinates (53.2°) is consistent with the latitude of the site, whereas it drops to 26.5° after tectonic correction, which is inconsistent with the expected value for the site. Notably, the α95 value of inclination data shows minimal change (from 3.2 to 4.7°) between geographic and tectonically corrected coordinates, indicating no significant tectonic influence. The inclination discrepancy, most likely attributable to the shallowing effect commonly observed in laminated lake sediments (e.g. Philippe et al., 2023), thus reflects a systematic bias rather than increased data scatter.

The observed angle between the mean normal and reversed (corrected) directions is γ = 21.0°, whereas the critical angle at the 95 % confidence level is γc = 13°, resulting in a formal failure of the reversal test. The datasets are strongly unbalanced (N = 85 vs. N = 21), and the reversed population shows greater dispersion (k = 7.51; α95 = 12.4°) than the normal population (k = 17.49; α95 = 3.8°; see Fig. S5). This statistical failure likely reflects the smaller sample size and higher scatter of the reversed dataset rather than a true departure from geomagnetic antipodality, and the angular difference between the means remains moderate. Given that this study is based on continental sedimentary rocks, such angular differences between mean directions remain within acceptable limits (Tauxe et al., 2010). Fisher statistics indicate a mean inclination of 28.7° for the normal samples and 34.7° for the reversed (corrected) samples (see Fig. S5). These relatively low and comparable values further support the presence of inclination shallowing. Consequently, the tilting observed in the section indicates that the magnetization was acquired prior to tilting rather than post-tilting.

The primary character of the magnetization is further supported by the presence of several magnetic reversals along the section. Furthermore, samples containing both magnetite and hematite exhibit a single-component magnetization, suggesting that hematite recorded the ambient magnetic field during or shortly after deposition, rather than during later mineralogical alteration, and thus supporting a primary nature for the magnetization.

Four distinct magnetic chrons have been identified within the Sekili Member interval of the studied section. These chrons are defined based on the presence of at least five consecutive, reliable samples exhibiting consistent magnetic polarity. In contrast, the lower part of the section, corresponding to the Incik Formation, displays less well-defined magnetic chrons, possibly reflecting either lower sampling resolution and/or enhanced diagenetic overprinting. However, the dominance of normal polarity intervals in this portion of the section led us to interpret it as representing a normal chron.

4.3 Stable and Clumped isotopes

The stable isotope composition of calcite samples ranges from −6.5 ‰ to −4.9 ‰ for δ¹⁸O (on average: −5.7 ‰ V-PDB) and from −7.8 ‰ to −5.2 ‰ for δ¹³C (‰ V-PDB). Calcite δ¹⁸O values do not exhibit any clear stratigraphic trend throughout the section. In contrast, δ¹³C values for calcite are slightly more variable. Below hiatus 1, δ¹³C values average −6.8 ± 0.3 ‰. Between hiatus 1 and hiatus 2, these values increase slightly, with a mean of −6 ± 0.6 ‰. Above hiatus 2, δ¹³C values display a decreasing trend, ranging from −5.9 ‰ to −7.8 ‰ with an average value of −6.9 ± 0.8 ‰. Among the samples above hiatus 2, one contains both calcite and dolomite in approximately equal proportions. Given its comparable clumped isotope values to those of other samples, it was retained in the results and discussion (see Table S3). Clumped isotope-derived temperatures (Δ₄₇) from pedogenic calcite carbonates fall within the range of Earth's surface temperatures (< 40 °C). The Δ₄₇ results can be divided into three distinct groups. The first group (all samples below hiatus 1, but the very last sample; 4 samples in total) shows relatively stable values, with a mean temperature of 24.4 °C ± 0.7 (2 s). A significant increase in surface temperature is observed in the second group, corresponding to the last sample below hiatus 1 and all the samples situated between hiatus 1 and hiatus 2 (5 samples in total). This interval yields a higher and relatively consistent mean temperature of 34 °C ± 3.5 (2 s), representing a ∼ 9 °C increase compared to the interval below hiatus 1. Then the third group shows a decreasing trend in temperatures. The third group (5 samples), located above hiatus 2, displays decreasing temperatures of ∼ 7–8 °C. This group yields a mean Δ₄₇-derived temperature of 27.4 °C ±  3.6 (2 s), with respectively the highest to the lowest values from 31.9 ± 5 to 21.6 ± 4.2 °C. Clumped isotope results from the sample with a mixture of dolomite and calcite (see Table S2 sample BTCARB16) sample yield similar Δ₄₇ values than for pure calcite samples of the same interval, suggesting that the CO₂ extracted during the acid reaction at 70 °C is prominently calcitic (e.g. Li and Guan, 2023). The exclusion of the sample from our dataset does not alter the average temperature of the interval (e.g. 27.5 ± 4.1 2 s without the sample instead of 27.4 ± 3.6 2 s). Student's t tests were performed between the first and second groups, and between the second and third groups. The three groups show statistically significant differences. The comparison between Group 1 and Group 2 yielded a p value of 0.005 and an effect size (Cohen's d) of $\left|d\right|$ = 3.33, indicating a strong statistical difference between these two groups. Similarly, the comparison between Group 2 and Group 3 produced a t test p value of 0.030 and an effect size of $\left|d\right|$ = 1.668, also suggesting a substantial difference between the two groups (see Table S4).

Based on clumped and stable isotope data, δ¹⁸O values of soil water (δ¹⁸O_water) values (V-SMOW) have been reconstructed according to Kim and O'Neil (1997) recalculation (see Fig. 6). These values can be divided into three groups, consistent with the grouping defined from clumped isotope temperatures. Group 1 displays a mean δ¹⁸O_water value of −3.6 ‰ (V-SMOW). A marked increase is observed at the final sample just before hiatus 1, with Group 2 showing a mean value of −1.9 ‰. This is followed by a slight decrease in Group 3, which exhibits a mean δ¹⁸O_water value of −2.7 ‰ (see Table S4).

Figure 6 Scaled stratigraphic log highlighting the pedogenic carbonate sampling interval. From left to right are plotted δ¹⁸O, δ¹³C, δ¹⁸O_water, and Δ₄₇ (±1 se) values from pedogenic carbonates. Coloured fields (light green, light orange, light blue) indicate distinct groups, interpreted as reflecting differences separated by the hiatuses. Horizontal grey dotted lines delineate the hiatus 1 and hiatus 2.

5 Discussion

5.1 Refining the Age of the Büyükteflek Section Across the Eocene–Oligocene Transition

The paleomagnetic record of the Büyükteflek section was established by integrating the chrons identified here with the magnetostratigraphic framework published by Licht et al. (2022a). This former work was correlated to the Geological Paleomagnetic Time Scale through biostratigraphy and U-Pb constraints. A total of six distinct chrons were recognized, comprising four normal polarity intervals and two reversed polarity intervals. Magnetozones were identified and correlated with the Paleogene period time scale of Speijer et al. (2020) from the GTS 2020.

Following the work of Licht et al. (2022a), we initially correlated the uppermost reversed-polarity magnetozone (R1) with Chron C12r (33.21–30.98 Ma), and the underlying normal polarity magnetozone (N1) with Chron C13n (33.73–33.31 Ma) based on the age of two volcanic layers dated in the section at 355 and 420 m. The limestone beds beneath our lowermost paleomagnetic sample have been dated by biostratigraphy to the Priabonian stage (Licht et al., 2022a), providing a maximum age of 37.71 Ma. This dating is strengthened by the sandstones underlying the limestone beds, which yield a U-Pb maximum depositional age of 38.7 ± 0.5 Ma (2σ) (Licht et al., 2022a). The overlying strata, which include the normal polarity magnetozone N3, therefore postdate this age. This interval contains mostly normal polarity samples, with only a few isolated reversed ones, suggesting that magnetic reversals may not have been preserved due to erosion, non-deposition, or limited sampling resolution. Several correlation scenarios were considered using these constraints (see Figs. S3, S4 and Table S5). Before correlating our magnetozones, we first considered all the deposits below the second hiatus as a single, continuous normal magnetozone N1bis (thus excluding the reverse samples defining R2). If correlated to chron C13n (as N1), this approach yields a minimum sedimentation rate of 58.6 cm ka⁻¹ for N1bis, which is unrealistically high compared to the upper part of the section and thus geologically implausible. Consequently, these rates do not support an Oligocene age for the underlying section and ensure that the section below hiatus 2 is partly Priabonian. Below hiatus 2, the reverse polarity magnetozone R2 has to be correlated to one of the five reverse chrons between C13r and the base of the Priabonian (C16r to C13r). The correlation for R2 that yields the most consistent accumulation rates when compared with the rates of the higher part of the section (4.5 cm kyr⁻¹ for R1) is with Chron C16n.1r (35.77–35.72 Ma; accumulation rate: 5.4 cm kyr⁻¹ for R2). We acknowledge that other hypotheses of correlation are possible (see Figs. S3 and S4 for a second and third hypothesis), but all other correlations yield much lower accumulation rates (e.g. < 1.8 cm kyr⁻¹ for R2). The magnetozone N3 below the reverse magnetozone R2 could not be confidently correlated with a specific chron due to the scarcity of reversed-polarity samples, suggesting that one or more magnetic reversals may have gone unrecorded or unsampled in this lower interval (see Fig. 7). Consequently, the interval extending from the marine fossil assemblage of Priabonian age (ca. 60 m; see Fig. 7) up to the top of N3 should fall within the interval between the base of the Priabonian and the base of Chron C16n.1r (e.g. within the 37.71–35.77 Ma window).

Figure 7 Stratigraphic log of the Büyükteflek section showing the magnetostratigraphic correlation with associated accumulation rate (R1; R2) and minimum accumulation rates (N1; N2; N3; N4), calibrated with the Paleogene Geomagnetic Polarity Time Scale (Gradstein et al., 2020).

This magnetostratigraphic framework allows for temporal subdivision of the Büyükteflek section. The basal unit, below the first unconformity, corresponds to the early to middle Priabonian. The layer at 215 m, which yielded fossils of Brontotheriidae (Embolotherium aff. andrewsi) and Hyracodontidae (Prohyracodon sp.), belong to magnetozone N3; this magnetostratigraphic position constrains the age of the fossil-bearing horizon to the 37.71 to 35.77 Ma window. The interval between the two unconformities is constrained by the correlation of magnetozone R2 to chron C16n.1r; it likely spans Chrons C16n.1n, C16n.1r and parts of the middle Priabonian (see Fig. 7). The upper part of the section, above the second unconformity, is attributed to chrons C13n and C12r. Chron C13n coincides with the Eocene Oligocene Glacial Maximum (EOGM; Hutchinson et al., 2021), while the uppermost part corresponds to the transition to the later, warmer times of the Rupelian.

5.2 Evolution of the depositional system and origin of the unconformities

Previous work by Gülyüz et al. (2013) suggested that episodes of folding, thrusting, and uplift along the Çiçekdağı anticline, located directly to the south of our section, began between ∼ 38 and 35 Ma. They also identified a reversal of the direction of paleocurrents from predominantly southward to northward flow direction within the Sekili Member, suggesting a local reorganization of paleodrainages. It is thus very likely that the two unconformities recorded in our section are at least partly related to on-going tectonic uplift and deformation of the southern margin of the basin. Our magnetostratigraphic correlation does not allow us to determine the exact duration of the first unconformity (hiatus 1), which spans part of the early-mid Priabonian. The second unconformity (hiatus 2) spans part of the early C13n, the whole of C13r, C15n, C15r, and the upper part of the C16n.1n. This corresponds to a minimum interval ranging from 35.6 to 33.8 Ma, including the latest Priabonian and most of the EOT (34.44 to 33.65 Ma, sensu Hutchinson et al., 2021).

The first unconformity is associated with a marked shift in prominent sedimentary facies and is associated with a change of depositional system. Deposits below hiatus 1 only display one facies association (FA1), which is interpreted as reflecting fine-grained floodplain or alluvial plain deposits (Table 1; Miall, 2010; Leeder, 2011). The deposits above hiatus 1 are dominated by coarser, more energetic, and more pedogenized/carbonated facies grouped in Facies Association 2 (FA2), which is interpreted as a streamflow-dominated alluvial fan/delta with subaqueous portions extending into the lake. This more subaqueous or distal part of the fan-delta is supported by the presence of the facies Smp (see Fig. 4 and Tables 1, 2; Miall, 2010; Renaut and Gierlowski-Kordesch, 2010).

According to Gülyüz et al. (2013), the local reversal in paleocurrents around that time corresponds to a shift from southward to northward and is interpreted as marking the uplift of the Çiçekdağı anticline directly south of our section. This uplift is corroborated by the reworking of nummulitic limestones, starting with hiatus 1. We thus suggest that the change of sedimentary facies relates to a change in the fluvial network and the regional hydrosystem. This reorganization is consistent with the uplift of the Çiçekdağı anticline and the associated re-routing of regional drainage systems (Gülyüz et al., 2013; Licht et al., 2022a). During the early parts of the Priabonian, the area was the locus of a southward-oriented, medium-distance drainage network exiting to the southern edge of the Çankırı basin, where the youngest marine deposits are found (Licht et al., 2022a). Hiatus 1 marks the uplift of the southern margin of the basin, the disconnection of this river network and the setup of a local, proximal-sourced northward-oriented drainage exiting to a newly formed lacustrine system occupying the Çankırı Basin.

Pedogenesis and carbonate cementation are progressively observed into coarser sediments through the section. At the base, pedogenic features of the Fmp facies (Table 1) are developed within clayey deposits in a floodplain environment. In the upper part of the succession, beginning above hiatus 1, pedogenesis and carbonate cementation are instead mainly associated with coarser gravity-flow deposits, represented by sandstones of the Smp facies (Table 1). These sandstones accumulated in the distal portions of a fan-delta system, near the lake margin or even within its subaqueous domain. Owing to this transitional position, such deposits were highly sensitive to lake-level fluctuations. The occurrence of well-developed caliches or carbonate-cemented horizons within these sandstones likely indicates phases of lake retreat, during which pedogenesis could develop on exposed surfaces.

The second unconformity is not associated with any major change of sedimentary facies but the appearance of dolomite. From the hiatus 2 onward, almost all of our carbonate-rich sandy facies (Smp) consist predominantly of dolomite. When not related to diagenetic processes, dolomite can precipitate in hypersaline environments such as lagoons or saline lakes (García-Ruiz, 2023) depending on local salinity and pH fluctuations (Kim et al., 2023). We thus interpret the appearance of dolomite as reflecting an increase of lake salinity and a shift toward higher evaporation relative to precipitation, where water loss concentrates dissolved salts faster than they are diluted by freshwater inputs (Guo et al., 2023). This interpretation is further supported by the occurrence of Oligocene evaporites in other parts of the Çankırı Basin (Kaymakcı et al., 2003; Karadenizli, 2011). Together, these observations provide multiple independent sedimentological indicators of lake-level fluctuations under increasingly arid conditions.

The synchronicity between the second unconformity and the latest Priabonian and the EOT, a period of documented increase of aridity in many places in Eurasia (Abels et al., 2011; Page et al., 2019; Lettéron et al., 2022; Li et al., 2016; Sun et al., 2020; Semmani et al., 2024), also supports this interpretation. The sedimentary hiatus at the second unconformity was thus likely the combination of a favorable tectonic regime and a period of climate-controlled lake-level drop. Interestingly, the peak of lake retreat starts during the latest Priabonian and covers the EOT, and thus predates the EOGM and the global low of eustatic level (Miller et al., 2008); lake levels increase again during the EOGM. The lake level drop at the second unconformity seems thus uncorrelated with glacio-eustatic levels, and rather related to an earlier, late Priabonian and EOT aridity crisis.

5.3 Preservation of the primary Δ₄₇ signal

Potential diagenetic processes, including burial recrystallization or hydrothermal alteration, could affect isotopic values in pedogenic carbonates. However, thin section observations (Fig. S2) show well-preserved micritic to microsparitic textures in all analyzed pedogenic carbonates. These carbonates occur as in situ nodules and caliche horizons within paleosols (facies Smp) or as diffuse cements within floodplain deposits (facies Fmp, Gmm), with no evidence of coarse sparite recrystallization, blocky calcite cement, or hydrothermal mineral assemblages, features that would typically accompany significant burial diagenesis or hydrothermal overprinting (Henkes et al., 2014; Stolper and Eiler, 2015). Moreover, the Δ₄₇-derived temperatures are coherent with formation in the vadose zone, where temperature variations reflect near-surface soil conditions (see Fig. 6; Quade et al., 2013), further supporting an in-situ origin. Together with the consistency of our temperature trends with other Eurasian clumped isotope records (e.g. Page et al., 2019; Semmani et al., 2024), this petrographic, stratigraphic, and geochemical evidence supports the primary nature of the Δ₄₇ signal and the robustness of the reconstructed paleotemperatures.

5.4 A first record of the Eocene–Oligocene Transition in Anatolia

Our Δ₄₇-derived temperatures indicate a drop of ∼ 7 °C in soil carbonate growth temperature across the second unconformity, spanning the late Priabonian and the EOT, and persisting through the onset of the EOGM. This temperature drop is similar to what is observed in several other clumped isotope records of the EOT world-wide. A comparable cooling of ∼ 7 °C in surface temperatures across the Eocene–Oligocene Transition was documented by Fan et al. (2018) in eastern Wyoming, inferred by clumped isotopes on pedogenic carbonates. In Western Europe, Δ₄₇ measurements on freshwater gastropods from the Hampshire Basin (England) reveal a ∼ 4–6 °C decrease (Hren et al., 2013). In northeast Tibet, analyses from the Xining Basin indicate an apparent drop of ∼ 20 °C, which has been interpreted as resulting from both a shift in the carbonate growth season and an actual temperature decline, leading to an estimated ∼ 9 °C decrease in surface temperatures (Page et al., 2019; see Table S6 for additional continental temperatures across the EOT.

Studies have shown that modern soil carbonates typically form during periods of soil dewatering during warm and dry seasons, and are commonly found in mid-latitude seasonal climate (Methner et al., 2016). This leads to a Δ₄₇ bias toward Warm Month Mean Temperatures (WMMT) rather than MAT (Methner et al., 2016; Hough et al., 2014; Passey et al., 2010; Peters et al., 2013; Quade et al., 2013). In some instances, particularly under sub-humid climates with alternating wet and dry cycles, carbonate formation may instead reflect values closer to the mean annual temperature (Breecker et al., 2009; Peters et al., 2013). Under humid tropical monsoonal regimes, the bias can even shift toward recording winter temperatures, more closely reflecting Cold Month Mean Temperatures (CMT; Licht et al., 2022b). We compared our reconstructed surface temperatures with Eocene temperature estimates from pollen assemblages in Anatolia. There are no paleotemperature estimates for the Bartonian and Priabonian of central Anatolia, but Raynaud et al. (2026) report paleotemperatures for the Lutetian (48 to 41 Ma), stage that is globally warmer than the later stages of the Eocene (Westerhold et al., 2020). Lutetian Mean Annual Temperatures (MAT) range from 17.2 to 21.9 °C, while Warm Month Mean Temperatures (WMMT) range from 24.7 to 28.1 °C. Our Δ₄₇-derived temperatures thus likely reflect carbonate growth during the warm season (see Fig. 8). We thus interpret the drop of ∼ 7 °C in our record as a drop in summer temperatures.

Figure 8 Clumped isotope (Δ₄₇) paleotemperatures from pedogenic carbonates in the Büyükteflek section (Çiçekdağı Basin, Turkey), spanning the Early to Late Priabonian (∼ 37.8–35 Ma) to the Early Rupelian (∼ 33.4 Ma). The orange line shows the moving average; error bars = 1 standard error (SE). Temperature estimates for the Lutetian (48–41 Ma), derived from paleobotanical data reported by Raynaud et al. (2026). (WMT: Warm Month Temperatures; MAT = Mean Annual Temperatures; CMT = Cold Month Temperatures)

Stable isotope data from pedogenic carbonates reveal no clear change in δ¹⁸O values across the second unconformity, and a slight decrease of approximately 0.7 ‰ is observed in δ¹³C values. Reconstructed δ¹⁸O_water values exhibit marked variability both before and after the second unconformity, ranging from −1 ‰ to −4 ‰. These values contrast with previous reconstructions for the same region and time interval, which yielded consistently lower values. For example, Licht et al. (2017) reported δ¹⁸O_water values between −5 ‰ and −7 ‰ from late Lutetian lacustrine carbonates, while Robert et al. (2026) obtained values of −5.3 ‰ to −8 ‰ from late Lutetian Balkanatolian mammal tooth enamel. The latter are also consistent with modern δ¹⁸O values of meteoric waters in coastal Mediterranean Anatolia (Schemmel et al., 2013), after applying a ∼ 1 ‰ correction for differences in seawater isotopic composition during the Eocene (Tindall et al., 2010), and with Oligocene meteoric water estimates (Lüdecke et al., 2013). The discrepancy between our reconstructed δ¹⁸O_water values and those reported in earlier studies may reflect a proxy-related bias. Reconstructions are expected to capture primary meteoric water signatures, yet the δ¹⁸O_water is influenced both by meteoric input and by evaporation. As emphasized by Kelson et al. (2023), evaporation in arid environments can enrich soil water δ¹⁸O_water by 1 ‰–5 ‰, even at depths > 40 cm. Such an effect could account for the offset between our estimates and previous reconstructions, and can explain the great variability in reconstructed soil water δ¹⁸O_water values.

The slight decrease in δ¹⁸O_water values recorded across the unconformity, from mean values of −1.9 ‰ to −2.7 ‰ is thus difficult to interpret considering the significant imprint of evaporative effects in our record. We note however that a ∼ 7 °C air temperature drop would likely decrease rainwater isotopic composition by ca. −2 ‰ at these latitudes, and could thus partly explain this decline (Bowen, 2008).

Our record thus indicates that the Eocene–Oligocene Transition and the EOGM in central Anatolia are associated with summer temperature cooling of ∼ 7 °C and an increase of aridity marked first by a temporary lake retreat during the EOT, followed by a lake rewatering during the EOGM, though marked with a persistent increase in lake salinity compared to the late Priabonian lake system, as evidenced by the appearance of dolomite beds.

5.5 A Late Eocene Warming?

The portion of our section surrounding the first hiatus exhibits more pronounced shifts in isotopic data. Clumped isotope-derived surface temperature estimates reveal a ∼ 9 °C increase beginning with the last paleosol below hiatus 1, and warm temperatures persist within all paleosols between hiatus 1 and 2. Previous studies (e.g., Bohaty and Zachos, 2003) have documented a “Late Eocene Warming” around ∼ 37–36 Ma seen in benthic foraminiferal records from the Southern Hemisphere. Their reconstructions suggest a ∼ 3–5 °C temperature increase under ice-free conditions. The more recent benthic foraminiferal compilation of Westerhold et al. (2020) dates this warming episode to ∼ 37.5 to ∼ 36.5 Ma. Tremblin et al. (2016) suggest a significant warming of approximately 4 °C in equatorial and North Atlantic regions starting during this time window and lasting until at least 35 Ma, and possibly up to the EOT. The Recent pCO₂ reconstruction of CenCO2PIP (The Cenozoic CO₂ Proxy Integration Project (CenCO₂PIP) Consortium et al., 2023) suggests an increase of pCO₂ during the latest Eocene before the EOT, though the amplitude and chronology of this increase are still poorly constrained. By contrast, Tremblin et al. (2016) suggests that this warming episode could result from a reorganization of ocean circulation in the North Atlantic. Numerous paleoceanographic records in the North Atlantic suggest an onset of a regional subtropical gyre or an early version of the Atlantic meridional overturning circulation (AMOC) around 36 Ma, which could lead to regional warming (Coxall et al., 2018; Liu et al., 2018).

Other records of the late Eocene warming on land are sparse. A recent study made in the Paris Basin records a decrease of δ¹⁸O values on bulk carbonate, also interpreted as a record of the late Eocene Warming during the C16r chron (Le Callonnec et al., 2025). Clumped isotope data from Tibetan soil carbonates highlight a warming trend across the late Priabonian, culminating between 35 and 34 Ma (Page et al., 2019), and starting after a significant aridification step near the top of chron C17n.1n (Abels et al., 2011), coherent with the beginning of the warming trend observed here.

In our record, this warming is accompanied by a mean δ¹⁸O_water increase of 1.7 ‰, along with enhanced variability. Given the proximity of the study area to both the Neotethys and Paratethys domains, major changes in moisture source would require substantial paleogeographic reorganization, and can be thus disregarded as a cause for the δ¹⁸O_water change (Kayseri-Özer, 2013). The observed increase in δ¹⁸O_water values is more consistent with progressive evaporative enrichment under increasingly arid conditions, as frequently seen in soil carbonates (Kelson et al., 2023). The pronounced variability in δ¹⁸O_water values is also consistent with increased soil evaporation. This change is associated with a moderate rise in δ¹³C values of approximately 0.7 ‰. This shift could be explained by a slight increase in plant water stress, also pointing to a trend toward increased aridity (Cerling and Quade, 1993; Kohn, 2010).

The synchronicity of this apparent increase in temperature and aridity with the first unconformity suggests that it might have partly controlled the chronology of the hiatus. Though the uplift of the basin's southern margin is likely the primary cause for the discordance, the unconformity might have been enhanced by a base-level drop related to an increase of aridity. Overall, our record highlights an important climatic transition during the middle Priabonian and preceding the EOT, marking the onset of an irreversible shift toward more arid conditions, coupled with elevated surface temperatures. This event can be hypothetically attributed to a “late Eocene Warming” impacting Eurasia, though its chronology and origin remain to be more precisely constrained. The subsequent cooling during the EOT brought back clumped isotope temperatures to their earlier, early Priabonian values, but accentuated a long-term trajectory toward increasing aridity.

5.6 Regional climatic response and implications for Asian terrestrial fauna dispersal

As many previous Eurasian records, our terrestrial section highlights the existence of a profound aridity crisis during the Priabonian, marked by two aridification steps during the middle (hiatus 1) and latest (hiatus 2) Priabonian. We show that the first aridification step is associated with a significant warming event, that we attribute to the Late Eocene Warming. The first step likely initiated prior to the base of the Chron C16n.1r (35.77 Ma).

In Western Europe, a comparable aridification phase is observed: Lettéron et al. (2022) describe an increase in aridity within the ASCI lake system (Alès, Saint-Chaptes, and Issirac; France) between 37.1 and 35.2 Ma (Lettéron et al., 2018). Other records from Central Asia indicate a similar increase in aridity during this interval. Wang et al. (2020) document enhanced aridity in the Tajik Basin, inferred from pedogenic carbonates showing reduced plant-derived soil respiration, correlated with the retreat of the Paratethys Sea during chron C17n.1n. Studies from the Xining Basin likewise record an early phase of aridification occurring during chron C17 (Page et al., 2019; Yang et al., 2022).

The link between this mid Priabonian aridity crisis and the late Eocene warming event is yet difficult to draw. Warming in the North Atlantic related to a proto-AMOC or a regional gyre, the only current explanation for the late Eocene warming should lead to more regional evaporation and more precipitation on land (Elsworth et al., 2017). It thus cannot explain a Eurasian-wide aridity crisis. A tectonically driven Paratethys sea retreat episode in central Asia, as suggested by Wang et al. (2020), could partly explain a regional increase in aridity concomitant with the warmer temperatures of the late Eocene Warming, though its farflung effects up to Balkanatolia and Europe are difficult to predict.

Regardless of its origin, the warming phase observed in our record covers mammalian biohorizons MP18 and MP19 and the first arrival of selected Asia-derived taxa in western Europe such as anthracotheriids, amphicyonids, and possibly gelocids, all of which likely dispersed through Balkanatolia (Métais et al., 2023). This synchronicity suggests that the late Eocene warming could have favored the westward dispersal of several Asia-derived taxa, which were already present on Balkanatolia since the Bartonian (Licht et al., 2022a), and remained endemic to this biogeographic province during the early Oligocene (Van de Weerd et al., 2023; Maridet et al., 2025). The environmental stress associated with the warming event and the aridity crisis could have exacerbated the decline of endemic Balkanatolian taxa and further favored the penetration of Asia-derived taxa (Raynaud et al., 2026). The quasi-complete disappearance of Balkanatolian endemic taxa during the Priabonian (Licht et al., 2022a) further supports this interpretation, albeit Priabonian Balkanatolian faunas remain poorly documented (fewer than ten Priabonian sites; Licht et al., 2022a). We note however that the MP18 and MP19 interval corresponds to the peak of endemic mammalian diversity for large clades in western Europe like artiodactyls (Weppe et al., 2023). The late Eocene warming seems thus to have contrasted impacts: in Balkanatolia, it may have stressed endemic taxa and facilitated Asia-derived dispersal, whereas in western Europe, endemic diversity peaked before collapsing at the Grande-Coupure. The second aridity step is synchronous with the latest Priabonian and EOT and starts in our record at ∼ 35.2 Ma. It is associated with a ∼ 7 °C cooling in carbonate growth temperatures leading to the EOGM, which resembles other terrestrial clumped isotope records covering this time interval in Eurasia (Page et al., 2019; Hren et al., 2013), and is compatible in amplitude with the rare European paleobotanical records that show some sensitivity to the EOT (Teodoridis and Kvaček, 2015; Suc et al., 2025). It is yet unclear if cooling started with the beginning of the second aridity step associated with the lake retreat at ∼ 35.2 Ma, or, more likely, with the end of the event through the end of the EOT and early EOGM, just before the lake rewatering.

Interestingly, the second aridity step predates the EOGM by at least 1.5 Myr, suggesting that its origin is not directly related to the Antarctic glaciation. Rather, we suggest that it reflects the hydroclimatic response to an earlier cooling event in the North Atlantic. Indeed, cooling in the North Atlantic seems to have preceded the EOT and to have been heterogeneously distributed both spatially and temporally, starting at least 1 Myr before the Eocene–Oligocene boundary (Śliwińska et al., 2023). This heterogeneity is so far explained by a complex response of the North Atlantic gyres to decreasing pCO₂, which changes the location, strength and structure of these gyres (Śliwińska et al., 2023). Regardless of its origin, the second aridity step covers mammalian biohorizons MP20 and the onset of biodiversity decline for western European endemic artiodactyls, which culminates in the earliest Oligocene with a decline of 77 % in artiodactyl diversity (Weppe et al., 2023). We thus suggest that the environmental stress associated with the second aridity step could have enhanced the decline of western European taxa. While it is unclear how this environmental stress impacted Balkanatolian fauna due to the sparse paleontological record of the Priabonian and Rupelian in the area, the last known occurrence of an endemic Balkanatolian taxon (the embrithopod Axainamasia sandersi, Métais et al., 2024) is recorded around the EOT (Sanders et al., 2014) and suggests that local ecological or paleogeographic factors allowed some endemic Balkanatolian species to persist well beyond the colonization of Balkanatolia by Asian predators and competitors. The aridification step and the consecutive EOGM correspond to periods without faunal exchanges between Balkanatolia and western Europe (Mennecart et al., 2021; Métais et al., 2023). This timing is counter-intuitive because the eustatic drop at the EOGM fully connected western Balkanatolia to western Europe (Licht et al., 2022a; Montheil et al., 2025). This suggests that this climatic event resulted in the set-up of a persistent environmental barrier in central and southern Europe, blocking the westward expansion of Asia-derived mammals. This paleoenvironmental barrier may have been enhanced by orographic effects on regional precipitation, caused by the nascent Alps (Kocsis et al., 2014).

6 Conclusion

Our sedimentological and stable isotope records document two major phases of aridification during the late Eocene. The first, occurring prior to the base of chron C16n.1r (35.77 Ma), is associated with a marked increase of ∼ 9 °C in summer surface temperatures recorded by clumped isotope data, representing the first evidence of the Late Eocene Warming in western Eurasia. The second phase took place during the latest Priabonian, extending through the EOT and into the EOGM, and is characterized by a subsequent decline of ∼ 7 °C in summer surface temperatures during the EOGM.

These climatic transitions likely reflect large-scale reorganizations in atmospheric and oceanic circulation in the North Atlantic preceding the EOGM. Such environmental changes may have contributed to the decline of endemic fauna in Balkanatolia and western Europe during the late Eocene. However, the precise chronology and extent of these events require further refinement through integration with additional continental records to strengthen this interpretation.