Impact of terrestrial biosphere on the atmospheric CO₂ concentration across Termination V

Among the 100 kyr climatic cycles of the Late Pleistocene, Termination V (TV, ∼ 433–404 kyr BP), the fifth last deglaciation, stands out for its minimum in astronomical forcing associated paradoxically with maxima in sea level, Antarctic temperature and atmospheric CO₂ concentration. However, the driving mechanisms explaining TV remain only partially understood. For instance, climate models cannot fully represent the atmospheric CO₂ variation observed in paleoclimate data. Aside from essential oceanic circulation processes, there is increasing evidence that terrestrial biosphere may have played a key role in the global carbon cycle. This study proposes a three-step integrated approach, combining regional and global vegetation records with modelling results, to unveil the evolution of terrestrial biosphere and its contribution to the carbon cycle during TV. First, we provide a new high-resolution (∼ 700 years) deep-sea pollen record from the Gulf of Cádiz (site U1386, 36^∘49.680^′ N; 7^∘45.320^′ W) for TV, which shows a moderate expansion of the Mediterranean forest. We then construct the first global forest pollen database for this period. Our compilation features distinct evolutions for different types of forest, highlighting a strong development of temperate and boreal forest which might have delayed the atmospheric CO₂ increase during TV. Finally, the direct comparison of global simulated forests (iLOVECLIM model) to our pollen database reveals overall consistent temperate and boreal forest evolutions despite model biases, thereby supporting the hypothesis of a significant CO₂ sequestration by middle and high-latitude forests of the Northern Hemisphere shortly after the onset of TV.

1 Introduction

Among the five major short-term carbon reservoirs of the climate system (ocean, atmosphere, terrestrial biosphere, surface sediments, permafrost), terrestrial biosphere (including vegetation and soils) is regarded as an essential component of current anthropogenic climate change mitigation strategies (Harris et al., 2021; Shukla et al., 2019). Yet, natural land sinks are also marked by large internal climate variability resulting in uncertainties in the global carbon budget (Loughran et al., 2021). In this respect, the study of past climate changes and associated terrestrial biosphere responses remains crucial to reduce uncertainty related to Earth's sensitivity to climate forcings (Overpeck et al., 2003; Masson-Delmotte et al., 2021). Paleoenvironmental data also allow us to investigate terrestrial biosphere feedbacks on climate at multi-millennial timescales (Overpeck et al., 2003). However, while extensive studies cover the last glacial cycle (Hoogakker et al., 2016; Prentice et al., 2011) and the Holocene (Bartlein et al., 2010), there are only few global descriptions of terrestrial biosphere for older periods of the Quaternary which are marked by different climatic boundary conditions (e.g. Kleinen et al., 2014).

Termination V (TV, ∼ 433–404 kyr BP), the fifth last deglaciation, is defined as the period between the highest and the lowest values of the benthic foraminifera oxygen isotope ratio (δ¹⁸O_b) following seminal work by Broecker and van Donk (1970) and by Broecker (1984). TV begins during glacial Marine Isotopic Stage (MIS) 12 and extends to MIS11, thus covering the longest deglaciation period (∼ 29 kyr) in the last 650 kyr (Sarnthein and Tiedemann, 1990) while orbital precession is at the lowest. Indeed, Berger and Wefer (2003) refer to the “Stage-11 Paradox” in the following terms: “The amplitude of climate variation is at maximum in the vicinity of Stage 11, at the very time when astronomical forcing is at minimum.” Such a conflicting observation suggests that the deglaciation process embodied complex multi-scale mechanisms and emphasises the necessity of addressing the problem within the framework of the entire climate system (cryosphere, ocean, atmosphere, biosphere) together with the Milankovitch theory (Berger and Wefer, 2003). It also highlights the relevance of MIS 11 as an astronomical analogue to the present Holocene interglacial (MIS1) featuring low eccentricity (Loutre and Berger, 2003; Berger and Wefer, 2003; Candy et al., 2014).

Parallel to the Stage-11 Paradox, TV marks a change within the carbon cycle of the 100 kyr climate cycles. As of MIS11, younger interglacials are characterised by higher atmospheric CO₂ concentrations (by about 40 ppm; Lüthi et al., 2008) and temperatures (Jouzel et al., 2007) compared to older ones. This climatic transition is known as the Mid-Brunhes Event (MBE, ∼ 430 kyr BP). First model simulations by Köhler and Fischer (2006) have shown that low surface temperature in the Southern Ocean together with a reduction in the Atlantic Meridional Overturning Circulation (AMOC) accounted for low atmospheric CO₂ concentrations during pre-MBE interglacials. Simulation experiments also indicate that the pre-MBE period is marked by intense Antarctic Bottom Water (AABW) formation and significant ventilation of the Southern Ocean contributing to carbon storage in the ocean (Yin, 2013). However, recent simulations including a representation of the carbon cycle reveal that the Southern Ocean contribution is too weak to fully explain the carbon uptake during pre-MBE interglacials (Bouttes et al., 2018, 2020). This result suggests a misrepresentation of the magnitude of oceanic processes and/or a lack of crucial mechanisms involved in the carbon cycle such as permafrost or biosphere feedbacks (Bouttes et al., 2018; Barth et al., 2018). Antarctic ice core measurements also indicate that an increase in atmospheric δ¹⁸O during the MBE cannot be fully explained by ocean ventilation (Landais et al., 2010), calling for other explaining mechanisms such as a change in biosphere productivity. Yet, a recent study by Brandon et al. (2020) demonstrates that the ocean primary productivity is much weaker than the terrestrial one during TV. Additionally, the observed strong increase in carbonate production is not associated with an increase in atmospheric CO₂ as expected but rather with a dampened signal suggesting a compensating process (Brandon et al., 2020).

Terrestrial biosphere interacts with the global climate through the photosynthesis–respiration feedback involving CO₂ exchanges. So far, the hypothesis of terrestrial biosphere contribution to the varying atmospheric CO₂ across the MBE through changes in productivity or geographical distribution remains open. Brandon et al. (2020) suggest that terrestrial biosphere productivity (a major carbon flux component) could have delayed the atmospheric CO₂ increase through photosynthesis-driven carbon fixation during the 425–415 kyr BP period. Several modelling studies already address the question of climate–carbon–vegetation feedbacks over MIS11 using either time slice (Kleinen et al., 2014; Rachmayani et al., 2016) or transient experiments (Kleinen et al., 2016; Ganopolski and Brovkin, 2017). However, these studies either do not fully cover TV or use simple models leaving room for a global model assessment of climate–vegetation interactions during the whole deglaciation period. Additional regional and global observations of paleo-biosphere evolution (Brandon et al., 2020) together with a better representation of biosphere in the coupled climate–vegetation models (Bouttes et al., 2018) are required to identify the impact of the terrestrial biosphere on atmospheric CO₂ concentration across TV. The present study is a first step to fill in this gap, laying the groundwork for future evaluation of the full carbon cycle. It follows a three-step data–model integrated approach from a regional pollen record analysis to a global assessment of forest evolutions.

First, we analyse a new pollen record (U1386) from the Iberian margin to clarify the vegetation evolution during TV in southern Iberia. This region is a key location in terms of climate because it is affected by low- and high-latitude processes.

We then establish the first global pollen database for TV in order to understand terrestrial biosphere evolution at a global scale. The pollen records (including U1386) are analysed in a bioclimatic framework, therefore allowing us to capture the spatial variations in terrestrial biosphere evolution across TV.

Finally, we run transient simulations over TV to evaluate the quality of terrestrial biosphere representation by comparison to the database and to improve our understanding of the role of terrestrial biosphere in the carbon cycle.

2 Material and methods

2.1 Pollen analysis of IODP site U1386

2.1.1 Regional setting

The Integrated Ocean Drilling Program (IODP) site U1386 (36^∘49.685^′ N, 7^∘45.321^′ W; 561 m water depth) was drilled in 2011 during the Mediterranean Outflow Expedition no 339. The site is located in the Faro Drift (Fig. 1), on the upper slope of the Iberian margin (Gulf of Cádiz). As of today, no full record of TV exists for southwestern Iberia. While site MD01-2447 only describes temperate forests in northern Iberia (Desprat et al., 2005), southern Iberian Mediterranean forest evolution remains partially unknown because of a sedimentary hiatus at site U1385 (Oliveira et al., 2016). The U1386 sequence covers a quasi-continuous period from the late Miocene to Holocene (Stow et al., 2013), crossing hemipelagite, turbidite and contourite depositional systems. The core section of interest for TV (Unit IA) is embedded in a unique contourite structure shaped by the historical variations in the Mediterranean Outflow Water (MOW) originating from the Strait of Gibraltar. The local lithology is composed of nannofossil muds, calcareous silty muds and silty bioclastic sands with a dominance of silty muds (Stow et al., 2013; Moal-Darrigade et al., 2022a) and is overall homogeneous across TV. The size of pollen grains (10–100 µm) falls within the silt fraction and a recent study has shown that fine sediment fraction, including pollen, found at site U1386 essentially originates from the Guadalquivir and the Guadiana rivers (Alonso et al., 2016; Moal-Darrigade et al., 2022b). Moreover, modern pollen samples from the southern Iberian margin deep-sea floor have proven to accurately represent the vegetation of the adjacent continent (Naughton et al., 2007; Morales-Molino et al., 2020). Therefore, the pollen record from site U1386 reflects the regional vegetation evolution of southwestern Iberia. This vegetation presently belongs to the thermo-Mediterranean and meso-Mediterranean bioclimatic stages across the Guadalquivir and the Guadiana catchments (Rivas-Martinez, 1987). Deciduous oaks are found at middle elevation while evergreen oaks, olive trees, Pistacia, Phillyrea and rockroses (Cistus) dominate the lower elevations (Castro et al., 1997). The present-day climate of the southern Iberian Peninsula is directly controlled by the North Atlantic pressure systems dominated by cyclonic cells modulated by the North Atlantic Oscillation (NAO) in winter and by the Azores subtropical high during summer (Lionello et al., 2006). These strong seasonal atmospheric patterns result in mild winters and hot and dry summers (Peinado-Lorca and Martínez-Parras, 1987).

Figure 1 Study area featuring main ocean currents (arrows), vegetation (shading) and sediment cores (dots). Atlantic currents are shown in blue (North Atlantic Current: NAC; Atlantic Current: AC; Portuguese Current: PC) and the Mediterranean Outflow Water (MOW) in red. Colour shading indicates the modern vegetation distribution in western and central Europe. The present pollen record originates from site U1386 (purple dot). The Guadalquivir catchment is delimited by a dotted black line. Courtesy of Vincent Hanquiez.

2.1.2 Experimental setup

A total of 42 pollen samples from site U1386, spanning about 7 m between 145 and 138 corrected metre composite depth (cmcd), were analysed. The composite core was sampled approximately every 12 cm between 143.39 and 139.81 cmcd and more loosely at the edges of the sequence (∼ 50 cm). Sample volumes ranged from 2.5 to 5 mL with an average value of 4.1 mL. Each sample was prepared for pollen analysis following the standard palynological procedure for marine samples in use at UMR EPOC, University of Bordeaux (https://ludovicdevaux.wixsite.com/ephe-paleoclimat/about-the-lab, last access: 23 June 2022). First, a coarse sieving (150 µm mesh) allowed the separation of the very coarse fraction of the sediment from the finer elements including pollen, spores and dinoflagellates. The resulting fraction was treated with gradually increasing concentrations of HCl and HF to remove carbonates and silicates respectively and then sieved through a finer mesh (10 µm). Known quantities of Lycopodium spores were added to the residue to allow reconstruction of absolute pollen concentrations. Finally, the sample was mounted on a glass slide in glycerol, a mobile medium enabling pollen grain rotation during observation. Pollen identification was performed at 400× and 1000× magnifications with a Zeiss AXIO Imager A1 microscope.

In order to obtain a reliable representation of the sample compositions at site U1386, we observe both a quantity and a diversity criterion following the methodology by Janssen (1981) and Mc Andrews and King (1976). For each sample, at least 100 terrestrial pollen excluding Pinus, aquatics and spores were counted (except for one sample with only 93 pollen counts). Thus, the main pollen sum including Pinus ranges from 128 to 770. Additionally, the majority of samples is composed of at least 20 different pollen morphotypes to accurately depict the floristic diversity of the source vegetation (Mc Andrews and King, 1976). Rull (1987) shows that a pollen sum of 200 grains, including Pinus, is sufficient to produce reliable estimates. Indeed, for higher values there are no significant variations in the confidence interval width. It is further argued that the pollen percentages derived from 100 grain counts (excluding the dominant taxon – Pinus here) yield a reliable picture of the considered sample, within the 0.95 confidence limits (Maher, 1981). Pollen concentrations range between ∼ 2500 and 52 000 grains cm⁻³ (Fig. S2 in the Supplement), which are typical values for deep-sea sedimentary sequences (1–2 orders of magnitude lower than for lacustrine sequences) and remain above the reliability threshold (∼ 2500 grains cm⁻³) for pollen assemblages containing deteriorated grains (Hall, 1981). We observe both strong pollen concentration and MOW (reflected by a high $\mathrm{Zr}/\mathrm{Al}$ ratio and high sortable silt values) at the end of TV but low pollen concentration and high MOW during MIS11, suggesting that pollen concentration is relatively independent from MOW strength (Fig. S2). Counting results are expressed as pollen percentages calculated against the main sum (which discards Pinus) for terrestrial taxa, against the main sum plus Pinus for the Pinus taxon, and against the total sum (pollen + spores + indeterminables + unknowns) for the aquatic pollen and spores. Such a distinction avoids flattening of the pollen signal by the overrepresented Pinus in marine sediments (Heusser and Balsam, 1977).

Table 1 Taxonomic composition of the six ecological groups.

Each taxon was assigned to an ecological group representative of specific climate conditions. In Table 1, we define six ecological groups which are relevant for the Iberian vegetation following previous studies on European vegetation (Polunin and Walters, 1985; Suc, 1984). The Mediterranean taxa group (MTss) is defined by typical sclerophyllous Mediterranean taxa. This group indicates a seasonal climate with cool and wet winters and hot and dry summers (Polunin and Walters, 1985; Van Campo, 1984). The Mediterranean forest (MF) is a broader group including the Mediterranean taxa group together with all temperate and moisture-loving tree (Populus and Salix are included here following Polunin and Walters, 1985 and Castro et al., 1997) and shrub taxa excluding Pinus, Cedrus, Hippophae, Helianthemum and Cupressaceae. The pioneer group encompasses all plant species able to rapidly colonise an ecosystem with moderate increase in temperature and precipitation and poor-nutrient soils. The semi-desert (SD) group is composed of shrubs developing in arid conditions (Amaranthaceae, Ephedra distachya-type, Ephedra fragilis-type and Artemisia), following Van Campo (1984), while Ericaceae and Calluna, typical of annual wet climate, form the heathland group (Polunin and Walters, 1985). Finally, taxa which can be found in multiple climates fall in the ubiquist group.

2.1.3 From pollen assemblage to vegetation and climate reconstruction

In order to associate the depth of our composite core with a timescale, we use an improved version of the age model by Kaboth et al. (2017). The age correspondence is based on direct visual correlation between the δ¹⁸O_b record at site U1386 and the global mean benthic isotope stack LR04 (Lisiecki and Raymo, 2005), for which the chronology is well established. We construct our age model by linear interpolation based on the calibration points (Table S1 in the Supplement) of Kaboth et al. (2017) and taking into account the splice correction of Moal-Darrigade et al. (2022a) (+0.17 m on specific cmcd points). The sequence covers a total of about 28 kyr assuming a constant sedimentation rate of 0.2 m kyr⁻¹ (Kaboth et al., 2017) and features a highest resolution of approximately 500 years (Fig. S1). We estimate the uncertainty on the chronology to be ∼ 1.2 kyr by computing the quadratic sum of the different age error sources (pollen record resolution, δ¹⁸O record resolution) following the methodology by Govin et al. (2012).

The onset of Termination V, as other Middle and Late Pleistocene terminations, is defined by a drop in the δ¹⁸O_b below the 3.5‰ threshold (McManus et al., 1999). Averaging over the highest and lowest δ¹⁸O_b (Skinner and Shackleton, 2005; Shackleton et al., 2003) allows positioning the transition 12/11 (i.e. the end of MIS12 glacial period and the beginning of MIS11 interglacial period) at ∼ 428 kyr BP, which falls in the range of Termination V midpoint values found by Raymo (1997). The same method applied to local δ¹⁸O_b minima and maxima defines the onset of the substages MIS11d at ∼ 422 kyr BP and 11c at ∼ 418 kyr BP.

For the cluster analysis, the pollen diagram computed with Psimpoll software (Bennett, 2000) was divided into a subset of pollen zones defined by the variation in the percentages of at least two pollen taxa with different ecological affinities (Birks and Birks, 1980). This visual interpretation is supported by a hierarchical clustering analysis based on the minimisation of the Euclidean distance and constrained by the number of sample depths (chclust in R environment, Juggins, 2009). The resulting dendrogram defines six pollen zones across TV.

2.2 Global pollen database

2.2.1 Pollen record selection

We conducted an extensive literature review to compile the available marine and lacustrine pollen records covering TV. Only records spanning at least the 423–405 kyr BP period were selected. At a global scale, a total of 19 pollen records (including the present U1386 record) with variable resolution were made available by the scientific community (Pangaea database, European pollen database, Neotoma database). However, we note the existence of a few more pollen records covering TV that are not included in this work due to a lack of data accessibility (Prokopenko et al., 2002; Tarasov et al., 2011; Litt et al., 2014; Cheng et al., 2018) or the use of different metrics such as pollen concentration instead of pollen percentages (de Vernal and Hillaire-Marcel, 2008). The deep-sea pollen sequence located off Greenland (de Vernal and Hillaire-Marcel, 2008) only provides pollen concentration values that cannot be taken into account for the construction of the composite pollen percentage records. As each pollen record featured different pollen assemblages according to their location, we extracted (or defined and computed when it was lacking) the arboreal pollen percentage from each pollen assemblage (“Forest proxy” column, Table 2) in order to build the most representative metric for forest evolution globally. In the specific case of site GIK16415-2 record, we use fern spores as a forest proxy, for their percentage accurately reflects the forest development (Dupont and Agwu, 1992). As marine and large lacustrine sites collect sediments from surrounding watersheds, each pollen record depicts the regional vegetation setting (Bradshaw and Webb III, 1985; Morales-Molino et al., 2020). Table 2 summarises essential information for each pollen record.

Table 2 Database summary table. Pollen record altitudes (Alt) are given in metres below sea level (b.s.l.) or above sea level (a.s.l.). Pollen selected to build the arboreal pollen percentage are described in the Forest proxy column. Each pollen record is allocated to a group (Grp column): ES-BA (euro-Siberian–boreo-American), Subtropical-N (Subtropical Northern Hemisphere), T (Tropical) or Subtropical-S (Subtropical Southern Hemisphere). The penultimate column (Res) features the mean temporal resolution (kyr) of the sequence. The location of pollen records is displayed in Fig. 5a.

2.2.2 Data analysis

Each pollen record is resampled with a 2 kyr time step in order to build consistent composites. The time step value was chosen in order to optimise the database resolution and minimise oversampling of coarse-resolution records (the majority of pollen records features a mean resolution below 2 kyr, Table 2). Pollen records are classified into four forest groups defined according to their main bioclimatic specificity. The euro-Siberian–boreo-American group (ES-BA), as defined by Braun-Blanquet (1930), covers boreal and temperate forests. Therefore, both Tibetan Plateau pollen records (Heqing and ZB13-C2) are counted in the ES-BA group, the high altitude compensating for their relatively low latitudes. Although Tenaghi Philippon and Lake Ohrid records are on the geographical boundary with the Mediterranean region, they are also classified as the ES-BA group for they reflect a temperate forest cover. The Northern Hemisphere Subtropical group (Subtropical-N) encompasses warm temperate forests of the Northern Hemisphere and therefore extends beyond the geographical boundaries of the Mediterranean region (Mediterranean sensu lato). The Tropical group (T) refers to tropical and Afromontane forests. Finally, the Southern Hemisphere Subtropical group (Subtropical-S) only includes the MD96-2048 record, which is at the limit between tropical and subtropical forests (Table 2). Forest pollen composites are then obtained by averaging the forest pollen percentage for each group. It should be noted that all groups (i) are made of very few records and (ii) do not contain the same number of records. As a consequence, it is difficult to compare the composites between them and to estimate a statistical robustness of their signal. However, a sensitivity analysis (not shown) performed on the composites for varying resampling time step values (between 1 and 4 kyr) revealed no significant change. Besides, as each record has its own chronology (and associated uncertainty), the present analysis only allows us to compare long-term trends in the forest pollen signal.

2.3 Model and experimental setup

2.3.1 Model description

We use the three-dimensional Earth system model of intermediate complexity iLOVECLIM, a newer version of the LOVECLIM version 1.2 (Goosse et al., 2010) resolving the carbon cycle (Bouttes et al., 2015). iLOVECLIM includes an ocean module (3^∘ × 3^∘ horizontal grid, 20 vertical levels), a simplified atmospheric module (∼ 5.6^∘ × 5.6^∘ horizontal grid, on three vertical layers), a biogeochemical module (Bouttes et al., 2015, 2020) and a terrestrial biosphere module (VECODE; Brovkin et al., 1997). For each terrestrial grid point, the vegetation module computes two main plant functional types (PFTs), herbaceous (g) and trees (t), and a desert fraction (d) such that their sum equals 1 (i.e. $g+t+d=\mathrm{1}$). For example, an ice sheet grid point will feature null tree and herbaceous fractions and a desert fraction of 1 (Roche et al., 2007). We estimate the terrestrial biosphere carbon stock as the sum of above- and below-ground carbon (e.g. green biomass + structural biomass + slow soil organic matter (SOM) + fast organic matter SOM) derived from the net productivity computed for each continental grid following Brovkin et al. (2002).

2.3.2 Model simulations

Using iLOVECLIM, we run a long-time experiment covering TV (Degla-V, 436–404 kyr). While the present study focuses on TV, we also run a simulation of TI (Degla-I, 21–0 kyr) to check consistency with previous studies of the last deglaciation (Fig. S6, Roche et al., 2011; Shakun et al., 2012). In order to obtain stable initial conditions, a spin-up (stationary) simulation is run over 5 kyr for each experiment (∼ 3 kyr are required to reach an equilibrium in the deep ocean) by fixing the orbital parameters (Berger, 1978), the greenhouse gases (CO₂ and CH₄) and an ice sheet model output (Ganopolski and Calov, 2011) at the time of maximum ice sheet (436 kyr BP for Degla-V and 21 kyr BP for Degla-I). The two transient simulations are then run with the forcings (orbital parameters, greenhouse gases and ice sheet), now varying with time.

3 Results

3.1 Pollen record and inferred southwestern Iberia vegetation and climate

We identify five major pollen zones corresponding to major shifts in vegetation cover and composition (Fig. 2). A summary of the following vegetation description for each pollen zone is provided in Table 3.

Figure 2 Pollen percentage diagram of selected morphotypes and ecological groups (described in Table 1) at site U1386 versus depth (cmcd). Lithology (full grey: muddy clay; dotted black: clayey silt; from Moal-Darrigade et al., 2022a) and sporo-pollen (S-P) concentrations (grains cm³⁻³) are displayed on both sides of the diagram (left and right respectively). Alternating grey and white shading with corresponding numbers indicates pollen zones, and the black line sets the limit between sub-zones U1386-2a and U1386-2b. The dendrogram showing the results from the hierarchical clustering analysis is represented on the right. The purple dashed line indicates the onset of the SINES regional terrestrial interglacial. The same colours are used for the ecological groups and the main morphotype components thereof.

Table 3 Description and interpretation of the TV pollen record at site U1386 by pollen zones.

The first major zone U1386-1 (∼ 145.1–143.7 cmcd) is dominated by ubiquists (∼ 30 %), SD (∼ 40 %) and pioneer taxa (∼ 20 %), while the MF remains low (<10 %). The important presence of SD such as Amaranthaceae, Ephedra and Artemisia suggests a dry and cold climate consistent with the MIS12 glacial period. The large proportion of pioneer taxa (mainly Cupressaceae) and low MF further confirms that climate conditions were too harsh for large temperate trees to grow, leaving a shrub-dominated land cover. The end of zone U1386-1 is characterised by a decrease in SD and pioneer percentages, while the ubiquist proportion increases drastically mainly driven by Taraxacum and Poaceae expansion. Interestingly, this rapid shift occurs immediately after an extremely high pollen concentration peak (> 50 000 pollen cm⁻³, Fig. S2 in the Supplement) highlighting a significant increase in the ubiquist pollen number and a real transition from an SD-dominated assemblage to a ubiquist-dominated one. Such an important pollen concentration change can be caused by a wide range of factors including changes in pollen productivity and dispersability, in the source area and in the distance to the sampling site but also in the rate of sediment deposition associated with sea level rise.

The second major zone U1386-2 (∼ 143.7–142.6 cmcd) is characterised by a dominance of ubiquists (∼ 70 %) mainly represented by Taraxacum and Poaceae. The SD and pioneer groups are much weaker than in zone U1386-1 (∼ 15 % and <10 % respectively), while the heathland group grows up to 10 % across the zone providing evidence for increasing humidity. The beginning of a continuous Isoetes spore record – albeit its weak intensity – also describes a moist environment. However, the dominance of Taraxacum indicates a limited moisture availability. In modern pollen assemblages, the highest Taraxacum percentages are found in Greece and Morocco, where mean annual temperatures span from 12 to 17 ^∘C and precipitation from 0 to 250 mm (Leroy, 1997). We divide zone U1386-2 into two sub-zones (U1386-2a and U1386-2b), considering the emergence of MF taxa (mainly Olea and Pistacia) at around 143 cmcd. Sub-zone U1386-2a corresponds to a relatively humid climate, but it is still too mild to bear forest development. Conversely, sub-zone U1386-2b is both humid and warm enough to allow a first development of the Mediterranean forest. The overrun of the 20 % threshold in MF at the limit between U1386-2b and U1386-3 determines the onset of the regional southern Iberian terrestrial interglacial (Sánchez Goñi et al., 1999; Shackleton et al., 2003) – the SINES forest phase previously defined by Oliveira et al. (2016). This settled and developing Mediterranean forest indicates atmospheric warmth and moisture availability but also a seasonal climate (Sánchez Goñi et al., 2016).

The MF (mainly Quercus deciduous and Pistacia) and the ubiquist groups are still present throughout Zone U1386-3 (∼ 142.6–141.4 cmcd), oscillating around 20 % and 50 % respectively, indicating continuous development of forested area. Zone U1386-3 is also characterised by a higher development of heathlands (10 % to 20 %) than in zone U1386-2 mainly driven by increased Ericaceae proportions. Moreover, the significant rise of Isoetes spores (∼ 10 %) starting around 142 cmcd is a strong signal for deglaciation. In fact, Isoetes grows in marshlands, which gradually expand during an increase in precipitation associated with the transition between glacial and interglacial climates (Sánchez Goñi et al., 1999). The end of Zone U1386-3 is marked by a strong drop in pollen concentration, which could be caused by the abrupt and concomitant sea level rise (Fig. S2), which increases the distance between site U1386 and land, therefore providing further evidence for ice sheet collapse during the transition between Zone U1386-3 and U1386-4. The beginning of Zone U1386-4 (∼ 141.4–140.3 cmcd) is defined by a strong decrease in the heathland group (down to less than 10 %) mainly owing to a decline in Ericaceae. Yet, Zone U1386-4 features one of the highest Ericaceae percentages of the whole sequence.

The onset of U1386-4 is also marked by an increase in the ubiquist group (up to 70 %), while the SD group remains at relatively low levels (∼ 10 %), changes which could be attributed to a wetter or warmer climate. However, the MF also decreases and remains low (mostly below 20 %) throughout the period while Isoetes spores grow over 10 %. These variations describe a decreasing atmospheric temperature and still wet conditions. Moreover, the very stable pollen concentration (∼ 4000 pollen cm⁻³, Fig. S2) reflects an absence of major changes in deposition factors, which are mainly related to sea level and oceanic currents. Therefore, we hypothesise that atmospheric temperature decrease is the major environmental change influencing vegetation in Zone U1386-4.

The transition to Zone U1386-5 (140.3–138 cmcd) is characterised by a surge of the MF group (up to 40 % mainly owing to increasing deciduous Quercus) and a decrease in the ubiquist group (down to 35 %), while the SD group and Isoetes spores remain relatively constant and MTss pollen percentages increase. As pollen concentrations do not vary much over this period, we can infer a strong expansion of the Mediterranean forest. In Zone U1386-5, the MF and ubiquist group percentages follow opposed variations while the other groups remain relatively stable, depicting variations in the percentage of forested and open vegetation land cover. A second MF peak is reached around 139 cmcd, suggesting a stable development of the Mediterranean forest that indicates the warmest climate of the considered interval and wet winter conditions.

3.2 Southern Iberian terrestrial biosphere evolution across TV

We now interpret the U1386 pollen record in a temporal framework by comparing reconstructed vegetation to other climate indices (Fig. 3). Given that all proxies originate from the same sediment core, a direct comparison is possible without chronology issues.

Figure 3 U1386-retrieved environmental indicators. From bottom to top: pollen percentages (this study), warm and cold foraminifera abundance (Moal-Darrigade et al., 2022a), δ¹⁸O_b (from U1386 in blue line (Kaboth et al., 2017) and global in black dashed line (Lisiecki and Raymo, 2005)) and astronomical parameters (source: http://vo.imcce.fr/insola/earth/online/earth/online/index.php, last access: 23 June 2022). The dashed lines represent the 95 % confidence intervals for MF, heathland and SD (computed with the “exactci” package, R). TV period defined by δ¹⁸O_b variations is represented by the blue box, and the SINES regional terrestrial interglacial is indicated by the black hatched box.

The first important result is that the SINES regional terrestrial interglacial (∼ 424 kyr BP, > 20 % MF) begins ∼ 4000 years after the marine isotopic interglacial (∼ 428 kyr BP, δ¹⁸O_b). The onset of MIS11e is characterised by a minimum of precession and increasing insolation at 65 ^∘N responsible for a significant decrease in the ice volume visible in the drop in δ¹⁸O_b. This important melting period is also supported by records from nearby cores (Fig. S3). The increase in the $\mathrm{Si}/\mathrm{Sr}$ values from site U1308 in the central North Atlantic (Fig. 1) indicates the presence of ice-rafted debris (IRD) carried by icebergs originating from the North American ice caps, while the increase in C_37:4 values from an Iberian margin core (MD03-2699) located slightly north of site U1386 (39^∘ N, Fig. 1) depicts freshwater pulses in the southwestern Iberian margin at the end of MIS 12 (Fig. S3). Such a cold water input in the northern North Atlantic from the end of MIS12 until the beginning of MIS11e results in a weakening of AMOC highlighted by low δ¹³C values, a decrease in sea surface temperature supported by the high abundance of Neogloboquadrina pachyderma (s), and cold and dry conditions in southwestern Iberia revealed by the SD taxa expansion. The onset of MIS11e coincides with a severe drop in SD taxa and a rise in ubiquists. These opposite evolutions suggest an increase in atmospheric humidity. However, the MF only develops several millennia later, towards the end of MIS11e. This first moderate forest development is driven by the growth in sclerophyllous Mediterranean taxa (MTss) particularly adapted to the low-precession-induced high seasonality. A decoupling between the evolution of ice volume and the Iberian vegetation is observed during TV similarly to the MIS 5e and the Eemian terrestrial interglacial (Sánchez Goñi et al., 2013). The present U1386 high-resolution pollen record allows repositioning the beginning of the regional terrestrial interglacial at ∼ 424 kyr BP, about 4 kyr after the onset of the marine interglacial (MIS11).

TV further stands out owing to a two-step forest development, absent from more recent terminations (see Fig. 5 in Sánchez Goñi et al., 2018). The first increase in Mediterranean forest (∼ 424 kyr BP) is relatively limited (∼ 20 %) compared to the second (∼ 40 %), which is more representative of a true forest expansion in the southern Iberian Peninsula. The reduced Mediterranean forest develops across MIS11d synchronously with heavier δ¹⁸O_b, suggesting a temporary ice sheet growth. The increasing obliquity causes enhanced insolation seasonality with higher precipitation in winter while high precession results in increased summer precipitation over the Iberian region (Bosmans et al., 2015), accounting for the growth of heathland taxa which require year-round precipitation. The second forest development starting around 414 kyr BP occurs with slightly different orbital parameter settings: high obliquity but decreasing precession such that the highest MF expansion coincides with the insolation maximum. A combination of low precession and high obliquity is known to be even more favourable to winter precipitation (Bosmans et al., 2015). Moreover, the high abundance of warm-loving foraminifera indicates a warmer ocean than at the beginning of the SINES interglacial. This warming is supported by the reconstructed atmospheric temperature featuring a difference of about 2 ^∘C between 424 and 414 kyr BP (Fig. S4). Such environmental conditions are consistent with a larger growth of the Mediterranean forest. Finally, we observe an abnormal contraction of the MF at the beginning of MIS11d (∼ 418 kyr BP), which can be explained by the combination of maximum obliquity and precession accounting for low seasonality. This relatively mild and humid climate is however adapted for the development of Ericaceae, which are responsible for the peaks in heathland taxa. Overall, even the highest Mediterranean forest pollen percentages of MIS11 are relatively low in regard to recent interglacials such as MIS1 and MIS5e (see Fig. 5 in Sánchez Goñi et al., 2018), suggesting a moderate development of the terrestrial biosphere in southern Iberia during TV.

Figure 4 Comparison of Mediterranean forest (%) from sites U1386 (our study) and U1385 (by Oliveira et al., 2016), see map.

Our new high-resolution pollen record (U1386) sheds light on the southern Iberian terrestrial vegetation across TV, a poorly documented period owing to a sedimentary hiatus in the previous IODP site U1385 pollen record (Oliveira et al., 2016). A comparison of U1386 and U1385 records in Fig. 4 reveals a more complex Mediterranean forest evolution than previously thought. While the overall trend of both series is consistent, we observe a slight shift between U1385 and U1386 over the last part of the sequence (U1385 leading U1386 by about 2 kyr, with a maximum at 415 and 413 kyr respectively (Fig. 4)), which we attribute to differing age models. The U1386 record allows us to define a more precise onset of the Mediterranean forest at ∼ 424 kyr BP. As previously mentioned, the U1386 record unveils a two-step forest development absent from the U1385 record, with a first moderate expansion between 424 and 418 kyr BP (onset of the SINES Mediterranean forest phase) and a second larger expansion starting at 414 kyr BP (Fig. 4). Finally, both U1386 and U1385 feature a maximum Mediterranean forest percentage of ∼ 50 %, therefore agreeing on a moderate Mediterranean forest development during the SINES forest phase.

3.3 Global terrestrial biosphere evolution across TV

3.3.1 Global pollen composite records

In this section, we take a step back from our regional analysis to investigate the evolution of the terrestrial biosphere at a global scale owing to the new pollen database. Altogether, forest ecosystems represent the largest biomass reservoir during the preindustrial Holocene period (Prentice et al., 2011) and the Last Glacial Maximum (Crowley, 1995; Prentice et al., 2011). Additionally, the carbon uptake by modern forests accounts for the large majority of the terrestrial carbon sink (Pan et al., 2011). In this respect, forests are reliable proxies for the terrestrial biosphere. Therefore, we postulate that the evolution of the pollen-based forest percentage is representative of terrestrial biosphere evolution (Brandon et al., 2020). Figure 5a features the pollen record database locations and Fig. 5b the corresponding forest pollen percentages across TV. At first sight one notices a large percentage difference between pollen records (see also Fig. S5), some records averaging around 10 % (e.g. ODP 108-658) and others rising over 90 % (e.g. Tenaghi Philippon). This reveals a strong diversity of environments between the sampled regions. Building composites allows comparing records in presumably similar global climatic conditions to extract information concerning the long-term evolution of forests in a specific location.

Figure 5 (a) Forest pollen record locations (dots, numbering refers to Table 2) classified into four bioclimatic groups defined in the text (red: Subtropical-S; purple: Tropical; yellow: Subtropical-N; blue: euro-Siberian–boreo-American, ES-BA). Stars refer to the pollen records covering TV and unavailable at the time of the study: numbers 20 to 24 correspond to sites ODP646 (de Vernal and Hillaire-Marcel, 2008), GeoB16602 (Cheng et al., 2018), Lake Van (Litt et al., 2014), Lake Baikal (Prokopenko et al., 2002) and Lake Biwa (Tarasov et al., 2011) respectively. Colour shading represents the modern biome distribution modelled by multinomial logistic regression (from Levavasseur et al., 2012) on the BIOME 6000 database (Prentice et al., 1996; Pickett et al., 2004) and the black line delineates the ES-BA megabiome southern boundary (courtesy of Vincent Hanquiez). (b) Record time series with colours referring to the same bioclimatic groups as in (a). Pale lines represent the individual forest pollen records and bold lines the forest pollen composite for each group.

The Tropical composite features a relatively constant forest pollen percentage (between ∼ 40 and ∼ 65 %) across Termination V. However, this constant envelope is marked by variations such as a large drop between 432 and 427 kyr BP followed by an increase until 423 kyr BP and a smoother trough centred at 413 kyr BP. The Suptropical-N composite is characterised by an increasing trend in forest pollen percentages from 14 % at 433 kyr BP to 42 % around 419 kyr BP, mainly accounted for by the limited but rapid rise in the Valles Caldera pollen percentages around 425 kyr BP and later on by the pollen surge in the Ioannina record from 423 to 419 kyr BP. The Subtropical-N forests then develop relatively slowly, reaching a peak of 47 % towards the end of the termination at 409 kyr BP. Two additional pollen records support this Subtropical-N forest pattern. Lake Van pollen sequence (Litt et al., 2014), located on the eastern Anatolian high plateau (no. 22 on map Fig. 5a), depicts a progressive and moderate forest development reaching a ∼ 45 % maximum around 415 kyr BP. Similarly, the southeastern China sea record (no. 21, Fig. 5a, Cheng et al., 2018) features a reduced warm temperate forest development peaking around 410–400 kyr BP. The ES-BA group includes records from very different locations (Tibetan Plateau, Siberia, continental Europe, Japan). The resulting composite is first marked by a sharp and substantial increase in pollen percentages from 38 % at 429 kyr BP to 70 % at 419 kyr BP and then by a stagnating period (419–409 kyr BP), while the last part of the sequence features a slight drop at 407 kyr BP. The ES-BA composite increase is mostly driven by the strong development of Praclaux, Tenaghi Philipon and Lake Ohrid temperate forests and then by the growing Siberian forest (ICDP5011). Finally, the Subtropical-S forest pollen percentages decrease by over 13 % across TV. Overall, the present composite analysis highlights two interesting features. The warm temperate forests (Subtropical-N) continuously develop across the termination and reach a maximum around 409 kyr BP, while the temperate and boreal forests (ES-BA) feature a strong expansion over a short period (only ∼ 10 kyr) starting early in the termination (∼ 429 kyr BP). However, the limited number of pollen records also calls for both more observational data and for climate–vegetation simulations to gain further insight into the global terrestrial biosphere evolution.

3.3.2 Model–data comparison

The iLOVECLIM climate simulation of TV (Degla-V) provides a complementary view to that of the pollen composites, allowing a better evaluation of the forest development and its possible impact on atmospheric CO₂ changes recorded during TV. Figure 6 features the simulated tree fraction (Degla-V) and the forest pollen records at 433 kyr BP (glacial period, Fig. 6a) and 405 kyr BP (end of deglaciation, Fig. 6b). Based on numerous studies on the present-day pollen–vegetation relationships (Huntley, 1983, 1990; Zanon et al., 2018; Morales-Molino et al., 2020), we adopt pollen forest percentage as an indicator of the tree fraction.

Figure 6 Global comparison of simulated tree fraction (square patches, %) and forest pollen from records (circles, %) at 433 kyr BP (a) and at 405 kyr BP (b). Panels (c) and (d) feature the tree fraction difference between 405 and 433 kyr BP and 405 and 423 kyr BP respectively. Pollen record percentages should be compared to the simulated tree fractions of the closest adjacent land.

The absence of simulated trees at high latitudes (> 60^∘ N) at 433 kyr in Fig. 6a is consistent with the presence of a large ice sheet over North America, Greenland and Siberia during MIS12 (Batchelor et al., 2019). As expected, the high simulated tree fractions are located at lower latitudes and specifically around the Equator. However, there are two important and known model biases concerning simulated tree fractions (Goosse et al., 2010). The model features high Amazon forest along the eastern, northern and western South American coastline, while there is no reason why trees should be absent from inland regions of northern South America. This can be explained by a tube-like modelled Amazon river hindering moisture fluxes from the river to the atmosphere and resulting in a desert environment in the tropics. On the African continent, equatorial forest is shifted to the east with desert-like sub-Saharan regions. This anomaly derives from a poorly resolved African monsoon (Goosse et al., 2010).

The forest pollen observations overall support the simulated tree fraction while they also underline the aforementioned model biases both at the onset and the end of TV (Fig. 6a and b respectively). The northward gradient of simulated forest fraction over the Iberian Peninsula at 405 kyr BP (Fig. 6a) is consistent with lower forest pollen percentages in southern Iberia compared to northern Iberia supporting the idea that temperate forest is more developed than Mediterranean forest at the end of TV. At 433 kyr BP the simulated tree fraction over the Balkan region (<20 %) agrees with the low pollen percentages, and the simulated latitudinal forest gradient over the eastern Tibetan Plateau is in line with the two pollen records (Heqing and ZB13-C2). The northernmost African record (ODP 108-658) featuring a nearly absent Saharan forest is also well captured by the model both at 433 and 405 kyr BP. However, the model and observations disagree at 433kyr BP over equatorial Africa where the records describe a strongly rooted forest (forest pollen > 80 %), while the simulated tree fraction is very low, and over southern Africa where the model overestimates the forest proportions (for both 433 and 405 kyr BP), thus highlighting the aforementioned model biases. Finally, Japanese, Siberian and Amazonian simulated forests are relatively consistent with the corresponding pollen percentages.

In Fig. 6c we compute the difference between maps (b) and (a) in order to cancel out the systematic model biases and to illustrate the forest expansion/retreat over the whole deglaciation (405–433 kyr BP). As two pollen records (ICDP5011 and MD01-2447) do not cover the whole Termination, a second difference map (Fig. 6d, 405–423 kyr BP) complementary to Fig. 6c and including all sequences is proposed, although we mainly focus on Fig. 6c to reveal the full deglaciation forest trends. The impact of ice sheet melt is marked by strong forest development north of 60^∘ N in Siberia, surrounding the Bering Strait (confirmed by ICDP5011 record on a shorter timescale, Fig. 6d), and at 50^∘ N over the northwestern American continent. Model and pollen data agree over the eastern Tibetan Plateau and over Japan where the forest remains relatively stable. The simulation underestimates the forest expansion in the Iberian, the western Europe, the Balkan and the Amazonian regions and the forest retreat in southern Africa. Finally, we observe contrasting results over tropical Africa. The increasing simulated forest is in line with two west African deep-sea records (ODP 108-658 and GIK16867-2) but disagrees with recorded forest decline on the western (GIK16415-2) and the eastern (Lake Magadi) adjacent landmasses.

4 Discussion

Our integrated research approach from regional pollen record analysis to global climate simulation opens multiple discussion levels for both the observation and the modelling communities.

4.1 Assessing regional forest development across Termination V

In this section, we use th4e Degla-V simulation (in addition to the sparse, yet essential pollen records) to provide a finer temporal analysis of forest developments for selected regions of the world (Fig. 7). We choose eight regions based on special features emerging from Fig. 6c: either regional uniform trends (e.g. Siberia) or regional heterogeneous trends (e.g. central Asia). Because the selected regions are much larger than the catchment areas represented by marine and lacustrine records, the associated simulated forest features smoother variations than the pollen records. Therefore, it is more relevant to compare the simulated forest to the record mean when there is more than one record by region. For each region, we focus on the forest expansion trend and timing of the maximum. As already noticed for pollen composites (Fig. 4), we also observe regional discrepancies in the simulated tree fraction.

Figure 7 Comparison of simulated tree fraction and forest pollen percentage by region across TV. Each plot shows the simulated tree fraction average (%, solid line) over the corresponding regional box (coloured contours) and associated forest pollen record (%, pale dotted line) and record average (%, bold dotted line) when available. Locations of pollen records are marked by coloured stars.

The simulated Amazon tree fraction is relatively constant (∼ 60 %) and overall in line with the single record (Funza), although it does not reproduce the significant forest increase observed from 433 to 423 kyr BP. In tropical Africa, we observe a constant simulated forest growth across TV reaching a maximum of nearly 50 % at 411 kyr BP before decreasing. This region is characterised by contrasting pollen record signals. Although the record average follows a decreasing trend contradicting the simulated trend, it features a local forest expansion peak synchronous with the simulated one. Aside from a temporary development at 435 kyr BP, the North America (S) simulated tree fraction remains stable across the termination with increasing trends (east coast) compensating for decreasing trends (central North America, see Fig. 6c). The latter is confirmed by the pollen sequence of Valles Caldera in SW North America, which shows a relatively low development of forest cover with a decreasing forest trend from 420 to 405 kyr BP. The small Mediterranean region features an increasing forest trend both for the simulation and record average; however, this is much weaker in the simulation which can partly be explained by the null contribution of the northern Sahara in the model. Forest expansion reaches a maximum at ∼ 410 kyr BP both for the model and the data before shrinking at the end of the sequence for the modelled forest. Europe is characterised by a very weak increase in simulated forest coverage following a step-like pattern at ∼ 433, ∼ 427 and ∼ 420 kyr BP. The simulated tree fraction is relatively consistent with the pollen record average over the second part of the deglaciation (although it is ∼ 15 points below) in the sense that both reach a plateau around ∼ 420 kyr BP. However, they strongly differ at the beginning of the sequence. We propose two complementary explanations for this discrepancy. First, the observed strong temperate forest development at the onset of TV originates from southern Europe records only, yet Koutsodendris et al. (2010) argue that the central European forest expanded from ∼ 415 kyr BP only (based on the Dethlingen pollen record). Therefore, the record average is not fully representative of the European forest dynamics. Second, given the low horizontal resolution of the atmospheric module, we suspect that precipitation inaccuracies might account for a misrepresentation of growing conditions and an associated underestimation of simulated forest growth over Europe. For central Asia, the simulated tree fraction features a smooth and weak increase over TV with a slightly stronger trend from 435 to 430 kyr BP, consistent with the strong rise in the record average at 430 kyr BP. Finally, the two northernmost regions (North America (N) and Siberia) are totally unforested during the MIS12 glacial stage and both highlight a pronounced two-stage simulated forest development during the deglaciation. A first local maximum is reached at ∼ 423 kyr BP in both regions and a second and final peak arrives at ∼ 409 kyr BP in the North American (N) region and at ∼ 413 kyr BP in the Siberian region. Despite the lack of a precise pollen record to constrain the simulated forest evolution over North America (N), de Vernal and Hillaire-Marcel (2008) show, in a deep-sea sedimentary sequence collected off Greenland, that the highest pollen concentrations of Picea over the last 1 Myr occur between 410 and 390 kyr BP, thus supporting the idea of a strong forest expansion in the North American (N) region during TV. Concerning Siberia, the only available record (ICDP5011) features two maxima (∼ 423 and ∼ 413 kyr BP) consistent with the simulated development stages and a similar trend, although 20 to 30 points above the simulated one. We explain such a difference by the variable influence of ice retreat on forests depending on the latitude. The northernmost grid points of the region experience widespread growth (similarly to the ICDP5011 record) because of a net increase in available land, whereas southernmost grid points are mostly affected by changes in atmospheric temperatures, CO₂ concentrations and precipitation, resulting in lower forest expansion (see Fig. 6c).

This regional comparison of simulated forests and pollen record averages suggests different forest dynamics across TV. Over tropical Africa, the Amazon and North America (S), we find no evidence for clear forest trends. The Mediterranean region shows a significant forest growth mainly after ∼ 425 kyr BP and reaching a peak at ∼ 410 kyr BP. Finally, Europe, central Asia, North America (N) and Siberia share similar features. The forest expansion begins early (from ∼ 435 kyr BP) and reaches a first maximum around 425 kyr BP. The forest growth is particularly strong for the North American (N) and the Siberian regions (+40 % to +50 % over TV), supporting the argument for a massive boreal and temperate forest development. Such a significant forest development suggests a strong warming across TV and during MIS11 consistent with the reconstructed high temperatures, the reason for which Melles et al. (2012) refer to this period as a “superinterglacial”.

4.2 The potential role of boreal and temperate forests on atmospheric CO₂ levelling during TV

At a global scale, both pollen observations (ES-BA composite in Fig. 5b) and model simulations (North American (N) and Siberian tree fraction in Fig. 7) outline that TV is characterised by a strong increase in boreal and temperate forest between ∼ 429 and ∼ 423 kyr BP. We suggest that this result can support recent findings by Brandon et al. (2020). Brandon et al. (2020) show an abrupt increase in the global biosphere productivity from ∼ 429 and ∼ 427 kyr BP resulting in anomalously high values of productivity (10 % to 30 % higher than preindustrial), which would explain the atmospheric CO₂ levelling from ∼ 425 to ∼ 415 kyr BP (Fig. 8). This global biosphere productivity extreme is unique among the last five interglacials, and it is mainly attributed to terrestrial biosphere (Brandon et al., 2020). With forests as terrestrial productivity indicators, we seek to specify the geographical origin of such important carbon dioxide stabilisation. Among the four pollen composites, the Subtropical-N and the ES-BA feature the most important forest expansion (i.e. the strongest increasing trend, Fig. 8). However, the Subtropical-N composite features a lower trend than the ES-BA, which is also concomitant with the surge in the productivity signal (429–425 kyr BP). Moreover, it is important to note that the ES-BA composite represents a broader area than the Subtropical-N because it is composed of seven pollen records (versus four for Subtropical-N) covering a larger diversity of regions. The model simulations also confirm that the ES-BA pattern is present over the Siberian and North American (N) regions and the European and central Asian regions to a smaller extent (Fig. 7). Additionally, this significant development is synchronous with the surge in the productivity signal (429–425 kyr BP). Therefore, we propose that the temperate and boreal forests of the Northern Hemisphere are the major drivers of the high terrestrial biosphere productivity signal and the moderate atmospheric CO₂ values at the beginning of MIS11, while the contribution of the warm temperate forests of the Northern Hemisphere is weaker and arrives later (from ∼ 420 kyr BP).

Figure 8 Comparison of forest pollen composites and observed global biosphere productivity (inferred from Δ¹⁷O of O₂, solid green line, Brandon et al., 2020) and CO₂ (solid black line, Nehrbass-Ahles et al., 2020) over TV. Forest pollen composites from Fig. 5 are shown in dashed blue (ES-BA), yellow (Subtropical-N), purple (Tropical), red (Subtropical-S) and grey (record average) lines.

Figure 9 Zonal averages of simulated terrestrial biosphere (a) carbon stock (Gt C) and (c) fraction of total carbon stocks (%) across TV. Panel (b) features the total carbon stock of the simulated terrestrial biosphere and (d) the atmospheric CO₂ concentration from Nehrbass-Ahles et al. (2020). The simulated carbon stock is derived from the net productivity computed over each continental grid (Brovkin et al., 2002; see Sect. 2.3.1). The linear regression coefficient (Gt C kyr⁻¹) is provided for each zonally averaged time series (a and b). Panel (c) presents the total simulated terrestrial biosphere carbon stock (Gt C) across TV.

To further assess the carbon sink potential of temperate and boreal forests we analyse the variation in the simulated total carbon stock (i.e. carbon fluxes between the atmosphere and the vegetation + soil system), both a proxy for terrestrial biosphere productivity and a direct measure of atmospheric CO₂ removal (Fig. 9). Overall, the simulated total carbon stock (Fig. 9b) increases across TV, peaking around 410 kyr consistently with the average pollen-based forest maximum (grey line, Fig. 8) and in line with the simulated vegetation carbon maximum by Kleinen et al. (2016). Figure 9a represents the total simulated terrestrial carbon stock for four different zonally averaged vegetation types following the preindustrial Holocene biomes simulated by Prentice et al. (2011): temperate parkland and sclerophyll woodland ([−90, −30]^∘ N), tropical forest ([−30, 15]^∘ N), warm-temperate forest and sclerophyll woodland ([15, 40]^∘ N), and high-latitude boreal and temperate forest ([40, 90]^∘ N). The southern forest group features very low and slightly growing carbon stocks (brown, Fig. 9a). By contrast, the Tropical forest is the largest terrestrial carbon stock, accounting for approximately 50 % of the total carbon stock (dark green, Fig. 9c) throughout the termination. The relatively strong increasing trend in carbon stocks (∼ 8.1 Gt C kyr⁻¹, Fig. 9a) is segmented into two parts: a steep carbon stock rise at the onset of MIS11 from 436 to 425 kyr BP and a smoother rise from 425 to 410 kyr BP. The midlatitude temperate forest represents around 20 % of the total carbon stock (yellow, Fig. 9c). The slowly increasing stock becomes very sharp at 416 kyr BP and peaks at 500 Gt C around 411 kyr BP (Fig. 9a). Finally, the high-latitude boreal and temperate forest features the strongest mean increase in carbon stock (+13.7 Gt C kyr⁻¹, light green, Fig. 9a) among all regions. We observe a two-stage pattern (similarly to the North American (N) and the Siberian tree fractions, Fig. 7), with a steep rise until 425 kyr BP, along with the Tropical forest, giving way to a smoother trend. A second strong increase from 416 kyr BP leads to a peak value of 700 Gt C at 409 kyr BP. The first carbon stock maximum is synchronous with the biosphere productivity surge (Fig. 8).

While the total carbon stock increases by ∼ 900 Gt C until 410 kyr BP (Fig. 9c), Fig. 9b shows a redistribution of the carbon sinks mainly from the Tropical forests (∼ 54 % of total carbon stock at 436 kyr BP to ∼ 45 % at 404 kyr BP) to the high-latitude boreal and temperate forest ([40, 90]^∘ N) over the 436–420 kyr BP period (∼ 20 % of total carbon stock at 436 kyr BP to ∼ 28 % at 420 kyr BP) and to the larger area above 15 ^∘N over the 405–415 kyr BP period, indicating a northward increase in carbon fluxes. Therefore, we suggest that the high-latitude boreal and temperate forests account for the overall biosphere productivity increase leading to atmospheric CO₂ sequestration between 425 and 415 kyr BP. The warm temperate forests and sclerophyll woodlands, which include the Mediterranean forest, contribute more to the biosphere productivity over the 405–415 kyr BP period. The specific role of tropical forests remains unclear: the decrease in tropical forest percentages – mainly driven by GIK16415-2, Lake Malawi and Lake Magadi records – near 430 kyr BP leading to a local minimum around 428 kyr BP (Fig. 5) suggests that shrinking tropical forests could have temporarily contributed to the increase in CO₂ at 430 kyr BP, but we do not have any confirmation in the carbon stock analysis (Fig. 9).

Finally, the Tropical forest, the midlatitude temperate forest, and the high-latitude boreal and temperate forest are marked by a decreasing carbon stock from ∼ 408 kyr BP resulting in a total carbon loss of 200 Gt C in 404 kyr BP (Fig. 9b). This carbon stock decline is both consistent with the shrinking simulated forest in the North American (N) and Siberian regions (Fig. 7) and with the forest simulations by Kleinen et al. (2014). Additionally, we observe a similar decreasing trend in the biosphere productivity starting ∼ 3 kyr earlier (Fig. 7) concomitant with increasing atmospheric CO₂ concentrations.

5 Conclusions

Our new high-resolution pollen record (IODP U1386) allows us to resolve millennial-scale changes in southern Iberian vegetation cover over TV. This sequence shows a moderate Mediterranean forest expansion during TV (maximum MF pollen percentage of ∼ 50 %) in line with seminal work by Oliveira et al. (2016) and also unveils short timescale forest variations at the onset of the interglacial. TV is characterised by a two-step Mediterranean forest growth beginning at ∼ 424 kyr BP lagging behind the marine isotopic interglacial (defined by the global ice volume decrease) by approximately 4 kyr.

At a global scale, we assemble the first pollen record database for TV. Despite sparse pollen data, we observe different forest developments varying with latitude and highlighting the strong increase in temperate and boreal forests across the termination. We expect this new pollen database to open the way for extensive knowledge of terrestrial biosphere–climate interactions during TV. Yet, model–data comparisons remain fundamental to improve simulated forest accuracy in order to achieve a global cover. Here we observe that iLOVECLIM simulated forests feature overall similar but weaker trends than the pollen-based ones. This first assessment of global forest evolution during TV points out different driving mechanisms: a strong impact of ice sheet retreat on forests is observed at high latitudes, while lower latitude forests are probably more influenced by precipitation and temperature changes (Yin and Berger, 2012).

Based on the presented global pollen-based vegetation records and simulation analysis we suggest the following scenario:

1. The strong warming at the onset of MIS11 results both in a strong increase in high-latitude terrestrial biosphere productivity and in a rise in atmospheric CO₂ concentration probably driven by ocean physical and chemical degassing processes.

2. Because of the exceptionally high ice sheet melt characterising MIS11 (Dutton et al., 2015), forests can develop northwards. This strong expansion, mainly observed in the temperate and boreal forests, would allow carbon removal by terrestrial biosphere to partially compensate for the oceanic carbon losses, resulting in a long and unique CO₂ plateau from ∼ 425 to ∼ 415 kyr BP. Although the warm temperate forests of the Northern Hemisphere subtropical regions also respond to this warming, they might only play a significant role in the carbon cycle change around 410 kyr BP.

While there are yet too few data to depict the full impact of vegetation on the carbon cycle, we expect the aforementioned conclusions to be regarded as a first benchmark case for further assessment of global terrestrial biosphere during TV.

In that respect, we identify three major future research perspectives. First, a necessary step in order to achieve an accurate picture of terrestrial biosphere evolution during TV is to increase the pollen record coverage, especially in poorly documented regions such as the American continent, Siberia, and central and southern Asia (tropical forests). Second, an integrated understanding of the carbon cycle during TV will require both improving terrestrial biosphere simulations and taking into account the permafrost. Increasing the number of PFTs and applying spatial downscaling methods should be considered to get precise terrestrial biosphere representations. Similarly, investigating vegetation changes in terms of forest biomes will allow an understanding of forest replacement processes. Finally, comparative studies of terrestrial biosphere–climate interactions over contrasted terminations such as Termination II (TII) may be successful in determining the specificity of TV. In this view, cooperation between the modelling and observation communities remains essential.