Rising atmospheric CO₂ concentrations: the overlooked factor promoting SW Iberian Forest development across the LGM and the last deglaciation?

Across the last deglaciation, the atmospheric CO₂ concentration (CO₂) increased substantially from ∼ 180 to ∼ 280 ppm, yet its impact on vegetation dynamics across this major climatic transition remains insufficiently understood. In particular, Iberian pollen records reveal an intriguing feature that can be related to an often overlooked role of CO₂ in shaping vegetation responses during the last deglaciation. These records reveal the near disappearance of forests during the cold Last Glacial Maximum (LGM) and Heinrich Stadial 1 (HS1) phases and an unexpected recovery during the Younger Dryas (YD) cold phase when CO₂ increased. Here, we present high-resolution tracers of terrestrial (pollen, C₂₉ : C₃₁ organic biomarker) and marine (alkenone-derived Sea Surface Temperature, C_37 : 4 %, and long-chain n-alkanes ratios) conditions from the southwestern (SW) Iberian margin Integrated Ocean Drilling Program Site U1385 (“Shackleton site”) for the last 22 cal kyr BP. This direct land-sea comparison approach allows us to investigate how the Iberian Peninsula vegetation responded to major global CO₂ changes during the last deglaciation.

Our results show that cool and moderately humid conditions of the LGM supported a grassland-heathland mosaic ecosystem, but low CO₂ likely caused physiological drought and suppressed forest development. HS1, the coldest and most arid period, combined with sustained low CO₂ values, almost suppressed forest growth in favour of Mediterranean steppe. In contrast, the warmer Bølling-Allerød, characterised by a temperature optimum and variable but generally wetter conditions, along with the rise of CO₂ above 225 ppm at ∼ 15 cal kyr BP, contributed to substantial forest development. During the YD, sufficient moisture combined with increasing CO₂ enabled the persistence of a mixed grassland-forest mosaic despite cooler temperatures. Our study suggests that during cool and humid periods (LGM and YD) different pCO₂ values led to contrasting SW Iberian vegetation responses. In contrast, during periods of relatively high CO₂, temperature and precipitation played the main role in shaping the distribution and composition of the vegetation.

1 Introduction

The last deglaciation, spanning 20–19 cal kyr BP (e.g., Denton and Hughes, 1981; Toucanne et al., 2008; Denton et al., 2010) to ∼ 7 cal kyr BP (e.g., Dyke and Prest, 1987; Carlson et al., 2008) was marked by a global annual mean surface air temperature increase of 5 °C, (Annan et al., 2022), driving progressive melting of Northern Hemisphere glaciers. This interval was interrupted by an alternation of cold and warm phases (Alley and Clark, 1999): the warmer Bølling-Allerød (BA, 15–12.5 cal kyr BP) was bracketed by two major cold phases, Heinrich Stadial 1 (HS1, 18.5–15 cal kyr BP) and the Younger Dryas (YD, 12.9–11.6 cal kyr BP), particularly in the North Atlantic region. This period was also characterized by an increase in atmospheric carbon dioxide concentration (CO₂) from ∼ 180 to ∼ 280 ppmv (Monnin et al., 2001; Shakun et al., 2012; Marcott et al., 2014), one of the largest shifts of the last 800 000 years (Lüthi et al., 2008). This rise was not gradual, rather there were three main rapid (< 200 years) CO₂ rises, of ∼ 10 to 15 ppmv, at the end of HS1, within the BA and at the onset of the YD, as recorded in the West Antarctic Ice Sheet Divide ice core (Marcott et al., 2014).

Beyond its role in shaping global climate, pCO₂ directly influences plant physiology and the nature of vegetation responses to environmental change. The annual atmospheric-biospheric exchange of pCO₂ from photosynthetic activity is more than one-third of the total stored in the atmosphere (Farquhar and Lloyd, 1993). Atmospheric pCO₂ plays a critical role in plant physiology; plants absorb CO₂ through their stomata (small leaf pores) and lose water. At lower pCO₂, such as during glacial periods, plants must open these pores wider or increase their number to capture enough CO₂ (Royer, 2001), which whilst enhancing gas exchange, also increases water loss through transpiration, reducing water-use efficiency (WUE) and inducing physiological drought stress, even under moderate climatic conditions (Street-Perrot et al., 1997; Körner, 2000). Plasticity in stomatal conductance and density is considered an adaptive trait that evolved under declining Cenozoic CO₂ levels, enabling plants to sustain carbon uptake as concentrations approached glacial minima (∼ 180–190 ppmv), though at the cost of greater water loss (Wagner et al.,1996). These effects are especially pronounced in semi-arid environments. Under low pCO₂ conditions, species better adapted to drought and nutrient stress – such as those typical of steppes – are more likely to dominate and are typically observed in colder periods. Conversely, higher pCO₂ promotes forest expansion and higher plant productivity, particularly in trees that benefit from improved WUE (Huang et al., 2007; Randall et al., 2013). However, the response to pCO₂ is not globally uniform, with differences in water and nutrient availability, alongside other environmental constraints, mediating the response (e.g., Tognetti et al., 2008).While many reconstructions of past vegetation focus only on temperature and precipitation, the importance of pCO₂ as a limiting factor in plant productivity, coverage, and WUE is now widely supported by both empirical and model-based studies (e.g., Cowling and Sykes, 1999; Harrison and Prentice, 2003; Claussen et al., 2013; Piao et al., 2020). The BIOME4 model and a biome-scale reconstruction compiled from pollen records covering the last 40 kyr across the Northern Hemisphere (> 30° N) reveal a level of unexplained variability in vegetation patterns across both space and time (Cao et al., 2019), and provide key insights into other factors than temperature, precipitation and potential evapotranspiration that drive changes in vegetation dynamics and composition, such as pCO₂ (Ludwig et al., 2018; Cao et al., 2019). Similarly, coupled vegetation-climate modelling and multiproxy reconstructions demonstrate that CO₂ significantly impacts regional vegetation extent and productivity across glacial-interglacial transitions (Wu et al., 2007; Wei et al., 2021; Koutsodendris et al., 2023; Clément et al., 2024). These findings underscore the need to include CO₂ changes when interpreting pollen data or evaluating biome shifts (Prentice et al., 2017; Cao et al., 2019). While current day pCO₂ fertilization receives considerable attention (Piao et al., 2020), studying the effects of low CO₂, and major transitions from low to high CO₂, are equally critical.

Growing evidence shows that many climate reconstructions for glacial periods based on vegetation records may be biased as they neglect the influence of low pCO₂ on WUE. Neglecting this influence may contribute to the underestimation of past precipitation under full glacial conditions (Jolly and Haxeltine, 1997; Cowling and Skyes, 1999; Gerhart and Ward, 2010; Prentice et al., 2017; Cleator et al., 2020; Izumi and Bartlein, 2016; Chevalier et al., 2021), a concern still highlighted by recent studies (Chevalier et al., 2020; Wei et al., 2021; Prentice et al., 2022; Cruz-Silva et al., 2023). By contrast, under interglacial conditions with higher CO₂ levels, model experiments suggest that forest expansion in SW Iberia is mostly controlled by precipitation rather than by CO₂ levels (Oliveira et al., 2018, 2020).

Recognising the role of pCO₂ is crucial, not only to interpret the drivers of past ecosystems accurately, but also for anticipating how semi-arid landscapes – particularly in the Iberian Peninsula and Mediterranean region, which are predicted to experience significantly increased aridity – will respond to ongoing climate change. The terrestrial and marine climatic indicators from the SW Iberian margin sedimentary sequences provide an extremely valuable record of environmental change, and last deglaciation vegetation changes have been widely studied here using palaeoecological records (e.g., Peyron et al., 1998; Carrión, 2002; Chabaud et al., 2014; Combourieu Nebout et al., 2009; Dormoy et al., 2009; Fletcher et al., 2010a, b; Aranbarri et al., 2014; Bartlein et al., 2011; Naughton et al., 2011, 2019; Tarroso et al., 2016) and ecological niche modelling (Casas-Gallego et al., 2025). However, these are traditionally interpreted as a result of the combined effects of temperature, precipitation and evaporation changes, rather than pCO₂. Furthermore, few of the existing records span the entire deglaciation, or offer resolution or chronological precision to detect short-term vegetation and climate shifts in detail.

Here we present a new multiproxy study of IODP Site U1385 that allows a direct comparison between terrestrial and marine climatic indicators across the LGM and deglaciation at high (centennial-scale) temporal resolution. Hence, this record provides a detailed reconstruction of vegetation changes in SW Iberia along with sea surface temperature (SST) trends in its margin during the LGM, HS1, BA and the YD. Our new paleoenvironmental record enables the exploration of the main factors driving forest development during the LGM and the last deglaciation, and an evaluation of potential CO₂ thresholds for western Mediterranean forest development.

2 Materials and environmental setting

IODP Site U1385 is a composite record of five drillings in the SW Iberian margin (37°34.285^′ N; 10°7.562^′ W, 2587 m below sea level) located on a spur at the continental slope of the Promontório dos Príncipes de Avis, which is elevated above the abyssal plain and free from turbidite influence (Hodell et al., 2015) (Fig. 1). This work focuses on Hole A, a continuous record of 10 corrected revised meter composite depth (crmcd) mainly composed of hemipelagic silt alternating with clay (Hodell et al., 2015). For this study, Hole A was sampled from 3.84 to 1.08 crmcd, which corresponds to the period between ∼ 21.5 and 6.4 cal kyr BP. The sediment supply, including pollen grains, to Site 1385 is mainly derived via fluvial transport from the Tagus and Sado hydrographic basins, providing a reliable signature of the vegetation of the adjacent continent (Naughton et al., 2007; Morales-Molino et al., 2020).

Figure 1 Location of IODP Site U1385 and of the marine and terrestrial pollen records discussed in the text. Marine sedimentary records: MD03-2697 (Naughton et al., 2016); MD99-2331 (Naughton et al., 2007); MD95-2039 (Roucoux et al., 2005); D13882 (Gomes et al., 2020); MD95-2043 (Fletcher and Sánchez Goñi, 2008); MD95-2042 (Chabaud et al., 2014); SU81-18 (Turon et al., 2003); ODP Site 976 (Combourieu Nebout et al., 2009); Continental sedimentary records: Lake de Banyoles (Pèrez-Obiol and Julià, 1994); Quintanar de la Sierra II (Peñalba et al., 1997); La Roya (Allen et al., 1996); Ayoo de vidriales (Morales-Molino and García-Anton, 2014); Charco da Candieira (Van der Knaap and van Leeuwen, 1997); Padul 15-05 (Camuera et al., 2019). Black arrows represent the surface water circulation (MOW, Portugal and Azores Current). Note that coastline boundaries are for the present day. From: Schlitzer, Reiner, Ocean Data View, https://odv.awi.de (last access: 9 May 2020), 2019.

The present-day climate of southwestern Iberia is characterised by a Mediterranean climate strongly influenced by the Atlantic Ocean, Köppen classification Csa with warm summers (average ∼ 22  °C in the warmest month) mean annual temperatures 12.5–17.5 °C, and mean annual precipitation 400–1000 mm yr⁻¹. The rainy season peaks in the winter between November and January and drought occurs in the summer generally from June to September (Agencia Estatal de Meteorología (AEMET) and Instituto de Meteorologia (IM, Portugal), 2011).

The present-day vegetation of southwestern Iberia reflects a transitional biogeographical zone between temperate and Mediterranean climates (Rivas-Martínez et al., 2017). Coastal areas, influenced by oceanic humidity and milder winters, support thermophilous evergreen species such as Quercus suber, Olea europaea var. sylvestris, Myrtus communis, and Pistacia lentiscus (Asensi and Díez-Garretas et al., 2017). Inland, as elevation increases and oceanic influence diminishes, Mesomediterranean forests dominate, composed of both evergreen (Q. suber, Q. rotundifolia, Q. coccifera) and deciduous oaks (Q. faginea, Q. robur), often combined with heathlands or aromatic scrublands (e.g., Cistus spp.). Distinctive oak–juniper woodlands appear in drier zones, and pine forests (Pinus pinaster, P. pinea) are common on sandy coastal soils. Riparian zones feature Alnus glutinosa and Salix spp., while widespread Cistus and Erica shrublands reflect the area's susceptibility to fire.

3 Methods

3.1 Chronological framework

Sixteen AMS ¹⁴C dates were used to generate a new age-model for the last deglaciation at Site U1385 (Table 1 and Fig. 2). Five of these were published by Oliveira et al. (2018), based on monospecific Globigerina bulloides samples and analysed at the Vienna Environmental Research Accelerator (VERA), University of Vienna, Austria. A new set of eleven samples for AMS ¹⁴C analysis was selected primarily from monospecific assemblages of G. bulloides. When sample size requirements could not be met, a mixed assemblage of G. bulloides and G. inflata was used. All samples were processed at the Keck Carbon Cycle AMS Facility, University of California, Irvine (Table 1). The new age model, using solely the available radiocarbon dates for sedimentary record U1385 without any tuning, was calculated using a Bayesian approach through the software Bacon implemented in R (Blaauw and Christen, 2011; R Development Core Team, 2020) using the Marine20 calibration curve (Heaton et al., 2020). The studied interval encompasses the period from ∼ 22 to 6 kyr, as shown by the radiocarbon age model (Fig. 2). The average temporal resolution for the pollen and organic biomarkers across the deglaciation is 110 and 104 years, respectively, or slightly lower (174 and 135 years, respectively) when including the Holocene section (Figs. 3 and S1 in the Supplement).

Table 1 AMS radiocarbon ages of IODP Site U1385.

^* AMS from Oliveira et al. (2018).

Figure 2 Age-depth model for IODP Site U1385 using a Bayesian approach with Bacon v.2.3.9.1 (Blaauw and Christen, 2011). The original dates were calibrated using Marine20 (Heaton et al., 2020). Grey stippled lines show 95 % confidence intervals; red curve shows single “best” model based on the mean age for each depth. Upper graphs show from left to right: Markov Chain Monte Carlo (MCMC) iterations and priors (green line) and posteriors (dark grey line with a grey fill) for the accumulation rate and variability/memory. Note: the depth (y axis) was converted to cm from the corrected revised meter composite depth (crmcd).

3.2 Pollen analysis

A total of 97 samples (including 25 published by Oliveira et al., 2018) were analysed from 3.84 to 1.08 crmcd in Hole A, and prepared at the University of Bordeaux, France, using the standard protocol of the UMR EPOC laboratory (Georget et al., 2025). The sediment was firstly separated using coarse sieving at 150 µm, retaining the fine fraction. A sequence of chemical treatments, starting with cold HCl (hydrochloric acid) at increasing concentrations (10 %, 25 %, 50 %), eliminated calcium carbonate particles. Followed by cold HF (hydrofluoric acid) at increasing strength (45 % and 70 %) to eliminate the silicates. The remaining residue was micro-sieved (10 µm mesh), retaining the coarse fraction. Exotic Lycopodium spore tablets of known concentration were added to each sample to calculate pollen concentrations (Stockmarr, 1971). The obtained residue was mounted in a mobile medium composed of phenol and glycerol 1 % ($w/v$), to allow pollen/spore rotation and accurate identification. Samples were counted using a transmitted light microscope at 400X and 1000X (oil immersion) magnifications. To perform pollen identification, we used identification keys (Faegri et al., 1989; Moore et al., 1991), photographic atlases (Reille, 1992, 1995) and the SW Mediterranean modern reference collection.

The total count ranged from 198 to 1545 pollen and spores per sample, with a minimum of 100 terrestrial pollen grains and 20 pollen morphotypes to provide statistical reliability of the pollen spectra (McAndrews and King, 1976; Heusser and Balsam, 1977). The main pollen sum was calculated following previous palynological studies of Site U1385 (e.g., Oliveira et al., 2016) that excluded Pinus, Cedrus, aquatic plants, Pteridophyte and other spores, and indeterminable pollen. The pollen percentages are calculated against the main pollen sum; but the percentages of over-represented taxa were calculated based on the main sum plus the counts for that particular individual taxon; for example: 100× Pinus $/$ (Main sum + Pinus) and 100× Cedrus $/$ (main sum + Cedrus). Aquatic plants and spores were excluded because their abundant pollen originates in or near water bodies and can be transported far from their source, potentially overrepresenting regional vegetation. Pinus pollen, which is typically overrepresented in marine deposits, is transported by rivers from the Tagus and Sado's watersheds (Naughton et al., 2007). In contrast, the overrepresented Cedrus is transported by wind from the Atlas or Rif mountains in Morocco. Both taxa were also excluded from the main sum. PSIMPOLL 4.27 (Bennett, 2009) was used to plot percentages for selected taxa, grouped by ecological affinities (Gomes et al., 2020). Stratigraphically constrained cluster analysis by Sum of Squares determined the five statistically significant pollen assemblage zones (U1385-1 to 5 in Figs. 3, S1 and Table S1 in the Supplement) based on a dissimilarity matrix of Euclidean distances with pollen taxa ≥ 1 % (Grimm, 1987; Bennett, 2009).

Figure 3 Comparison of multiproxy records from Site U1385 with (a) 65° N July (W m⁻²) summer insolation (Berger and Loutre, 1991) (b) CO₂ (ppmv) composite from WAIS (Marcott et al., 2014). Principal pollen-based ecological groups percentages: (c) Temperate Mediterranean Forest from Site U1385 (solid black line) and compilation of Iberian Margin TMF records (D13882, MD03-2697; MD95-2042; MD95-2043; ODP-976; U1385) – GAM (grey curve), (d) Semi-desert taxa including Amaranthaceae (previously Chenopodiaceae), Artemisia, and Ephedra. (e) Heathland including members of the Ericaceae family (including various Erica spp.) and Calluna spp., (f) Mediterranean taxa and (g) Poaceae + Cyperaceae. Biomarkers: (h) C₂₉ $/$ C₃₁ ratio, (i) C_37:4 (%) and (j) SST ( °C). The different coloured shading corresponds to the pollen zones (Fig. S1 and Table S1 in the Supplement) which are related with the periods indicated above.

In addition to pollen-based ecological groups, we calculated the sum of Poaceae and Cyperaceae (Fig. 3g), to check the potential importance of C₄ plants in the Iberian Peninsula. While most of the present-day Poaceae and Cyperaceae in this region chiefly belong to the C₃ photosynthetic-pathway plant type (Casas-Gallego et al., 2025), it is possible that C₄ pathway plants were more important at other moments in recent Earth history. Pollen analysis is a core method in palaeoclimatology and palaeoecology, used to assess past climate conditions based on the ecological affinities of specific taxa grouped into pollen-based ecological groups. These groups reflect present-day vegetation–climate relationships, allowing inferences about dry, cold, warm, or moist conditions. As such, our pollen data reflect ecological responses rather than absolute quantitative climate parameters (Williams et al., 2001). A pollen diagram with clustering analysed (Fig. S1 in the Supplement) was produced revealing four main episodes over the LGM and the Last deglaciation (Fig. 3, further details in Fig. S1 and Table S1 in the Supplement).

3.3 Compilation of Iberian margin pollen records

In order to assess vegetation and climate changes more widely in the Iberian Peninsula region across the LGM and last deglaciation, we compiled available marine records along the Iberian margin covering the period from 23 to 6 cal kyr BP. Pollen count datasets from eight marine pollen records (D13882 – Gomes et al., 2020; MD03-2697 – Naughton et al., 2016; MD95-2039 – Roucoux et al., 2005; MD95-2043 Fletcher and Sánchez Goñi, 2008; MD95-2042 – Chabaud et al., 2014; ODP Site 976 – Combourieu Nebout et al., 2009; SU81-18 – Turon et al., 2003; Site U1385 – this study) were used with the original published chronologies, without any additional alignment or synchronization. Pollen percentages were recalculated against the main pollen sum. A uniform calculation of the pollen-based ecological group TMF (Temperate and Mediterranean forest) was made for each record, integrating the following taxa of (1) Temperate trees and shrubs: deciduous Quercus, Acer, Betula, Cannabis/Humulus, Carpinus, Castanea, Fraxinus excelsior-type, Hedera helix, Hippophae, Ilex, Juglans, Myrica and Vitis; and (2) Mediterranean taxa: evergreen Quercus, Quercus suber, Arbutus type, Buxus, Daphne, Jasminum, Ligustrum, Myrtus, Olea, Phillyrea, Pistacia, Rhamnus, Rhus.

To assess the general trend of vegetation patterns throughout the deglaciation, we applied a Generalised Additive Model (GAM), considered as a more robust statistical approach than loess curves (Wood, 2017; Simpson, 2018). The GAM model was fitted using the gam() function of the mgcv package (version 1.8.24; Wood, 2017) for R (version 3.6.3; R Development Core Team, 2020). We used a standard GAM with REML smoothness selection, specifying 30 basis functions (k=30) and a smoothing parameter of 0.0001 (sp = 0.0001). The relatively high K allowed the model to capture potential nonlinear patterns in the data without overfitting, while the small sp ensured sufficient smoothness; these values were chosen after exploratory analysis and diagnostic checks. To assess the validity of the smooth terms and confirm that the basis functions adequately captured the data wiggliness, we applied the gam.check() function of the mgcv package. The resulting k-index was greater than 1, and the p-value supported the hypothesis that enough basis functions were used. The fitted GAMs curves for TMF are presented along with approximate 95 % confidence intervals (Simpson, 2018).

3.4 Molecular biomarkers

Marine biomarker analyses were carried out in 123 levels, including 30 published by Oliveira et al. (2018). All analyses were performed following the extraction and analytical methods (Villanueva et al.,1997; Rodrigues et al., 2017). Marine coccolithophorid algae synthesise organic compounds including alkenones (Volkman et al., 1995) (Fig. 3i and j). Seawater temperature changes influence the amounts of di-, tri- and tetra-unsaturated alkenones produced by algae (Brassell et al., 1986). The use of organic solvents to separate the total lipid fraction from sediments allows the sea surface temperature alkenone-based reconstruction ($U_{\mathrm{37}}^{k^{\prime }}$ – SST) (e.g., Villanueva and Grimalt, 1997; Rodrigues et al., 2017). The $U_{\mathrm{37}}^{k^{\prime }}$ index (Prahl and Wakeham, 1987) was converted to temperature values using the global calibration equation defined by Müller et al. (1998) with an analytical uncertainty of 0.5 °C (Grimalt et al., 2001). Nevertheless, uncertainties remain, since $U_{\mathrm{37}}^{k^{\prime }}$ SST reconstructions may be affected by calibration biases, seasonal and ecological effects related to coccolithophorid production, and potential lateral transport or diagenetic alteration of alkenones (e.g., Conte et al., 2006; Ausín et al., 2022). Despite these caveats, the derived SSTs provide robust indicators of large-scale SST trends. Additionally, tetra-unsaturated alkenone (C_37:4) percentages were calculated due to their potential to identify the occurrence of cold freshwater pulses associated with iceberg discharges (Bard et al., 2000; Martrat et al., 2007; Rodrigues et al., 2011, 2017) and therefore, changes in the reorganisation of surface water masses in the North Atlantic (Rodrigues et al., 2017).

The ratio between C₂₉ and C₃₁ n-alkanes was also calculated to understand how epicuticular wax production in terrestrial plants varied through time. This index is generally considered to encompass the dynamic between woody plants and grasses plants of the adjacent continent (Cranwell 1973; Tareq et al., 2005, Bush and McInerney, 2013; Struck et al., 2020). This relation comprises the adaptation of plants by increasing the production of long-chain leaf wax, which reduces water loss during the photosynthetic processes and prevents desiccation promoted by harsh winds or more arid conditions (Bush and McInerney, 2013). Index values > 1 are typically considered to reflect higher quantities of C₂₉ n-alkanes produced by trees and shrubs, while values < 1 are generally considered to indicate higher quantities of C₃₁ n-alkanes by grasses and herbaceous plants (Cranwell, 1973; Rodrigues et al., 2010; Ortiz et al., 2010). However, the interpretation of this index may vary across biomes and depend on source vegetation types, and depositional processes (Carr et al., 2014; Diefendorf and Freimuth, 2017).

4 Results and discussion

4.1 The effect of CO₂ on biome changes during the LGM and deglaciation

Whilst a classic interpretation of ecosystem dynamics can be proposed solely considering variations in the main climatic parameters (temperature, precipitation), we hypothesise that changes in pCO₂ played an essential role in vegetation change, specifically in the deglacial forest expansion. Here, we evaluate the drivers of vegetation change by explicitly considering the evolution of CO₂ through the deglaciation. Our discussion is supported by the present-day environmental and climatic space, considering the temperature and precipitation in which different taxa exist in the Iberian Peninsula and characterising the TMF – Quercus sp., the Heathland (ERI) – Ericaceae family and the semi-desert (STE) landscapes (Fig. S2 in the Supplement).

4.1.1 Last Glacial Maximum (LGM, 23–19 cal kyr BP)

The pollen-based vegetation record from Site U1385 shows that during the LGM (pollen zone U1385-1: 21 500–17 990 cal yr BP, Fig. S1 in the Supplement) a grassland-heathland mosaic dominated the landscape, with semi-desert taxa (STE, ∼ 40 %) and heathland taxa (ERI, ∼ 10 %–20 %) (Figs. 3d, e, 4d), forming a distinctive non-analogue glacial vegetation cover.

Figure 4 Schematic representation of the relative change of climatic inferred parameters (precipitation and temperature) based on pollen-vegetation groups, biomarkers, SST as well as the physiological contribution of CO₂ for each period showing a schematic reconstruction of the potential ecosystem scenarios. The perceived temperature used the interpretation of pollen (TMF and STE groups), SST and n-alkanes; the plant available moisture (ERI, TMF and STE).

The prevalence of heath (Erica spp.) in Iberian pollen records underpins the classic view of the LGM in Iberia as a fairly humid interval, certainly compared with the extreme aridity of Heinrich stadials (Roucoux et al., 2005; Naughton et al., 2007; Fletcher and Sánchez-Goñi, 2008; Combourieu-Nebout et al., 2009; Sánchez-Goñi et al., 2009).

Nevertheless, the signals for moisture availability are somewhat complex. Semi-desert taxa, typically found in arid conditions, are abundant, while heathland taxa, associated with more humid environments, reach their maximum in the record (Figs. 3, S2c in the Supplement). Forest taxa are represented in low percentages (5 %–15 %) (Fig. 3c), suggesting cool and relatively dry conditions over the continent. The TMF values are consistent throughout the U1385 record and GAM-fitting to the data compilation (Fig. 3c). Similar patterns are observed across the marine records in southerly locations off the Iberian Peninsula (MD95-2043 – Fletcher and Sánchez Goñi, 2008 and ODP Site 976 – Combourieu Nebout et al., 2009 in the Mediterranean Sea, and SU81-18 – Turon et al., 2003 in the Atlantic Ocean), as well as further North off the Iberian Peninsula (MD99-2331 and MD03-2697 – Naughton et al., 2007, 2016).

Interestingly, the modern environmental space for the Ericaceae group (namely Erica arborea, E. australis, Calluna vulgaris) coincides with that occupied by the Quercus genus, the main constituent of the TMF group (Fig. S2b in the Supplement). This begs the question, if the environmental conditions that support heathland overlap with those for Quercus sp., then why were forests not thriving during the LGM? A possible explanation could be associated with cold atmospheric temperatures (SSTs average ∼ 14.5 °C, Fig. 3j), even if during the LGM the temperatures were not as extreme as the ones observed during the HS1 (Bond et al., 1993; Rasmussen et al., 1996). Hence, in addition to temperature, the lowest levels of CO₂ during the LGM (180–190 ppmv), could have been another important controlling factor, as they rank among the lowest concentrations recorded during the history of land plants (Pearson and Palmer, 2000; Tripati et al., 2009). Shao et al. (2018) modelled the global distribution of vegetation types as a function of temperature and precipitation under both modern climatic conditions and LGM CO₂ (185 ppm). Their analysis demonstrated pronounced qualitative differences in vegetation patterns, with reduced CO₂ levels grasslands were favoured to the detriment of evergreen broadleaf, evergreen and deciduous needle leaf forest. This study, however, did not include heathlands specifically, and it remain unclear whether this group has adaptations permitting better functioning under low CO₂ levels. We speculate that drought-adapted traits in Mediterranean Ericaceae especially E. arborea including thick cuticles, small leaf size, large photosynthetic thermal window and deep root system with large diameter and a massive underground lignotuber (Gratani and Varone, 2004) may have been advantageous under challenging trade-off between photosynthesis and water loss under very low pCO₂. Consequently, Ericaceae during LGM likely represented a vegetation group particularly resilient to such physiological constraints.

At the same time, the LGM coincides with a precession maximum, a configuration recognised to reduce seasonal contrasts (i.e., reduced summer dryness) and thus favour heathland development in the Iberian Peninsula, as documented in both glacial and interglacial contexts, including the Middle to Late Holocene (Fletcher and Sánchez Goñi, 2008; Margari et al., 2014; Oliveira et al., 2016, 2018; Chabaud et al., 2014; Gomes et al., 2020). Furthermore, heathland ecosystems thrive on acidic, low-nutrient soils, which can develop as a result of altered hydrological cycles during precession maxima. The ecological advantages of Erica also include less demanding edaphic requirements (low nutrient demand), more competitive re-sprouting strategy after disturbance, including fires, as well as a higher dispersal capacity compared with Quercus sp. for example (Pausas et al., 2008). However, these observations do not rule out a key impact of low CO₂ on vegetation composition during the LGM.

Diverse vegetation models have been used to understand the influence of climatic parameters and CO₂ during the LGM (e.g., Harrison and Prentice, 2003; Woillez et al., 2011; Izumi and Bartlein, 2016; Shao et al., 2018). However, there is a disagreement about the magnitude of the CO₂ influence, from being considered to have an equal influence (Izumi and Lézine, 2016) to being less critical than climatic parameters (Woillez et al., 2011; Shao et al., 2018; Chen et al., 2019). Harrison and Prentice (2003) also highlight model differences and the variable regional expression of the influence of CO₂ (with higher impact in tropical areas). However, these studies agree that low CO₂ had a negative physiological impact on forest development during the LGM in different continents (Jolly and Haxeltine, 1997; Cowling, 1999; Harrison and Prentice, 2003; Woillez et al., 2011; Shao et al., 2018; Chen et al., 2019). Jolly and Haxeltine (1997) used BIOME3 to simulate LGM vs pre-industrial CO₂ levels under different climatic conditions scenarios (temperature and precipitation) in tropical Africa; CO₂ was considered the primary driver of biome change from tropical montane forests to shrubby heathland ecosystems. This model included a photosynthetic scheme able to simulate plant response to different levels of CO₂ and its impact on stomatal conductance and water stress. This study showed that increasing CO₂ (above ∼ 190 ppmv), offsets the lower temperatures (changes of −4 to −6 °C), allowing the forest to thrive and replace heathland. However, plants with higher climatic demands (temperature and precipitation), which is the case of most temperate trees, are less competitive under low CO₂ conditions, compared with evergreen microphyllous species (e.g., Erica spp.).

Long-term studies considering CO₂ limitations on vegetation contrast in their perspectives; Gosling et al. (2022) argue that during the last 500 kyr, precipitation and fire exert the main controls on woody cover in tropical Africa while CO₂ effects were relatively small. In Asia, Clément et al. (2024) also emphasize the role of precipitation as the driver of vegetation distribution during interglacials, and that vegetation is not sensitive to CO₂ above 250 ppmv (value characterizing most of the interglacials); however, during glacial CO₂ conditions (< ∼ 185 ppmv), CO₂ is an important factor, favouring the increase of C₄ plants. The inclusion of CO₂ in climatic reconstructions for LGM for Africa and Europe yields a wetter LGM compared with reconstructions assuming CO₂ present-day concentrations (Wu et al., 2007). A similar impact is evident in the Last Glacial moisture reconstruction based on the pollen record of El Cañizar de Villarquemado in eastern Iberia; including a correction for the direct physiological effects of low CO₂, yields a wetter reconstruction of glacial climate (Wei et al., 2021). The implications of these experiments are important for the SW Iberian region and may help to solve the apparent contradiction between vegetation (abundance of semi-desertic plants and presence of heathland) and climate simulations, which indicate enhanced winter precipitation over southern Iberian and Northwest Africa due to southward shifting of the wintertime westerlies (Beghin et al., 2016). In the absence of CO₂ correction, temperature could also be misinterpreted; the LGM vegetation for Mediterranean sites was simulated and associated with warmer summers under LGM CO₂, instead of the colder conditions simulated with present-day levels of CO₂ (Guiot et al., 2000). In Europe, pollen reconstruction with steppe vegetation indicated warmer winter temperatures for LGM CO₂ compared with the modern CO₂ (Wu et al., 2007). The bias could extend to simulations of glacial vegetation; without the CO₂ effect, the cover of boreal and temperate forests is reduced, and evergreen forests are overestimated for the LGM (Woillez et al., 2011).

Experiments determining plant thresholds in response to low CO₂ have not received as much attention as research on the impact of high CO₂ levels (Gerhart and Ward, 2010; Dusenge et al., 2019), and to our knowledge no experimental work currently tests forest development under such low values. However, modelling approaches indicate that in C3 plants, photosynthetic capacity declines sharply once atmospheric CO₂ falls below ∼ 300 ppmv, making carbon assimilation increasingly limiting for plant growth (Wagner et al., 1996). When we assess the relationship between CO₂, SST and TMF across the LGM and deglaciation events we observe that the LGM (i) corresponds to SSTs below 15.5 °C and CO₂ below 225 ppmv, and (ii) that TMF values remain below 20 % (Fig. 5). In African mountain environments, a CO₂ threshold of approximately 220 ppmv has been suggested as the minimum above which forests could develop (Dupont et al., 2019). These results suggest that extremely low CO₂ below a critical threshold of ∼ 220–225 ppmv, played an important role in limiting forest development during the LGM. However, we emphasize that such thresholds should not be considered universal, as they may depend on plant taxa, edaphic conditions and microclimate. Nevertheless, these values, despite differences in baseline conditions such as insolation, are broadly consistent with other time intervals where Mediterranean forest expansion occurred, for example during MIS 13 at ∼ 216 ppmv (Oliveira et al., 2020) and MIS 18 at ∼ 215 ppmv under relatively high temperatures and increased winter rainfall (Sánchez-Goñi et al., 2023). Temperatures during the LGM in SW Iberia may have been sufficiently mild for forest development with sea surface temperatures of ∼ 15.5 °C (Fig. 3j) aligned with the broader threshold for forest development (Sánchez-Goñi et al., 2008). For this reason, one could speculate that a hypothetical increase in CO₂ above the observed critical threshold during the LGM could have permitted forest development in SW Iberia. Thus, while uncertainties remain, the convergence of multiple lines of evidence supports a key role for low CO₂ in constraining forest development specially during glacial periods.

Figure 5 Dispersion plot showing the relation between CO₂ (Marcott et al., 2014) and SST in relation to TMF % across the different intervals of the Last deglaciation, following the pollen zones boundaries.

4.1.2 Heinrich Stadial 1 (HS1, 18.5–15 cal kyr BP)

During HS1 (pollen zone U1385-2: 17 990–15 230 cal yr BP, Fig. S1 in the Supplement), a Mediterranean steppe landscape (Fig. 3d) with minimum arboreal development (Fig. 3c) corresponded to the lowest SSTs of the record (SST ∼ 12 °C, Fig. 3j). The dominance of semi-desert taxa and minimum TMF (Figs. 3 and 5c) indicate the highest levels of aridity. Additionally, high C_37:4 alkenone values (∼ 8.2 %, Fig. 3i) reflect major meltwater pulses, associated with severe cooling in the North Atlantic during HS1. Vegetation signals reinforce this picture of moisture stress. A sharp decrease in heaths (ERI, Fig.3e) and aquatic taxa such as Isoetes (Fig. S1 and Table S1 in the Supplement) suggests significant drying of terrestrial marshes and wetlands. The dominance of STE during HS1 is consistent with records across the Iberian Peninsula (Combourieu Nebout et al., 2002; Roucoux et al., 2005; Naughton et al., 2007, 2016; Fletcher and Sánchez Goñi, 2008) and is reflected also in a long-term minimum in reconstructed forest levels (Fig. 3c).

Throughout HS1, despite the gradual increase of CO₂ from ∼ 185 to ∼ 225 ppm between 18.1 to ∼ 16 cal kyr BP (Fig. 3b), this rise was insufficient to counteract the limiting effects of extreme cold and aridity. Simulations with the Weather and Research Forecast model that incorporate a CO₂ correction show a reduction in arboreal vegetation and an expansion of sparsely vegetated soils across Iberian region during HS1 compared with the LGM (Ludwig et al., 2018). The simulated precipitation values for SW Iberia (Tagus hydrographic basin catchment), remain below 700 mm yr⁻¹ for HS1, which agrees with the pollen evidence for widespread semi-desert taxa.

Interestingly, the differences between HS1 and LGM concerning temperature, precipitation and CO₂ are substantial. While CO₂ levels rose modestly, the climatic extremes of HS1 – marked by severe cooling and aridity – likely drove the observed loss of heathland and limited forest development across the peninsula. Moreover, based on pollen data from marine and terrestrial records we do not observe any significant (< 5 % TMF) latitudinal difference when comparing northern (e.g., Peñalba et al., 1997; Pèrez-Obiol and Julia, 1994; Roucoux et al., 2005; Naughton et al., 2007) with southern (e.g., this study; Comborieu Nebout et al., 2002; Fletcher and Sánchez Goñi, 2008) sites. Furthermore, the relationship between CO₂, SST and TMF across the HS1 shows scattered values of TMF (below 20 %) occurring at SST below 15.5 °C and CO₂ below 225 ppmv (Fig. 5), underscoring the combined climatic and physiological constraints on forest expansion.

4.1.3 Bølling-Allerød (BA, 15-12.5 cal kyr BP)

The BA (pollen zone U1385-3: 15230–12780 cal yr BP; Fig. S1 in the Supplement) had generally more favourable climatic conditions (higher temperatures, higher moisture availability) for TMF development (Fig. 3c) including a minor increase in thermophilous Mediterranean elements (Fig. 3c and f) and a reduction of STE (Fig. 3d). The combination of warming (SST above 16 °C, Fig. 3j) and a dry to wet trend are likely the primary drivers of progressive forest development during the BA. Additionally, the increase of CO₂ from ∼ 230 to 245 ppmv should have promoted a “fertilisation effect” during this time interval (Fig. 3b). BIOME3 simulations for African Biomes (Tropical Forest/Ericaceous scrub) with a present climate showed that above 190 ppmv, the increase of CO₂ at intervals < 20 ppmv, gradually offsets the negative effect of temperature changes. When CO₂ exceeds 250 ppmv with a temperature change of ∼ −6 °C the development of forests expands at expense of ericaceous scrubland (Jolly and Haxeltine, 1997).

Within age uncertainties, abrupt increases in CO₂ at 16.3 and 14.8 kyr (Marcott et al., 2014) (Fig. 3b) could be associated with the slight increase of forest at the onset of the BA and the subsequent highest peaks of forest development observed during the BA, respectively (Fig. 3c). Cao et al. (2019), using pollen-based biome reconstruction, suggested that worldwide expansion of forests was a consequence of the increasing CO₂ superimposed over the temperature increase between 21 and 14 kyr. Cao et al. (2019) further emphasise the role of CO₂ after the LGM driving a general northward expansion of forests and replacement of grassland by temperate forests in Europe. During the BA, considering that temperature and moisture availability in SW Iberia was favourable, increases in CO₂ levels (> 225 ppmv) may have amplified TMF expansion during this period (Figs. 4b and 5).

4.1.4 Younger Dryas (YD, 12.9–11.7 cal kyr BP)

The YD (pollen zone U1385-4: 12780–11190 cal yr BP, Fig. S1 in the Supplement) is characterised by an initial weak forest contraction followed by its progressive expansion (Fig. 3c). At the regional scale, the landscape likely consisted of a forest-grassland mosaic, as suggested by the relatively high presence of forest elements coexisting with semi-desert taxa (Figs. 3c, d and 4a). Strong SST cooling (Fig. 3j), (equivalent to LGM SSTs or even cooler), with a minimum of 13.2 °C, without significant freshwater pulses, may have been associated with cooler land surface temperatures. However, this impact may have been muted by the positive effect of higher moisture availability (based on the presence of TMF, Naughton et al., 2019) and/or the increasing trend of CO₂ (Fig. 3b). The weak reduction in TMF observed in our record and corroborated by the compiled records (Fig. 3c) contrasts with the steppe environment often described for this interval, especially in the southeast of the Iberian Peninsula (Carrión, 2002; Camuera et al., 2019). A more pronounced forest contraction is observed in the high-altitude terrestrial/lacustrine cores (Quintanar de la Sierra II – Peñalba et al., 1997; and La Roya – Allen et al., 1996) in which the near-disappearance of the forest might reflect the altitudinal adjustments in vegetation belts (Aranbarri et al., 2014). However, the U1385 and numerous Iberian records (e.g., Lake de Banyoles – Perez-Obiol and Julià, 1994; MD03-2697 – Naughton et al., 2007; MD95-2039 – Roucoux et al., 2005; Charco da Candieira – van der Knaap and van Leeuwen, 1997; MD95-2042 – Chabaud et al., 2014; D13882 – Naughton et al., 2019; MD95-2043 – Fletcher and Sánchez Goñi, 2008; ODP Site 976 – Combourieu-Nebout et al., 2002) show a relatively high percentage of TMF during the YD compared to HS1 (Fig. 3c).

Unfortunately, there is a lack of independent precipitation proxies for SW Iberia, and Denniston et al. (2018) highlight the limited reliability of speleothem proxies as indicators of precipitation in this region for this time interval. More widely in the Iberian Peninsula, a double hydrological structure with a drier first phase and wetter second phase was proposed, the latter favouring the expansion of mountain glaciers (García-Ruiz et al., 2016; Baldini et al., 2019). We observe that the notable YD forest development occurred, counterintuitively, in association with similar SSTs to those of the LGM and only slightly higher than those of HS1. Along with higher summer insolation, higher CO₂ (> 240 ppmv, Fig. 5) may have been a key factor in supporting forest development. A climate simulation from transient experiments using LOVECLIM, for the site SHAK06-5K/MD01-2444 located near U1385, obtained a weaker AMOC, colder winter temperature, and lower precipitation for the YD compared with the LGM (Cutmore, 2021; Cutmore et al., 2021). This supports the scrutiny of additional factors, notably pCO₂ influence on moisture availability for plants, to explain the substantial levels of TMF observed in the Iberian margin records (Fig. 3c). The increase in pCO₂ may have enhanced plant productivity and WUE (Cowling and Sykes, 1999; Ward, 2005) during the YD, partially compensating for the impact of atmospheric cooling and drying. Schenk et al. (2018) suggest pCO₂ may play an essential role in forest development if enough moisture is available. Tree cover may have been confined to suitable moist microhabitats and areas close to refugia; however, it was clearly less restricted than during earlier cold periods (Svenning et al., 2011), as indicated by the TMF abundances (Fig. 3c). Simulations from vegetation-climate models based on pollen records for biome reconstruction (Shao et al., 2018) and in a dynamic vegetation model (ORCHIDEE) driven by outputs from an AOGCM (Woillez et al., 2011) emphasise the role of increasing CO₂ as a critical factor for global forest development during the period including the YD (Shao et al., 2018). Underlying these changes, the increase in summer insolation (Fig. 3a), which contributed to the increase of summer temperatures and winter precipitation in the Mediterranean region (Meijer and Tuenter, 2007), cannot be neglected as a driver of forest development, at least where trees were not excessively water-stressed. However, disentangling the contribution of insolation vs CO₂ requires sensitivity experiments, not yet performed. In summary, the persistence of TMF during the YD, despite colder winters and drier summer conditions compared to the BA, is most plausibly explained by the combined interaction between precipitation variability, maximum insolation and increasing CO₂ (between ∼ 245 and 265 ppmv) (Fig. 4a).

4.1.5 Early to Middle Holocene (11.7–4.2 cal kyr BP)

Pollen zone U1385-5 (11190–4260 cal yr BP) corresponds to the Early to Middle Holocene. This interval is marked by the expansion of TMF and thermophilous Mediterranean elements, reflecting a regional increase in temperature and precipitation in parallel with warm SSTs (> 18 °C). Despite coarser temporal resolution for this interval, the U1385 record is consistent with nearby records showing a maximum forest development at ∼ 9000 cal yr BP (Fig. 3c). The specific timing of the Holocene forest maximum varied across the Iberian Peninsula along a gradient of regional moisture availability (Gomes et al., 2020). The Early Holocene CO₂ exceeded 260 ppmv, representing full interglacial conditions. The combination of coupled interglacial ocean-atmosphere conditions (reflected in high SSTs) and high CO₂ supported maximum forest development (Fig. 5). The impact on plant moisture availability compared to the preceding glacial conditions would have been profound, supporting high productivity and further increases in WUE. The progressive lifting of CO₂ constraints on photosynthesis throughout the Last Deglaciation may thus represent an important factor underlying forest development in SW Iberia.

4.2 C₂₉ $/$ C₃₁ ratio and C₃ $/$ C₄ dynamics: potential and limitations

Long-chain n-alkanes with odd-numbers, such as C₂₉, C₃₁, C₃₃, are epicuticular waxes produced by terrestrial plants. In lake sediments, higher abundance of C₂₉ occurs in catchments with more trees, while higher abundance of C₃₁ is observed in grassy catchments (Meyers, 2003). However, caution is required when interpreting C₂₉ $/$ C₃₁ in taxonomic terms because both woody plants (trees and shrubs) and grasses can produce C₂₉ and C₃₁ chain lengths (Ortiz et al., 2010; Bush and McInerney, 2013). Additionally, n-alkane chain-length patterns differ across species and environments, so C₂₉ $/$ C₃₁ ratios cannot be interpreted as strict woody vs grass markers (Bush and McInerney, 2013). Insights into the dominance of different plant physiological pathways in response to contrasting levels of CO₂ and humidity can be potentially gained by analysing C₂₉ $/$ C₃₁ n-alkanes from Site U1385. The C₂₉ $/$ C₃₁ ratio shows important variability between climatic phases, with increasing values during the LGM, high values during HS1 and the YD, and lower values during the BA and Holocene (Fig. 3h). The C₂₉ $/$ C₃₁ ratio is positively correlated (Pearson's correlation coefficient, r = 0.52, p-value = 2.473 × 10⁻⁸) with the STE pollen group and negatively correlated (r = −0.63, p-value = 2.821 × 10⁻¹²) with TMF (Fig. 3c, d and h), indicating a link between pollen-based vegetation changes and n-alkane chain-length distributions. Notably, the expected simple interpretation of C₂₉ $/$ C₃₁ as “trees vs grasses” does not appear to hold for this dataset across the different phases. We propose two possible interpretations, which may explain the observed C₂₉ $/$ C₃₁ variability. First a physiological stress response hypothesis: the C₂₉ $/$ C₃₁ ratio in this setting may reflect an adaptation of plants to aridity. Leaf-waxes n-alkanes protect plants against the loss of water during the photosynthesis (Post-Beittenmiller, 1996; Jetter et al., 2006). We could expect that arid, cold and windy conditions impose greater physiological stress on woody plants than on grasses. Consequently, increases in the C₂₉ $/$ C₃₁ values during HS1 and YD, may suggest that woody plants (TMF and ERI) responded to climatic stress by enhancing leaf wax C₂₉ production as a protective strategy (Fig. 3h). Second, changes in C₂₉ $/$ C₃₁ ratio may indicate compositional changes in vegetation, particularly between woody-dominated communities with diverse ecological tolerances – from semi-desert dwarf shrubs such as Artemisia to mesophyll broad-leaved trees. In this context, a prevailing “trees vs grasses” interpretative framework may not adequately represent the Iberian vegetation patterns. We stress that both hypotheses are plausible but require further validation. Nevertheless, the coherent climate signal observed in the U1385 record is encouraging for future studies aimed at linking leaf-wax chemistry of contributing species and vegetation dynamics in this region (Cutmore, 2021; Cutmore et al., 2021). Beyond taxonomic shifts, the link between n-alkane chain-length distributions and the C₃ $/$ C₄ plant dynamic is also relevant. African savannahs are dominated by C₄ plants, and C₃₁ n-alkanes have been associated with their presence in past landscapes (Dupont et al., 2019). Worldwide, 80 % of Poaceae (grasses) and Cyperaceae (sedges) use a C₄ photosynthetic pathway that is favoured by arid conditions (Sage, 2017). However, pollen analysis cannot discriminate between Poaceae and Cyperaceae pollen morphotypes exclusively or in its majority of C₄ plants. We have grouped the Poaceae and the Cyperaceae pollen taxa, noting the inherent limitations of this grouping to represent C₄ plants in Iberia as we know that C₃ are the dominant grasses across the region (Casas-Gallego et al., 2025) (Fig. 3g). Accordingly, the grouping of Poaceae + Cyperaceae pollen must be interpreted with caution. During the last deglaciation, this group (Poaceae + Cyperaceae) shows relatively high but fluctuating values between the LGM and the BA stabilising thereafter, without a clear correlation with other proxies (TMF, STE, or C₂₉ $/$ C₃₁). Thus, we do not find strong evidence for an increased role of grasses/sedges, or of C₄ plants specifically, during arid or low CO₂ intervals. Experimental studies, show that C₃ grasses outperform C₄ grasses when temperatures rise by 5 to 15 °C at a low CO₂ concentration of 200 ppm. Research on the quantum yield of photosynthesis identified a “crossover temperature” – the point at which C₃ and C₄ plants perform equally. This crossover depends on both temperature and CO₂ levels. Modelling across 0–45 °C and CO₂ levels from 150–700 ppm shows that whether C₃ or C₄ plants are favoured is determined by the interaction between these two factors; however, humidity was not considered in these studies (Ehleringer et al., 1997; Edwards et al., 2010). Most C₄ plants are confined to tropical grasslands and savannahs; where they are better adapted to higher temperatures, arid, nutrient poor conditions, and environments subject to intense disturbance caused by animals or fire regimes (Bond et al., 2005; Edwards et al., 2010). Likewise, one should expect that vegetation in SW Iberia after the LGM (Figs. 3 and 5) should be dominated by C₃ plants. This interpretation is consistent with the relatively cold SSTs (Fig. 5) and the high percentages of Artemisia spp. (C₃ plant) in the pollen record (Fig. S1 in the Supplement). However, it is not possible to completely rule out a role of C₄ plants in the glacial vegetation of SW Iberia, because pollen morphology does not allow the separation of these groups. Stable isotope studies on ancient grass pollen were able to discriminate C₃ $/$ C₄ grasses, but single-grain analyses remain technically challenging (Nelson et al., 2016). Further biomarker research is therefore needed to resolve C₃ $/$ C₄ dynamics in temperate – Mediterranean biomes. Current evidence suggests that the C₃ $/$ C₄ shifts documented in African savannahs (e.g., Dupont et al., 2019) are not directly applicable to our study area. Species-level biomarker fingerprinting will be essential to test whether C₄ plants played a significant role in Iberian Mediterranean ecosystems during the last deglaciation.

In summary, although future isotopic and biomarker approaches hold great promise for resolving C₃ $/$ C₄ dynamics, the current evidence strongly supports C₃ dominance in SW Iberia during the deglaciation. This interpretation is consistent with the modern distribution of plants in the region, where less than 10 % of grasses are C₄ (Casas-Gallego et al., 2025), and with the prevailing cool and relatively humid conditions of the LGM and YD, which favour C₃ over C₄ photosynthesis. Thus, while acknowledging the limitations of pollen-based proxies, the available data indicate that C₃ plants were the dominant contributors to the Iberian vegetation signal.

5 Conclusions

This study presents high-resolution pollen and biomarkers records from Site U1385 off the SW Iberian Margin, providing new insights into vegetation dynamics during key climate transitions of the last deglaciation and the associated CO₂ changes. We applied a biomarker proxy (leaf wax C₂₉ $/$ C₃₁ ratio), which is positively correlated with the semi-desert pollen curve and negatively correlated with TMF, demonstrating its potential as an indicator of aridity in marine cores of the Western Mediterranean region. The high temporal resolution of the record, combined with a robust radiocarbon chronology, enables consistent and more accurate comparisons with regional datasets, strengthening its contribution for future palaeoenvironmental reconstructions and model simulations.

Rather than simply interpreting our dataset in terms of past temperature and precipitation changes, we examine the U1385 record in the context of modern and palaeo observational and modelling studies that support a significant influence of CO₂ on past vegetation distribution and composition. During the LGM, cool temperatures, low seasonality, and physiological drought stress under low CO₂ restricted forest growth and favoured heathlands. Traits of Mediterranean Ericaceae, such as deep roots and thick waxy leaves, may have given these plants a competitive advantage. During HS1, woody vegetation was significantly suppressed due to cold and arid conditions, exacerbated by low CO₂. The subsequent expansion of temperate Mediterranean forests during the BA was promoted by warmer, wetter conditions and favoured by rising CO₂ concentrations. During the Younger Dryas, despite a return to cooler temperatures, forest-grassland mosaics persisted supported by adequate moisture availability and sustained higher CO₂ levels.

Furthermore, our study suggests a critical CO₂ threshold for forest expansion at ∼ 225 ppmv. Below this value (e.g., LGM and HS1), arboreal populations were generally restricted in their development with the impact of aridification and cooling being detrimental. Above this value, forests expanded (e.g., during the BA) and the effects of adverse climatic conditions (e.g., during the YD) was buffered. This threshold is consistent with observations from Mediterranean to tropical African environments (e.g., Dupont et al., 2019; Oliveira et al.,2020; Koutsodendris et al., 2023; Sánchez-Goñi et al., 2023). The concept should be further tested in regional vegetation models to determine the vegetation response to CO₂ fluctuations during past cold periods.

Finally, our findings highlight the importance of CO₂ as a key driver of vegetation change in the Mediterranean region through its control on plant moisture-availability and water-use efficiency (Koutsodendris et al., 2023). These palaeodata provide critical context for understanding vegetation responses under future climate scenarios with rising CO₂ and shifting precipitation regimes. At the same time, they underscore the need for further research on the relationship between long-chain n-alkanes, vegetation types, and C₃ $/$ C₄ plant dynamics, as the long-chain n-alkanes do not yet provide a reliable basis to disentangle the dynamics between woody plants and grasses in the Mediterranean ecosystems.