Evolution of continental temperature seasonality from the Eocene greenhouse to the Oligocene icehouse –a model–data comparison

. At the junction of greenhouse and icehouse climate states, the Eocene–Oligocene Transition (EOT) is a key moment in Cenozoic climate history. While it is associated with severe extinctions and biodiversity turnovers on land, the role of terrestrial climate evolution remains poorly resolved, especially the associated changes in seasonality. Some paleobotanical and geochemical continental records in parts of the Northern Hemisphere suggest the EOT is associated with a marked cooling in winter, leading to the development of more pronounced

Temperature, MATR) in parts of the Northern Hemisphere. However, this increase of the annual temperature range has been poorly studied by climate models; uncertainties remain about the geographical extent of this phenomenon and the associated climatic processes. Although other components of the climate system vary seasonally (e.g., precipitation, wind), we therefore focus on the seasonality of temperatures only. 20 In order to better understand and describe temperature seasonality patterns from the middle Eocene to the early Oligocene, we use the Earth System Model IPSL-CM5A2 and a set of simulations reconstructing the EOT through three major climate forcings: pCO2 decrease (1120/840 to 560 ppm), the Antarctic ice-sheet (AIS) formation, and the associated sea-level decrease (-70 m).
Our simulations suggest that seasonality changes across the EOT rely on the combined effects of the different tested 25 mechanisms which result in zonal to regional climate responses. Broad continental areas of increased MATR reflect a strengthening of seasonality (from 4°C to > 10°C increase of the MATR) across the EOT in agreement with MATR and seasonality increase is recorded in the decline of species characteristic of warm paratropical to temperate environments such as conifers Doliostrobus sp. (conifers), Nypa sp. (Palms), Rhodomyrtophyllum sp. (Myrtaceae), and some families with 65 tropical elements such as Annonaceae, Lauraceae, Cornaceae, Flacourtiaceae, Icacinaceae, Menispermaceae, and, depending on bioclimatic zones, the expansion of temperate to boreal vegetation through the increase of deciduous and / or coniferous species (Eldrett et al., 2009;Kunzmann et al., 2016;Kvaček, 2010;Kvaček et al., 2014;Mosbrugger et al., 2005;Utescher et al., 2015;Wolfe, 1992). These vegetation changes are associated with a decrease of the Coldest Month Mean Temperatures (CMMT) across the EOT (-2.2 to -10.3°C) and starting earlier at some localities, during the mid-to late Eocene (Eldrett et 70 al., 2009;Moraweck et al., 2019;Mosbrugger et al., 2005;Tanrattana et al., 2020;Tosal et al., 2019;Utescher et al., 2015;Wolfe, 1994). Isotopic analyses have also documented this seasonality increase in different continental localities between the Priabonian (37.8 to 33.9 Ma) and the Rupelian (33.9 to 27.82 Ma; Grimes et al., 2005;Hren et al., 2013;Zanazzi et al., 2015). Some of these changes are not directly quantifiable (e.g. the reduction of gastropods growing season length, United Kingdom; Hren et al., 2013), others are particularly strong, with MATR increase up to 15.6°C (Canada; Zanazzi et al., 75 2015). Seasonality increase is also documented for shallow waters of the Gulf of Mexico (increase of the MATR; Ivany et al., 2000;Wade et al., 2012). Seasonality strengthening and associated vegetation changes have been proposed as driving factors for dental morphological changes for European grazing perissodactyls (Joomun et al., 2010). Some studies have suggested a link between increased temperature seasonality and latitude (e.g. Eldrett et al., 2009;Zanazzi et al., 2007Zanazzi et al., , 2015 but data seem insufficient to validate this relationship and this trend has not been confirmed by recent palynological 80 compilation (Pound and Salzmann, 2017).
Southeast Europe and Anatolia show generally weaker and heterogeneous changes in temperature seasonality, with either no seasonality changes, slight seasonality lowering or slight seasonality strengthening from the mid Eocene to the Rupelian (Bozukov et al., 2009;Kayseri-Özer, 2013). This variable signal has been linked to a strong marine influence on this insular part of Eocene Europe (Kayseri-Özer, 2013). Conversely, North and East Asia temperature seasonality evolution is more 85 comparable to western Europe and North America trends (Quan et al., 2012;Utescher et al., 2015). Vegetation changes reflect an increase of the seasonal temperature range, mainly through the EOT (MATR increase of 2 to 2.5°C; CMMT decrease of ~2.2°C (Quan et al., 2012;Utescher et al., 2015). The appearance of tubers in lotus (Nelumbo sp.) during the Eocene suggests the establishment of a dormant phase in these plants and thus, of a period unfavorable to plant growth (Li et al., 2014). Fossils showing these structures have been described in Southern China (Hainan Province) and in the extreme 90 East of Russia (Kamchatka Peninsula); they could be favored by Eocene cooling and increased seasonality on the East Asian continent (Budantsev, 1997).
In the Southern Hemisphere, studies of Paleogene localities are rarer. Despite a record of late Eocene cooling in Australia, New Zealand and Patagonia, independent proxies (stable isotopes on teeth, bones and pedogenic carbonates, paleobotanical reconstructions) do not suggest a marked temperature seasonality during the Eocene Nott and Owen, 95 1992;Pocknall, 1989). In Australia, the presence of more pronounced wood rings suggests a late Paleogene increase in https://doi.org/10.5194/cp-2021-27 Preprint. Discussion started: 1 April 2021 c Author(s) 2021. CC BY 4.0 License. seasonality starting in the mid-Oligocene at the earliest (Bishop and Bamber, 1985;Nott and Owen, 1992). Finally, the environmental and climatic impact of the EOT in Africa remains virtually undocumented (Hutchinson et al., 2021).

Former model work
A first attempt to explain temperature seasonality change was made by Eldrett et al. (2009). In their palynological and 100 modelling study, Eldrett and coauthors explained high latitude (Greenland) seasonality strengthening by pCO2 drop and the consequent increase of sea-ice formation over the Arctic Ocean. Sea-ice extension induces a strong albedo feedback in their experiment, which results in a large decrease of atmospheric temperature during the winter. In addition, Li et al. (2018) showed an important contribution of geography changes to seasonality strengthening during the Oligocene. Apart from these two studies, EOT modelling experiments mostly focused on oceans through mean annual temperature changes (Baatsen et 105 al., 2020;Goldner et al., 2014;Hutchinson et al., 2018Hutchinson et al., , see 2020Kennedy et al., 2015;Kennedy-Asser et al., 2019, 2020Ladant et al., 2014b). The recent model intercomparison of Hutchinson et al. (2021) has shown a reasonable agreement between 4X to 2X modelling experiments and proxy-data atmospheric and surface ocean temperature reconstructions. The different simulations of the model intercomparison support once more the primary role of pCO2 drop. Hutchinson et al. (2021), also quantified the general effect of AIS formation and contemporary paleogeographic changes (opening of Southern 110 Ocean gateways, Antarctic geography or global geography). EOT modelling experiments yield mixed answers regarding the temperature feedback resulting from both AIS and paleogeographic changes (Goldner et al., 2014;Kennedy et al., 2015;Ladant et al., 2014aLadant et al., , 2014b. Among these models, only one tested the combined effect of both AIS and global paleogeographic changes (Kennedy et al., 2015) and suggested a moderate impact of AIS on global climate sensitivity, as previously suggested by other modelling work (Goldner et al., 2013;Kennedy et al., 2015). 115

Aim of the study
This study is an attempt to reconstruct temperature seasonality changes described by paleoenvironmental and paleoclimatic data between the middle Eocene and the early Oligocene. It includes both a synthesis of indicators of seasonality changes (Table S1) and paleoclimatic modelling. We reconstruct global patterns of temperature seasonality change and their main mechanisms through a set of five simulations taking into account the three-major climate forcings described through this 120 time interval: pCO 2 drawdown, AIS formation and concomitant sea-level lowering.

Material and methods
We used the IPSL-CM5A2 general circulation model, which is built upon the CMIP5 Earth system model developed at the Institut Pierre-Simon Laplace (IPSL), IPSL-CM5A-LR , Sepulchre et al., 2019. The IPSL-CM5A-LR Earth system model is composed of the LMDZ atmospheric model (Hourdin et al., 2013), the ORCHIDEE land surface and 125 vegetation model (Krinner et al., 2005), and the NEMO v3.6 ocean model which includes modules for ocean dynamics https://doi.org/10.5194/cp-2021-27 Preprint.  Madec, 2008), biochemistry (PISCES, Aumont et al. 2015) and sea-ice (LIM2, Fichefet and Morales-Maqueda, 1997). The atmospheric grid has a horizontal resolution of 3.75° longitude per 1.875° latitude (96×95 grid points), and is divided into 39 vertical levels. For a more detailed description of the model and its different components, the reader is referred to Sepulchre et al. (2020). 130 Five simulations were carried out to reconstruct the evolution of temperature seasonality from the middle Eocene to the early Oligocene (Table 1). The applied 40 Ma paleogeography framework is the map developed by Poblete et al. (2021) and already used in Tardif et al. (2020) and Toumoulin et al. (2020). It features common late Eocene geography characteristics such as: an open Panama Seaway, an open Tethys with submerged Arabic peninsula, a strongly maritime Europe, a Turgai Strait between Europe and Asia, and a narrow Southern Ocean gateways ( Figure 1). The orbital parameters were set to 135 preindustrial values and the solar constant was reduced accordingly to its Eocene value (1360.19 W/m 2 ; Gough, 1981).
Vegetation was implemented as a boundary condition, using a zonal band of PFTs using modern vegetation distribution patterns.  (Ladant et al., 2014b); % Land: total land surface (millions of km 2 : 10 6 km 2 ); MAT: Mean Annual global 2-meter air Temperature (°C); SST: Sea Surface Temperatures (°C). Simulations with an asterisk constitute realistic middle Eocene (Lutetian / Bartonian) and early Oligocene (Rupelian) scenarios, others are either sensitivity experiments (2X, 2X-ICE) or use the high value of the pCO2 range estimated for the time interval (4X). 145

Simulation
Simulations were compared in pairs to highlight differences between the late Eocene and the early Oligocene. The simulation set is composed of both realistic and idealized experiments (Table 1). Simulations 4X, 3X and 2X represent most of the pCO2 range described from the mid-Eocene (Lutetian) to the early Oligocene (Rupelian; Foster et al., 2017). These covered Antarctica on climate, while 2X-ICE-SL constitutes a realistic representation of the early Oligocene climate. In these experiments, the Antarctic ice cap was set to 32.5.10 6 km 3 according to Ladant et al. (2014b). The 70-m sea-level drop was defined following eustatic drop estimates for the EOT (Coxall et al., 2005;Katz et al., 2008;Lear et al., 2008). It is responsible for important geography changes related to an increase of land proportion, such as the emergence of the Arabic Peninsula and the retreat of the proto-Paratethys epicontinental sea. 160 All simulations are run for 4000 years when temperatures indicate a quasi-equilibrium with only negligible temperature drifts within the global mean ocean (< 0.1°C/century; Figure S1). These trends are consistent with most model studies and do not affect the quality of atmospheric change described in this study (e.g. Hutchinson et al., 2018;Lunt et al., 2016). The results considered here are averages of the last 100 years of the model runs.

Data compilation 170
Simulation results were compared to MATR changes (∆MATR) documented by proxy-data records (Table S1). We compiled published MATR and CMMT proxy-data from various research fields: paleobotany (macrofossils and palynology), geochemistry (isotopic measurements on various material), and paleontology. The data were selected to range from the Lutetian (47.8 Ma) to the end of the Rupelian (27.8 Ma). This large time interval allows the representation of seasonal temperature changes parallel to the long-term cooling of the Eocene. The inclusion of data from the middle Eocene allows a 175 comparison with simulations testing the effect of a pCO2 lowering alone, before AIS formation at the EOT. It is justified by the presence of paleobotanical records suggesting a strengthening of the seasons already from the Lutetian to the Priabonian (e.g. Li et al., 2014;Mosbrugger et al., 2005). Compiled MATR Eocene-Oligocene changes correspond either to the values given in original publications, when they were available and precise, or to values re-calculated from the original data. For https://doi.org/10.5194/cp-2021-27 Preprint. Discussion started: 1 April 2021 c Author(s) 2021. CC BY 4.0 License.
publications for which ∆MATR were recalculated, we proceeded by minimizing the paleogeographic distance between sites 180 (especially in terms of latitudes) and checked that the values obtained were consistent with the authors original interpretation of the paleoenvironmental context. Some previously published seasonality increases were not here associated with changes in MATR, because they were either (1) estimated from other parameters, such as the length of the growing season, which do not allow the calculation of the MATR (Hren et al., 2013), or (2) derived from qualitative data that cannot be specifically associated with temperature values 185 (e.g., organism morphological changes such as teeth shape or plant tuber appearance; Joomun et al., 2010;Li et al., 2014).
These sites are displayed on the maps, but are not analyzed quantitatively. In order to better estimate the impact of changes in temperature seasonality, the length of the plant growing season (i.e., the number of months with an average temperature above 10°C) was recalculated using the formula of Grein et al., (2013) for data from the approximate coexistence analyses (Table s1). Paleocoordinates for every locality were reconstructed using an online service of Gplate, the plate tectonic 190 reconstruction models of Matthews et al. (2016) and a 40 Ma geography.

Comparison of model and data ∆MATR
Different analyses were made to evaluate the data-model agreement for temperature seasonality changes from the Priabonian to Rupelian (Table 2). First, a general agreement percentage was calculated from the direction of seasonality changes alone.
For this metric, model predictions are considered "good" for an individual site if they have the same sign as the data (i.e., 195 positive model ∆MATR at the location of a data point showing seasonality increase). For data indicating null ∆MATR, a good agreement was considered with model values ranging from -0.5 to 0.5°C. This first metric enables us to assert the agreement between our simulation and qualitative data.
In addition, Priabonian to Rupelian seasonality changes were then compared to the model predictions, by (1) assessing their correlation and (2) calculating the root of the mean squared distance between their values. These two analyses were 200 performed using Rstudio software (version 4.0.3; R CoreTeam, 2020, Boston, USA). Given the limited number of quantitative Priabonian-Rupelian data (n = 15), the statistical correlation of data-model ∆MATR was assessed from average ∆MATR with the non-parametric Spearman rank test. Modelled MATR values were extracted from a 3° longitude by 3° latitude area surrounding each data locality. This test returns a coefficient, ρ (rho), which is a measure of the correlation (ranging from -1 to 1 indicating total negative and positive linear correlations respectively and 0 the absence of correlation). 205 The statistical significance of the correlation is described through the p-value. In this analysis, we used the common significance level, α, of 0.05 (i.e., p-values < 0.05 indicate significant correlations).
Data-model agreement for ∆MATR from Priabonian to Rupelian data was also assessed through the use of the Root Mean Squared Estimate (see Kennedy-Asser et al., 2020, and their figure S1 for a detailed presentation of the method). This method consists in calculating the root of the mean squared distance between model and data values for comparable points. 210 Conversely to the Spearman Rank Test for which mean ∆MATR estimations were used, the distance is here measured using the full range of estimates at each data locality (i.e., minimum and maximum ∆MATR). Note that because it considers the https://doi.org/10.5194/cp-2021-27 Preprint. Discussion started: 1 April 2021 c Author(s) 2021. CC BY 4.0 License. full range of ∆MATR, this method tends to underestimate the difference between model and data. The lower bound of modelled ∆MATR at each locality was calculated as the difference between the lower MATR value over the 3° x 3° area centered around the locality for an Oligocene-like cold simulation (2X, 2X-ICE or 2X-ICE-SL) and the highest MATR value 215 over the 3° x 3° area for an Eocene simulation (4X or 3X). For the upper bound, we used the difference between the higher MATR value over the same area for an Oligocene-like cold simulation and the lower MATR value for an Eocene simulation.
The RMSE adjusted to ∆MATR is written as follows: Where MATRdiff-data and MATRdiff-model are MATR Priabonian to Rupelian changes estimated by data and model respectively, 220 and n the total number of localities.  Atmospheric temperatures decrease strongly in our different experiments (Table 1). The halving of pCO2 alone (i.e. without AIS formation and sea level drop) induces a global cooling of 5.8°C and 5.0°C for the air temperature and the surface ocean respectively. Temperature changes are characterized by polar amplification, with a stronger winter cooling at high-latitudes ( Figure 2. a,b,e,f,i,j). This phenomenon is well explained by the combined effect of albedo and sea-ice feedback. The pCO2 drop enables sea-ice growth over the Arctic and, to a lesser extent, the subsistence of snow on the ground during the cold 235 season, which increases winter surface albedo (Figure 3). In addition, seasonal sea-ice expansion limits ocean-to-air heat transfer at highest Northern latitudes and contributes to further winter cooling.

MATR increase zones
Along with its effect on temperatures, pCO2 decrease induces zonal ∆MATR including (1) an increase of MATR at highlatitudes, (2) a decrease of MATR across most mid-latitudes, and (3)  preferential lowering of winter temperatures aforementioned. Furthermore, the areas of colder winters and broader MATR become larger as pCO2 decreases: MATR increases from 60°N poleward between 4X and 3X, to 50°N poleward between 4X and 2X (Figures 4 and S3). In contrast, Antarctica shows moderate changes (3°C from 3X to 2X) compared to high Northern latitudes (6°C MATR increase from 3X to 2X). This is because ice-free Antarctica is characterized by high variability of 245 albedo as the continent oscillates from snow-free to snow-covered soil within a year, resulting in highly variable intra-annual temperatures ( Figure 4).

250
Figure 4: Mean Annual Temperature Range changes (°C). Shadings are model differences calculated with a Student t-test over the last 100 years of comparative simulations (95% confidence); white areas indicate no significant change in MATR. Symbols correspond to ∆MATR from proxy-data for different time steps: Priabonian-Lutetian (triangles); Rupelian-Lutetian (deltas); Rupelian-Bartonian (circles); Rupelian-Priabonian (diamonds). Orange symbols indicate qualitative values describing a temperature seasonality increase. In case of proxies reconstructing a range of equally probable values (e.g., 255 Coexistence approach), values shown are mean values. References for every data site are displayed on Figure 7 and available in the data compilation provided in Table S1.

MATR decrease zones
Areas with decreased seasonality are characterized by summer cooling that exceeds winter cooling, which reduces the 260 MATR (Figure 2 and 4). The widest zones with decreasing MATR are located within the 30-50°N latitudinal band, across North America, Western Europe, Central Asia, and 30-50°S for South America and Australia (Figure 4). To a lesser extent, seasonality decreases in Amazonia, equatorial Africa and India (Figure 4). This phenomenon is well explained by two distinct chains of reaction. First, the pCO2 drop flattens the P-E (precipitation minus evaporation, or net precipitation) latitudinal gradient caused by a general increase of P-E values between 50°S and 50°N, notably in North America and East 265 Asia (especially in part of India, Mongolia and SE Asia), and an overall P-E decrease at high latitudes in both hemispheres  (Figure 5 g,h). Secondly, changes in ocean temperature 270 have a strong impact on continental MATR in several areas at mid latitudes, such as in Patagonia and Europe. In these areas, cold air masses originating from the nearest oceans (North Atlantic and South Pacific) are brought by westerly winds during summer. They favor a lowering of temperatures during this season and thus, a lowering of the MATR (Figure 6). These cold air masses result from the regional increase in low-level cloud cover during summer (5 to 15% higher cloud fraction between 40-60°; Figure 6 d,h, S4) which increases planetary albedo. In addition, the presence of sea-ice at the North Pole (Figure 3  275 b,c) limits heat loss via the atmosphere during winter and results in a greater SST cooling during summer which significantly impacts European MATR ( Figure S5).

Antarctic Ice Sheet alone
AIS formation is responsible for a supplementary 0.9°C and 0.3°C cooling of the air temperature and the surface ocean 280 respectively (Table 1). It induces striking mean annual temperature differences, reaching -35°C above the continent (Figures   2 and S6). Its effect on atmospheric temperature varies regionally and is more significant over the Southern Ocean and Australia. In contrast with Arctic sea-ice which increases seasonality at highest northern latitudes, the AIS decreases southern latitude temperature seasonality (Figures 4 and S6). Indeed, simulations with the AIS have year-long white soil and high stable surface albedo, which reduces seasonal temperature variability (Figure 3.g-l).

Sea-level drop
In addition, sea-level decrease is responsible for a 1.0°C mean annual global 2m-air cooling (0.7°C for surface oceans) and results in considerable regional temperature changes in areas with important land-sea distribution changes (Table 1 and Figure S6). The increased proportion of land over ocean enhances latent and sensible heat fluxes during boreal winter, which reduces CMMT and conversely reduces energy fluxes during boreal summer, which increases summer warmth (Figures 2,  300 S3 and S6). The numerous areas of winter cooling and summer warming resulting from land fraction increase provoke a wide spatial extension of the MATR reinforcement zone (Figure 4 d). The zonal distribution of seasonal temperature changes generated by pCO2 drop and AIS formation is dismantled. East Africa, Western Asia and Russia are the most impacted areas, following proto-paratethys retreat and the emergence of the Arabic Peninsula ( Figure 1). More moderate seasonality changes are also visible as a result of sea retreats of smaller extent, such as the emergence of the Florida platform and the 305 modification of the East Asian coastlines (Figure 4 d).  Table 2: RMSE for ∆MATR between the Priabonian and the Rupelian. The line "%" give the percentage of sites where the direction of ∆MATR is adequately modelled (e.g., the model described a MATR reinforcement in the zone where data 310 indicate MATR increase). Modelled MATR estimates were considered as positive when > 0.5, negative when < 0.5, or null when ranging [-0.5; 0.5]. "rho" indicates the strength of the correlation estimated with the Spearman Rank test, with associated p-values ("p"). Significant correlations are indicated in bold and with **.

pCO2 lowering
The ∆MATR described by the pCO2 drop experiments (from 3X or 4X to 2X) show neither a good agreement with middle to 315 late Eocene, nor with late Eocene to early Oligocene data estimates (Figure 4 a,b and 6 b). The simulations predict no change or a MATR decrease in areas where the Lutetian-Priabonian data points (n=6, triangles, Figure 4) describe increased seasonality (Figure 4 a,b). Priabonian-Rupelian ∆MATR modelled through 4X to 2X and 3X to 2X pCO2 drop (-1.8°C and -1.0°C) are lower than data estimates at similar localities (3.6°C, Figure 7 b). The use of the simulation 4X instead of 3X for the late Eocene stage slightly increases the percentage of agreement for the sign of the ∆MATR, but this result only depends on one high-latitude data-point for which a positive ∆MATR is predicted with a larger pCO2 drop (Table 2, line '%', datapoint 12 in Figure 7.a). In addition, pCO2 drop alone leads to zonal ∆MATR which do not transcribe the spatial heterogeneity visible in data. This misfit is visible through high RMSE scores and the absence of significant correlation between modelled ∆MATR resulting from pCO2 drop (simulations 2X -4X and 2X -3X) and ∆MATR described by proxydata (Table 2). Two data-model agreement patterns are nevertheless to be noted: (1) Regardless of their values (which are 325 higher in data than in our simulations), the northernmost data points are inside or surround the high-latitude seasonality strengthening zone we modelled (Figure 4, data-points 9,5,12 on Figure 7); (2) none of the Southern Hemisphere data localities showing no seasonality change are located within zones MATR increase zones (Figure 4 a,b).

Antarctic ice-sheet and sea-level
The formation of the AIS alone does not result in a better agreement between the modelled and Priabonian-Rupelian ∆MATR estimates. It is even slightly reduced ( Table 2) Table 2).
Geographic changes associated with sea-level drop result in a better agreement with data ∆MATR (Figure 4.d and Figure   7.d). The largest continental fraction changes affect the MATR on a broad geographic scale and allow for a better agreement, 345 even with several distant data points. For example, the retreat of the proto-paratethys reinforces the area of MATR increase at the high northern latitudes which spreads toward North America and creates a better fit with data estimates located on the Pacific Coast. Similarly, coastline changes along the eastern part of Africa and Western Asia cause an increase in seasonality in Anatolia and Central/Western Europe, improving the fit between model and data. These changes in temperature seasonality result in a reduction of 2 to 2.5 months in the duration of the plant growing season (as reconstructed with the 350 formula of Grein et al., 2013; Table S1). Smaller changes in coastlines such as in Florida, Kamchatka Peninsula or along the East Asian coast increase seasonality at a regional scale and improve the data-model fit (data-points 1, 3-5, 14, 17). The simulation still underestimates ∆MATR across the EOT (Figure 7.d; Table 2). Modelled MATR using simulation 2X-ICE-SL for the Rupelian stage are the only ones showing significant correlation with proxy-data ∆MATR (Table 2). Changes are equally correlated to proxy-data Priabonian-Rupelian ∆MATR, independently of the late Eocene simulation (4X or 3X) used 355 as initial stage (Spearman rank test: rho = 0.56 and 0.57, p-value < 0.05, Table 2). This better agreement is also visible in RMSE analysis, with the lowest data-model distance observed for the couple 2X-ICE-SL -3X (Table 2). A smaller datamodel difference is also described for 2X-ICE-SL -3X (as compared to 2X-ICE-SL -4X) by other indicators, notably (1) a 3.2°C mean difference, as compared to a 3.9°C offset for 2X-ICE-SL -4X (2X-ICE-SL -3X) and (2) the percentage of values for which a ∆MATR in the good direction is predicted (23.5% higher with 3X; Table 2). 360

Model climate sensitivity and climate response to EOT forcing
∆MATR across the EOT are better predicted with the lower pCO2 simulation 3X for the late Eocene stage. In addition, with a mean global cooling of 2.7°C between 3X and 2X simulations, surface temperature changes are also in agreement with the 365 mean changes described in marine proxies across the EOT (i.e., difference of 2.9°C between 38-34.2 and 33.7-30Ma, Hutchinson et al., 2021). The best fit with a limited drop of pCO2 reflects the high climate sensitivity of our model (i.e., the https://doi.org/10.5194/cp-2021-27 Preprint. Discussion started: 1 April 2021 c Author(s) 2021. CC BY 4.0 License. average temperature change per doubling of the pCO2 at model equilibrium; PALEOSENS, 2012). This high sensitivity is also highlighted in our experiments of pCO2 halving from 1120 to 560 ppm (2X-4X), which result in a dramatic mean annual global cooling (5.8°C for MAT, 5°C for SST). Such temperature difference is high compared to previous modelling studies 370 which describe a 3 to 4°C surface atmospheric temperature difference under similar forcing, and the 2.9°C cooling found in marine proxies across the EOT (Hutchinson et al., 2021). The pCO2 of 4X and 2X, are more generally used in simulations to represent the transition to an icehouse world (e.g., Baatsen et al., 2020;Goldner et al., 2014;Kennedy-Asser et al., 2019).
Although 4X is likely too high pCO2 for the late Eocene (values rather approximate ~ 800 ppm from the Lutetian to the Priabonian; Foster et al., 2017), the use of this value is justified to better reconstruct high-latitudes temperatures (Huber and 375 Caballero, 2011). A good agreement between warm conditions and Bartonian SST data has also been recently shown by other experiments using the model IPSL-CM5A2 with middle / late Eocene boundary conditions Toumoulin et al., 2020). We thus argue that the use of the 4X simulation is appropriate to study possible variations of pCO2 during the Eocene, but the use of the 3X simulation is better to study the changes between the Priabonian and the Rupelian.

Temperature seasonality changes through the late Eocene 380
The evolution of the different climate parameters responsible for ∆MATR is consistent with several aspects of former studies. First, earlier modelling experiments have described albedo and sea-ice increase resulting in polar amplification of the cooling (e.g. Baatsen et al., 2020;Kennedy-Asser et al., 2019) and a reinforcement of temperature seasonality (Eldrett et al., 2009). The strengthening and expansion of the high Northern latitudes MATR increase zone with pCO2 lowering is a good explanation for the dramatic seasonality increase at high latitudes suggested by some studies (Eldrett et al., 2009;385 Wolfe, 1992;Zanazzi et al., 2015). In addition, changes in the importance and distribution of net precipitation (i.e., precipitationevaporation) resulting from pCO2 lowering agree with former theoretical and modelling work suggesting an intensified hydrological cycle under higher pCO2 (e.g., Carmichael et al., 2016;Hutchinson et al., 2018). This phenomenon results from a greater capacity of the air to retain moisture and more intense atmospheric convection phenomena (Allen and Ingram, 2003;Carmichael et al., 2016;Held and Soden, 2006). Finally, the increase of low cloud cover is consistent with 390 former model studies describing higher fraction of low-level clouds under lower pCO2 (Baatsen et al., 2020;Caballero and Huber, 2013;Zhu et al., 2019). Such a trend could be explained by a higher efficiency of cloud water conversion into precipitation in warmer climates (Zhu et al., 2019). Nevertheless, although a low-level cloud cover increase due to pCO2 drop is consistent with increased air moisture in Western Europe at the EOT (Kocsis et al., 2014), this parameter remains poorly constrained in paleoclimate archives and modelling analysis Sagoo et al., 2013). 395 Despite these agreements, the MATR evolution resulting from pCO2 drop does not clearly match data estimates whether they correspond to both Lutetian-Priabonian or to Priabonian-Rupelian changes. This suggests that the temperature seasonality compiled can only be partly explained by a pCO2 drop. Since zonal ∆MATR patterns are visible with a pCO2 drop of 1PAL (either from 4X to 3X or from 3X to 2X) we hypothesize that they likely occurred before the AIS onset, and that high northern latitudes were among the first to observe seasonality strengthening. However, the major part of Lutetian-Priabonian 400 https://doi.org/10.5194/cp-2021-27 Preprint. Discussion started: 1 April 2021 c Author(s) 2021. CC BY 4.0 License. data is outside this area which prevents from testing this hypothesis (Figure 4 a,b; Supporting Table 1). Similarly, the presence of areas of decreased seasonality due to changes in the hydrological cycle cannot be confirmed because of a lack of data in these areas: although some of the data associated with a decrease in the MATR share the same latitudinal bands, none of them are directly located within a zone of MATR decrease. It might be of interest checking for additional seasonal temperature records to better trace such eventual early trends. This low fit was, to some extent, predictable since the 2X 405 simulation was set ice-free and thus, does not either represent the late Eocene nor the early Oligocene (see Material and Method section). Indeed, small scale glaciations (25-35% modern AIS) have been described already during the late Eocene, before the EOT, associated with a moderate sea-level decrease (Carter et al., 2017;Miller et al., 2020;Scher et al., 2014).
Interestingly, the combination of the three forcing mechanisms also lead to a better agreement of modelled ∆MATR and middle to late Eocene data, especially in coastal areas of Kamchatka, and South China (triangles, Figure 4). Although the 70-410 m sea level decrease from the 2X-ICE-SL simulation is too important for the late Eocene, the better data-model agreement when both AIS and sea-level decrease are considered suggests that small ice-sheet development before the EOT may have played a significant role in driving the middle to late Eocene ∆MATR. Additional sensitivity experiments, with lower sealevel drop, may enable quantifying the sensitivity of coastal localities to before EOT smaller sea-level variations.

Temperature seasonality changes through the EOT 415
The use of two simulations to set up the effect of the AIS onset (with or without a drop-in sea level) is interesting to unravel the mechanisms affecting temperatures. Temperature changes resulting from AIS alone are consistent with former model studies that simulate a highly regional effect on atmospheric temperature, although the changes we observe spread more widely over the Southern Ocean and Australia (see supporting information of Hutchinson et al., 2021 for model maps).
However, the decreasing seasonality zones modelled at high and mid-latitudes of the Southern Hemisphere are associated 420 with an absence of seasonality change in the data, which often display stable vegetation and biomes from the late Eocene to the Rupelian (Hutchinson et al., 2021;Kohn et al., 2015;Nott and Owen, 1992;Pocknall, 1989;Pound and Salzmann, 2017).
This apparent mismatch calls into question the capability of paleobotanical proxies to record temperature seasonality decrease in environments already characterized by low seasonality. Indeed, the decrease in the temperature seasonality is associated with a more pronounced drop in summer temperatures, which is a less limiting factor for flora distribution and 425 thus less constrained in the fossil record than winter temperatures (Huber and Caballero, 2011).
The best representation of the temperature seasonality evolution from Priabonian to Rupelian arises when sea level drop is taken into account. This consequence of the Antarctic glaciation has global repercussions and explains part of the heterogeneity documented in the data (Pound and Salzmann, 2017). Note that our results are very dependent on the paleogeography used in the simulations. Because our Rupelian simulations are constructed from Priabonian geography with 430 sea-level lowering, some of the paleogeographic changes that occurred between both periods are not represented, which may affect our seasonality reconstruction. The gradual northward migration of Australia is not considered; the Neotethys is gradually closed during the early Oligocene but a deep-sea passage remains present in our map, to the north of the Arabian https://doi.org/10.5194/cp-2021-27 Preprint. Discussion started: 1 April 2021 c Author(s) 2021. CC BY 4.0 License.
Plate (Barrier et al., 2018). Failure to take this difference into account could lead to a misrepresentation of ∆MATR in this area. Another source of error may come from fragmented continental areas such as seen in Europe at that time. In these 435 zones, temperature changes recorded through the EOT are heterogeneous with paleovegetation studies suggesting medium (1.8 -2.1°C; Moraweck et al., 2019;Teodoridis et Kvaček, 2015;Tosal et al., 2019) to strong (up to 8.3°C;Tanrattana et al., 2020) MATR increase. A late intensification of seasonality in Central and Eastern Europe, related to changes in geography induced by sea level decline and not by the pCO2 lowering, is consistent with fairly stable vegetation between the middle and late Eocene (Bozukov et al., 2009;Kvaček et al., 2014;Moraweck et al., 2015). This may result from the proximity with the 440 warm Tethys, which may have buffered the EOT cooling, as suggested by stable δ 18 O describing moderate temperature changes in this area (Kocsis et al., 2014). The heterogeneity shown in data might thus result from smaller scale paleogeographic changes through the EOT that are not represented in our simulations.
In addition, differences between our modelling results and data may also be related to the amplitude of the sea-level drop used in our simulation compared to its variability during the Rupelian. The EOT is generally described in two steps: a first 445 event at ~ 33.9 Ma with both a decrease in temperature and sea level (~ 25m) and a second event, the Early Oligocene Glacial Maximum (EOGM), between approximately 33.65 and 33.15 Ma, starting after a large oxygen isotope incursion (often referred to as "Oi-1"), which is characterized by an additional 50 m sea-level decrease (see Hutchinson et al., 2021 for synthesis andterminology, andMiller et al., 2020). The sum of these two steps corresponds to the boundary conditions of our simulation. However, important variations of the East Antarctic Ice-Sheet have been described until the early Miocene 450 (50-60 m sea-level equivalent; Miller et al., 2020). Directly after the EOGM phase, a decrease in ice-volume is visible between 33.15 and 32.8 Ma, before increasing again and remaining stable between 32.8 and 29 Ma (after the "Oi-1a" event; Galeotti et al., 2016). Due to the combined effects of the drop in CO2 and the development of the AIS (and the amplitude of the associated drop in sea level, 70 m), the important changes in seasonality reconstructed here (2X-ICE-SL minus 3X) were probably not in place throughout the Rupelian, but rather for shorter periods during the EOGM, or later between 32.  Most continental paleoclimate studies do not provide the resolution to distinguish these steps. Among the data points we compiled for this study only four sites have enough temporal resolution to be linked to the EOGM phase represented by our 2X-ICE-SL simulation (Bozukov et al., 2009;Eldrett et al., 2009;Hren et al., 2013;Kohn et al., 2015;Tosal et al., 2019).

Perspectives on environmental and biotic crisis
The EOT is associated with major extinction events, of which the best known are the "Grande Coupure" in Europe and the 460 "Mongolian Remodelling" in central Asia (Stehlin, 1909;Meng and McKenna, 1998;see Coxall and Pearson, 2007 for review). Although first recognized as a major biotic turnover for ungulates (Blondel, 2001;Stehlin, 1909), this biological crisis also affected other vertebrates, including rodents, primates, amphibians and squamates, and is also expressed by major changes in vegetation that vary regionally (e.g., Barbolini et al., 2020;Dolezych et al. in press;Eldrett et al., 2009;Marigó et al., 2014;Pound and Salzmann, 2017;Rage, 1986Rage, , 2013Rage and Roček, 2003). These changes have been linked to (1) 465 competitive interactions resulting from the dispersal of Asian taxa to Europe, and (2) EOT climate deterioration and https://doi.org/10.5194/cp-2021-27 Preprint. Discussion started: 1 April 2021 c Author(s) 2021. CC BY 4.0 License. selection processes through resource and/or habitat changes (e.g. Hooker et al., 2004;Kratz and Geisler, 2010;Marigó et al., 2014;Sun et al., 2015;Zhang et al., 2012). The latter mechanism is commonly related to irreversible cooling and/or aridification at the EOT (e.g. Blondel, 2001;Sun et al., 2015;Zhang et al., 2012). Climate cooling may have significantly reduced the habitat of well spread early Eocene tropical (and paratropical) species, characterized by narrow thermal 470 ecological niches (Hren et al., 2009;Huang et al., 2020;Jaramillo et al., 2006;Wing, 1987). Although the distribution of fauna and flora is based on a complex set of parameters, we discuss here how ∆MATR provides an additional interpretative key for understanding biotic turnover at the EOT.
While Asia and North America show comparable temperature changes, our simulations highlight significant differences in the evolution of their MATR, which increase and decrease respectively. Vegetation changes and the Mongolian Remodelling 475 are contemporaneous to AIS growth between 32.8-29 Ma and can be compared with our 2X-ICE-SL simulation (Galeotti et al., 2016;Kraatz and Geisler, 2010;Sun et al., 2015). The MATR strengthening modelled in Asia shows that cooling was particularly strong during winter. In addition to the aridification, this more pronounced winter-cooling may have contributed to the intensity of extinctions in this area (Barbolini et al., 2020;Dupont-Nivet et al., 2007). Conversely, the decrease in MATR in North America provides an explanation for the low impact of the EOT on fauna and vegetation in this area (Coxall 480 and Pearson, 2007;Prothero and Heaton, 1996;Stucky, 1992). Despite a similar decrease of the mean annual temperature, most of the temperature drop is in summer and isn't associated with the onset of cold winters ( Figure S3). We hypothesize that these patterns enabled a greater persistence of existing warm-temperate to paratropical vegetation and associated biota (Pound and Salzmann, 2017). The study of Tardif et al. (in preparation) (Kvaček et al., 2014;Hutchinson et al., 2021). Differences in the evolution of the MATR between North America, Europe and Asia support several studies suggesting different causes for EOT extinctions (Blondel, 2001;Hooker et al., 2004;Meng 495 and McKenna 1998;Sun et al., 2015). Finally, little is known about the EOT in Africa, but the rare data available suggest stable flora and fauna in northern Africa (Hutchinson et al., 2021;Pound and Salzmann, 2017;Rasmussen et al., 1992). The significant increase in seasonality in North Africa seen in our simulations differs from the one observed in other areas, since it is linked to an increase in summer temperatures (barely no change in winter). We propose that unchanged winter https://doi.org/10.5194/cp-2021-27 Preprint. Discussion started: 1 April 2021 c Author(s) 2021. CC BY 4.0 License.
temperatures might explain the apparent lack of biotic response through the EOT in North Africa., though more data are 500 needed on this continent to better document the EOT response.

Conclusion
This study investigates the changes in temperature seasonality during the middle to late Eocene and across the EOT. The combined effects of pCO2 drop, AIS formation and sea-level decrease impact the seasonality in the right direction when compared to the reconstruction of the late Eocene to early Oligocene changes in MATR. The decrease in pCO2 leads to a 505 marked strengthening of seasonality at high northern latitudes, which may have started earlier than the EOT, during the late Eocene. The formation of the AIS and the resulting sea-level drop lead to both an intensification and an extension of increasing seasonality areas. When all three parameters are combined, a statistically significant correlation between data and modeled MATR evolution throughout the EOT is reached. Accounting for sea level changes associated with Antarctic freeze-up appears to be the most important parameter to explain the heterogeneity of ∆MATR across EOT. A seasonality 510 increase is also visible in middle to late Eocene localities, which may reflect earlier sea level changes associated with the precursors of the Antarctic ice sheet. A discrepancy between data and model is present for MATR estimates across the EOT, with less marked seasonality changes in the model. This is mainly due to areas where a decrease in seasonality is predicted by the model while the vegetation proxies show stable vegetation. We hypothesize that this discrepancy can be explained by a low capacity of vegetation to register decreases in summer temperatures compared to winter temperatures. 515 Reconstructing changes in MATR brings additional constraints on the environmental pressure applied on land between the middle Eocene and the early Oligocene. The different mechanisms described here likely explain the heterogeneity in seasonality changes found in data across the greenhouse-icehouse transition, and provide insights into the diversity of continental paleoenvironments. The map of ∆MATR reconstructed here may provide new elements to better understand major extinction events at the EOT. This study only focused on the evolution of temperature seasonality; the variability of 520 other seasonality parameters, including rainfall seasonality, would also be worth investigating to better picture the evolution of terrestrial climates and biodiversity through the EOT. In addition, further work using higher resolution Rupelian paleogeography and regional models would be of great interest to better reconstruct temperature seasonality changes, especially in highly fragmented land areas such as Europe during the Eocene. Modelling Centre, 2021b) for a proper installation and compilation of the environment.

Data availability 540
The key climatological outputs of the simulations are stored in the PANGAEA database (data are currently under review and may be sent to reviewers upon request).

Author contribution
YD, AT and JBL conducted the modelling experiments. AT compiled proxy-data, analyzed the model results, generated figures and tables, and drafted the manuscript. AT, YD, DT and JBL discussed the paleoclimate results. LK contributed to 545 the paleobotanical context, and AL and GDN to the paleogeographic and geological context. All authors have provided critical feedback and contributed to the final manuscript.