Introducing relative pollen productivity estimates for Iberian taxa: methodological insights and implications for landscape modelling in the Western Mediterranean

Understanding the impact of ongoing global change on plant communities requires long-term quantitative reconstructions of past vegetation dynamics. Fossil pollen records offer one of the most powerful tools to reconstruct past landscapes, yet for their accurate interpretation it is important to take into account the differential pollen productivity of plant taxa. For southern Europe, and particularly for the Iberian Peninsula, estimates of pollen productivity remain scarce, limiting our ability to refine palaeoecological reconstructions.

Here we present the first relative pollen productivity estimates (RPPs) for 21 common taxa in continental Spain. For that purpose, we used 1113 modern pollen samples from our own surveys and the Eurasian Modern Pollen Database (EMPD2), and vegetation data from the Spanish Forestry Map (MFE) and the Iberian and Macaronesian Vegetation Information System (SIVIM). RPPs were derived by applying an optimisation algorithm with the REVEALS model (REgional VEgetation Estimates from Large Sites). To test the reliability of our RPPs, we validated 8 arboreal taxa in 27 present-day coretops across Spain. We also compared the obtained RPPs with different studies across Europe, using a bias-free comparison framework.

Our findings indicate that the dominant arboreal taxa (Pinus, evergreen and deciduous Quercus) are high pollen producers, whereas temperate forest, shrub and herbaceous taxa generally yielded medium to low estimates of pollen productivity. Validation of arboreal taxa from present-day coretops showed that REVEALS-based estimates perform better than raw pollen counts when compared with present-day vegetation cover. Comparison between different studies in Europe also showed that most of the Spanish RPPs are similar to those obtained in Europe, although notable differences emerged for some taxa.

This study calculates, validates and compares the first RPPs in the Western Mediterranean, highlighting the value of quantitative palaeoecological data for Holocene landscape reconstructions. The findings of this paper would support that the Iberian Peninsula could have been home to a heterogeneous mosaic of open areas, conifers and broadleaf trees, offering new frameworks to improve both palaeoecological reconstructions and contemporary forest management strategies. Note: nomenclature and authority of the taxa in this paper were selected according to Flora Iberica (Castroviejo, 1986–2012).

1 Introduction

Science for mitigation and adaptation to global change needs to quantify how much landscapes changed under the pressure of climate variability and human agency. Acquiring a numerically detailed understanding of changes in land use and vegetation cover through time is crucial to establish reliable environmental models of the impacts of past climate change on plant communities.

Likewise, land cover reconstructions are crucial to reconstruct changes in land use (Fyfe et al., 2015; Pearce et al., 2023; Roberts et al., 2018; Trondman et al., 2014; Woodbridge et al., 2019). A central challenge in palaeoecology is to determine not only how vegetation cover varied in space and time, but also how these changes interacted with the physical environment and disturbance regimes. This involves asking key palaeoecological questions shaped by some of the long-standing debates over the historical structure and dynamics of European landscapes. Some influential and controversial hypotheses (e.g., Svenning, 2002; Vera, 2000) dispute the traditional view of densely forested wilderness (Bradshaw and Mitchell, 1999; Ellenberg, 1988; Peterken, 2001). Instead, they propose that large herbivores and disturbance regimes maintained extensive areas of semi-open habitats across much of postglacial Europe. This debate continues today, supported by robust palaeoecological data used to evaluate these hypotheses (Pearce et al., 2023). Recent studies increasingly suggest that postglacial Europe had a mixed composition of both closed forests and open areas (Carrión et al., 2010a; Nikulina et al., 2024; Pearce et al., 2025a), highlighting that the drivers of vegetation openness extend beyond climatic and human influences and emphasising the significant roles of large herbivores (Pearce et al., 2025b).

Accurate reconstruction of forest structure and baseline conditions for ecological restoration largely depends on pollen-based land cover reconstructions, which are not only vital for unravelling the long-term interplay between ecological and human processes, but also for producing reliable regional and global climate models and biogeochemical cycles (Abrantes et al., 2012; Cheddadi et al., 1998; Li et al., 2011; Liu et al., 2023). Integrating quantitative pollen data with diachronic cartography and multi-proxy evidence is also crucial to reconstruct the effects of past disturbances (Githumbi et al., 2022; Pirzamanbein et al., 2014, 2020; Zanon et al., 2018) – including fire, human deforestation and herbivory – allowing us to quantify biomass affected over time and the spatial extent of disturbance regimes and to disentangle overlapping effects of natural and anthropogenic drivers (Ellis, 2021; Morrison et al., 2021; Nikulina et al., 2024; Pearce et al., 2025b).

In this regard, pollen records are the most widely used proxy records for past vegetation reconstruction (Andersen, 1970; Broström et al., 2008; Davis, 1963; Huntley, 1990; Sugita, 1994; Von Post, 1918). Indeed, one of the main goals of pollen analysis is to reconstruct past plant abundances. However, it has long been known that there is a lack of linearity between pollen presence and abundance of the producing plant taxa (Andersen, 1970; Davis, 1963; Prentice and Parsons, 1983; Sugita, 1994), resulting in some taxa being overrepresented in fossil pollen records due to their high productivity and effective dispersal, while others may be underrepresented owing to low productivity and limited dispersal capacity (Davis, 1963). This discrepancy can lead to biased reconstructions of past vegetation and land cover (Prentice and Webb, 2009; Sugita, 1994).

Thus, rigorous estimation of the composition of past vegetation relies on our ability to better comprehend and quantify the relationships between fossil pollen assemblages and the composition of the vegetation that produces them.

A first step in the quantitative reconstruction of Quaternary vegetation using fossil pollen records requires calculating the relative pollen productivity estimates (RPPs/PPEs) for the taxa whose land cover we aim to reconstruct. Pollen productivity is often defined as the number of pollen grains produced per unit relative abundance of a given taxon, and is usually expressed as a dimensionless ratio relative to a reference taxon, since absolute pollen production measurement is difficult to determine (Andersen, 1970). Relative pollen productivity estimates (hereafter referred to as RPPs) are one of the critical parameters required to produce a reliable model of past vegetation abundance, as they enable correction of under or over estimations of taxa abundance.

RPPs have been calculated for many regions of northern and central Europe in the last decades (Abraham and Kozáková, 2012; Baker et al., 2016; Broström et al., 2004; Bunting et al., 2005; Grindean et al., 2019; Hjelle, 1998; Kuneš et al., 2019; Mazier et al., 2008; Nielsen, 2004; Niemeyer et al., 2015; Poska et al., 2011; Soepboer et al., 2007; Theuerkauf et al., 2013). RPPs have also been calculated in the eastern Mediterranean (Ergin et al., 2024), as well as in North America (Calcote, 1995; Chaput and Gajewski, 2018; Commerford et al., 2013), Africa (Duffin and Bunting, 2008; Tabares et al., 2021), Asia (Han et al., 2017; He et al., 2016; Jiang et al., 2020) and Oceania (Mariani et al., 2016, 2022). Despite the abundance of RPP studies for different taxa in mid or high European latitudes, other regions of the world remain understudied. One example is the five Mediterranean-climate regions (MCRs) of the world. Despite occupying less than 5 % of the Earth's surface, they host about 48 250 known vascular plant species (Cowling et al., 1996) and yet, there are very few RPPs existing in these areas (Ergin et al., 2024; Githumbi et al., 2022). The MCRs also represent a critical biome for understanding long-term human-landscape relationships as the Mediterranean Basin in particular represents a very long history of human occupation and therefore ecosystem change, resilience or persistence, besides a particularly vulnerable scenario regarding current global change and future warming (IPCC, 2023). The lack of RPPs for MCRs precludes any quantitative land cover reconstruction from fossil records in these areas, although efforts have been recently made to obtain new RPPs for some Mediterranean areas (Ergin et al., 2024; Githumbi et al., 2022; Serge et al., 2023). The fact that all previous numerical approaches have been conducted in mixed temperate forests or in subtropical areas implies that the complexity of Mediterranean plant communities has rarely been considered in these attempts.

Some of the available RPPs for northern European taxa could potentially be of use in the Mediterranean Basin, but it is well-known that pollen productivity might be driven by a number of geographical factors, plant taxonomy and climate constraints (Baker et al., 2016; Broström et al., 2004, 2008). Often, RPPs for the same taxa may differ due to methodological issues at the data resolution and landscape characterisation level, and from a number of methodological assumptions (Liu et al., 2022). These inconsistencies challenge the transferability of RPPs beyond their original context and, consequently, their application is often restricted to localised settings (Liu et al., 2022). Essentially, having a robust estimate of pollen productivity for taxa of a particular region implies considering all these factors and thus obtaining new RPPs.

In the present work, our objectives are: (1) to produce RPPs for 21 woody and herbaceous taxa across the Spanish Territory of Iberia (STI, from now on), (2) to validate the obtained results using present-day pollen samples from coretops within the region and (3) to compare our results with those from other European RPP studies. STI holds one of the greatest ecosystem, habitat and plant species diversities of Europe (Maestre et al., 2021; Médail and Quézel, 1999; Mutke et al., 2010) and comprises two biogeographical domains: the Mediterranean region, which accounts for about 70 % of the STI, and the Eurosiberian region, located in the northernmost areas of the STI (∼ 30 %). Our work is the first comprehensive study conducted in Iberia to obtain relative pollen productivity estimates, and it represents the first step towards the quantitative reconstruction of past landscapes framework in the Western Mediterranean, and particularly in STI.

2 Material and methods

2.1 Study area

The study area covers the whole of continental Spain (504 782 km²) including the two biogeographical regions, Eurosiberian and Mediterranean.

2.1.1 Mediterranean biogeographic region

The Mediterranean region extends throughout most of the STI, covering all of the central and eastern region except for the mountain and alpine areas. STI exhibits remarkable climatic diversity, largely driven by its complex topography and geographic position. Mean annual temperatures (MAT) range from approximately 8 °C in the interior plateaus and mountainous regions to around 18 °C along the Mediterranean coast (Fig. 1). In alpine zones of the Mediterranean (Southern Spain), MAT often falls below 8 °C. Precipitation patterns are equally heterogeneous: while the national mean annual precipitation (MAP) is around 500 mm, values vary widely – from 1000–1500 mm in mountainous areas to as low as 200–600 mm in coastal and plateau regions (Chazarra et al., 2018) (Fig. 1). This climatic variability underpins the country's exceptional ecological and floristic diversity, resulting in a variety of habitats and landscapes.

Figure 1 (a) Mean annual temperature (MAT, °C) and (b) mean annual precipitation (MAP, mm) in continental Spain (Ninyerola et al., 2005). Black dots represent the locations of the surface samples used in this study; green diamonds represent the coretops used for validation. Black dashed line delimits the bioclimatic regions, Eurosiberian in the North, Mediterranean in the South.

Mediterranean forest ecosystems cover two thirds of the total wooded region in the STI (Costa et al., 1998). These woody communities are physiognomically diverse and vary from shrub to dense mature forests, and from thorny, macchia-like temperate steppes to cold semi-deserts (Costa et al., 1998; Gavilán et al., 2018). Mediterranean forests, structured along a marked altitudinal gradient, are predominantly monospecific, although they are occasionally mixed with other woody species. These are mainly evergreen sclerophyllous taxa, though sometimes deciduous taxa are also present (Fig. 2).

In coastal and lowland areas (0–400 m a.s.l.) macchia and garrigue-type shrublands, grasslands and forests of different species of Pinus (P. halepensis Mill., P. pinaster Aiton, P. pinea L.) are dominant. Some other pines can also be found in the Mediterranean foothills and mountain belts (400–1200 m a.s.l.): P. sylvestris L. (Scots pine), P. nigra subsp. salzmanii (Dunal) Franco, P. nigra subsp. nigra J.F. Arnold (black pine) and those from lowland areas, along with the main oak taxa of Mediterranean sclerophyllous forests, Quercus ilex L. and Q. suber L. (holm and cork oak woodlands, Fig. 2). These forests appear often combined with Juniperus spp. communities (J. thurifera L., J. phoenicea L., J. communis L., J. oxycedrus L., or J. sabina L.), which are quite characteristic of the plateau-continental STI, including open areas that rarely form continuous canopy forests. Oak, pine and juniper communities are all adapted to periods of aridity that may vary from two to nine months of the year. Some woodlands in the Mediterranean mountains also support deciduous taxa such as birch (Betula pendula Roth., B. alba L.) and ash (Fraxinus excelsior Vahl., F. angustifolia L.), as well as some semi-deciduous oak species such as Q. faginea Lam. and Q. pyrenaica Willd. Relics of Abies pinsapo Boiss. are also found in the mountain ranges of Southern Spain.

Mountain vegetation (>1200 m a.s.l.) often presents continuous arboreal cover as, over the last decades, forest recovery in previously managed montane regions has produced denser forest communities. Yet, these montane-subalpine regions also host patchy plant communities where forests blend with open ecosystems. Some species of juniper (J. phoenicea and J. thurifera) as well as mountain and Scots pine (P. sylvestris, P. uncinata Ramond ex DC. in Lam. & DC.) are present. This vegetation belt is also characterised by sparse shrubs and grasses.

Cultivated olive trees (Olea europaea L.) are very extensive in the southern half of Spain (see Fig. S1 in the Supplement), and often appear in the wild in shrubby habitat in the eastern half of Spain. The olive tree is a key Mediterranean taxon that existed in Iberia before domestication, represented since at least the Upper Pleistocene in continental records (i.e., Fernández et al., 2007; González-Sampériz et al., 2020) and even the Early Pleistocene in the marine cores of Portugal and Spain (Magri et al., 2017), but started to expand during the Early Holocene (Langgut et al., 2019), with cultivation beginning during the late Middle Holocene and intensifying during the last 4000 years (Carrión et al., 2010b; Martín-Puertas et al., 2008).

Figure 2 Main forest types in the Iberian Peninsula (MITECO, 2024). Black dots represent the locations of the surface samples used; red diamonds represent the coretops used for validation. Black dashed line delimits both bioclimatic regions.

2.1.2 Eurosiberian biogeographic region

The Eurosiberian or Atlantic region lies in the northern area of the STI, including the Cantabrian range and the Pyrenees. Climatic conditions are characterised by cold winters and mild summers, with MAT ranging from 5 °C in the mountainous areas to 14 °C in the coastal areas (Fig. 1). Drought periods are shorter with abundant and well-distributed rainfall throughout the year, with MAP ranging from 1000 to 2000 mm.

Forest composition differs considerably from that of the Mediterranean region. In coastal areas (0–300 m a.s.l.), vegetation has been severely affected by invasive alien species (“Non-native species communities” group in Fig. 2), especially by Pinus radiata D. Don and Eucalyptus spp. L'Hér., both of which were cultivated for reforestation and industrial timber production but are now naturalised. Forests dominate in montane areas (300–1000 m a.s.l. in the Cantabrian range; 1000–1600 m a.s.l. in the Pyrenees). These forests are deciduous or semi-deciduous mixed communities of oaks (Q. robur L., Q. petraea (Matt.) Liebl., Q. pyrenaica), beech (Fagus sylvatica L.), birch (Betula alba, B. pendula), ash (Fraxinus angustifolia, F. excelsior), hazel (Corylus avellana L.) and other mesic taxa (“Mixed broadleaf communities group in Fig. 2) that rarely form monospecific communities (Sorbus aria (L.) Crantz, S. aucuparia L., Acer monsspesulanum L., A. opalus Mill., A. campestre L., Tilia cordata Mill., T. platyphyllos Scop., Juglans nigra L., J. regia L., Castanea sativa Mill.). Atlantic mountains support conifers such as Pinus sylvestris or P. uncinata (Fig. 2), or other conifers such as Abies alba Mill. (silver fir), which often form mixed forests with beech in the Pyrenees (Fig. 2). Subalpine vegetation (1600–2400 m a.s.l.) is characterised by sparse shrub and grasses, which are heavily grazed by livestock, although conifers such as Scots and mountain pines can still be present, as well as the silver fir. Mountain pine marks the alpine tree line (>2400 m a.s.l.), above which only grasslands and cushion plant communities can resist the severe climatic stress found at these altitudes.

2.2 Quantitative pollen-vegetation relationships

Brief overview of existing methods

Pollen-vegetation models have expanded over the last decades from those based solely on pollen/vegetation ratios (Davis, 1963), linear regressions and extended R-values (ERV) (Parsons and Prentice, 1981; Prentice and Parsons, 1983; Sugita, 1994), to the Landscape Reconstruction Algorithm (Sugita et al., 2010) and the most recent Bayesian models (Dawson et al., 2016; Garreta et al., 2012; Liu et al., 2022; Veeken et al., 2022).

One of the first attempts to develop appropriate techniques that account for pollen productivity and dispersal was made by Davis (1963): the R-value model. This model assumes that pollen deposition rates are directly proportional to abundances, each taxon having a characteristic constant of proportionality, i.e., an R-value (Prentice and Parsons, 1983). The R-value would designate the ratio between pollen percentage and vegetation percentage for each taxon. Parsons and Prentice (1981) developed the Extended R-value (ERV) method, introducing two submodels (ERV-1 and 2) to overcome the statistical limitations of ratios, site-to-site variability, and the effects of long-distance pollen transport. ERV-1 expressed the background component (non-local pollen) as “a constant background pollen percentage for each taxon” while ERV-2 expressed it as “a constant proportion of total forest volume (or whatever measure of abundance is being used)”. A third submodel (ERV-3) assumes a constant absolute amount of background pollen deposition for each taxon at all sites (Sugita, 1994), since the correlation between pollen and vegetation will not improve further beyond a certain distance, introducing the concept of the Relevant Source Area of Pollen (RSAP). As a modification to ERV-3 to address site-to-site taxon variability, Theuerkauf and Couwenberg (2022) developed ERV-4, expressing the background component as a result of the pollen productivity multiplied by the distance-weighted regional plant abundance for each taxon. ERV-1 and 2 use pollen and vegetation percentages, whereas ERV-3 and 4 use pollen percentages and plant abundance data expressed in absolute abundances. In our study area, the application of ERV approaches, whilst possible, would seem logistically challenging due to the large spatial extent, rendering the need for more flexible methods.

Sugita (2007a, b) proposed a new framework for vegetation reconstruction, the Landscape Reconstruction Algorithm (LRA), consisting of two different steps: the REVEALS and LOVE models. The REVEALS (REgional VEgetation Estimates from Large Sites) model is designed to reconstruct regional vegetation composition over large spatial scales (typically >10⁶ ha) by correcting for biases in pollen representation due to differences in pollen productivity and dispersal. It uses pollen data from large lakes or multiple small sites to estimate the relative abundance of plant taxa in the surrounding landscape. This approach accounts for differential pollen production and transport, making it more robust than simple pollen percentage analyses. LOVE (LOcal VEgetation Estimates), on the other hand, focuses on reconstructing vegetation at smaller spatial scales (<10⁴ ha), integrating regional vegetation estimates from REVEALS with local pollen data to separate local vegetation signals from regional background. This two-step framework allows for a hierarchical understanding of vegetation patterns, from broad regional trends to fine-scale local dynamics. The LRA model was developed as an inverse process for the ERV model. Both hold the same assumptions, especially when the REVEALS model uses several small sites rather than one big lake. The most important one is the even (homogeneous) pattern of vegetation (Sugita et al., 1999). This implies that at all sites, on average, all species at all distances should be the same, including the regional vegetation. This can be fulfilled in simulated landscapes, but real landscapes potentially have very large vegetation patches produced by mountains or vegetation belts (Abraham et al., 2014; Fang et al., 2019; Kuneš et al., 2019). If their size is about half or one third the size of the region and if the sampling area falls uniquely into one of these patches, species representation across different distances produces a non-stationary pattern.

More recently, pollen-vegetation models have been parameterised by using Bayesian hierarchical models, which have the primary goal of accounting for the uncertainty of pollen dispersal and production (Dawson et al., 2016; Liu et al., 2022; Paciorek and McLachlan, 2009). These Bayesian approaches, contrary to the ERV, REVEALS and LOVE models, simultaneously estimate pollen productivity and dispersal by finding the parameter values that best explain the sediment pollen data given the known vegetation cover, adapting better to spatial complexity and making them more suitable for regions with diverse vegetation and topography. However, Bayesian approaches to estimating pollen productivity and dispersal are challenging to apply in regions devoid of nearly-continuous sampling strategies (Dawson et al., 2016; Liu et al., 2022; Trachsel et al., 2020). These models, as those from the ERV, typically assume relatively homogeneous forest structure, consistent vegetation composition, and access to fine-scale, spatially resolved vegetation data - conditions that are rarely met in complex landscapes like the Iberian Peninsula. In Iberia, vegetation is highly patchy and spans distinct bioclimatic zones, and available data sources vary in resolution and taxonomic detail, making it difficult to replicate the fine-scale separation between local and regional vegetation required by such models. Additionally, the small size and topographic complexity of most Iberian lakes, along with the presence of unique taxa with distinct dispersal traits, further complicate the direct application of these Bayesian methods. In short, these particularities limit the transferability of the standard Bayesian framework and highlight the need for regional recalibrations.

Our selection to compute RPPs in the Iberian Peninsula

In this study, we have obtained RPPs relative to Poaceae by using the inverted REVEALS method that finds optimal RPP values by comparing REVEALS vegetation estimates from recent pollen spectra with current regional vegetation (Fig. 3). This approach was developed by Abraham et al. (2014) as manual adjustments of RPP values in the performance of the REVEALS model. The next application of inverted REVEALS was made by Kuneš et al. (2019), by integrating the global optimisation algorithm DEoptim (Mullen et al., 2011), which finds the RPP value iteratively and gains higher precision than manual adjustments. The algorithm begins by generating initial candidate RPP values for each taxon, which are then used to compute REVEALS estimates from pollen and vegetation data. The distance between the REVEALS estimates and the observed vegetation is then calculated with a loss function (Eq. 1). To identify the best-fitting RPP values, we iteratively adjusted them to minimise the loss function, using the Generalised Simulated Annealing algorithm, as implemented in the GenSA package in R (Xiang et al., 2013) (Fig. 3). GenSA belongs to a class of stochastic global optimisation techniques particularly well-suited for navigating complex, multidimensional parameter spaces characterised by numerous local minima. Unlike traditional optimisation approaches that may become trapped in suboptimal solutions, Generalised Simulated Annealing leverages probabilistic transitions to explore the solution space more broadly and escape local optima. GenSA is implemented with a C++ core, ensuring computational efficiency and scalability for large ecological datasets.

Figure 3 The optimisation loop to obtain RPPs. The algorithm starts with trial RPPs for each taxon (1), then uses candidate RPPs to estimate REVEALS from pollen data (2), and calculates distance between modelled and observed vegetation with a loss function (3). By using GenSA, the algorithm iteratively adjusts RPPs to minimise the loss function (4) until optimal RPPs are found (5).

GenSA identifies the optimal set of taxon-specific RPP values that minimises the discrepancy between observed regional vegetation composition and vegetation proportions reconstructed by the REVEALS model. The objective function minimised during optimisation is the weighted sum of squared errors (WSSE) between modelled and observed vegetation proportions across all regions. This approach aligns with the weighted least squares (WLS) regression (Carroll and Ruppert, 1988), although we included a smoothing offset in the denominator to reduce the influence of low-abundance taxa and to avoid division by zero:

$$\begin{array}{}\text{(1)}& L\left(\mathit{a}\right)={\sum }_{r=\mathrm{1}}^R{\sum }_{i=\mathrm{1}}^m{\displaystyle \frac{{\left(V_{i,r}^{\mathrm{mod}}-V_{i,r}^{\mathrm{obs}}\right)}^{\mathrm{2}}}{V_{i,r}^{\mathrm{obs}}+\mathrm{1}}}\end{array}$$

where:

• a: vector of RPP to be estimated

• R: number of analysed grids

• m: number of taxa

• $V_{i,r}^{\mathrm{mod}}$: vegetation proportion for taxon i in region r, reconstructed by REVEALS

• $V_{i,r}^{\mathrm{obs}}$: = observed vegetation proportion of taxon i in region r.

Retrieving RPPs with the optimisation loop involved first setting the initial range of possible pollen productivity values to 0.1–50 (Fig. 3), meaning that the optimisation function would search for 500 possibilities before displaying the RPPs that gave the smallest distance between observed and reconstructed vegetation proportions. Moreover, we set 500 iterations to ensure the decrease between each run, and bootstrapped 100 resampled versions of all the sites, in order to obtain error estimates. Since such a setup of the parameters requires high computational efforts (500×100 runs per taxon are needed to retrieve a single RPP), we compiled the optimisation function for a better performance by using the “compiler” package, included in base R (R Core Team, 2025). We also parallelised the optimisation process using “foreach” (Folashade et al., 2009) and “doParallel” (Folashade et al., 2011) packages. Full code workflow is provided in the Supplement.

Regarding pollen dispersal, this parameter is crucial for reconstructing past vegetation abundances, although in practice it is not possible to reliably measure long-distance dispersal of airborne particles (Katul et al., 2005). While models for predicting airborne particle dispersal are needed (Kuparinen et al., 2007), the correct model selection according to the studied area is also key. In the present work, we compute the Gaussian Plume and the Lagrangian Stochastic models (GPM and LSM, respectively) to calculate the dispersal and deposition factor. GPMs are the most widely used dispersal models, which describe dispersal patterns based on observations of particles in the atmosphere (Jackson and Lyford, 1999). However, it is known that GPMs are based on observations of the dispersion of larger objects than pollen, and their application to pollen transport can provide an incomplete representation of the full, long-distance, dispersal pattern, underestimating dispersal of heavy pollen grain taxa like Abies (Abraham et al., 2014; Kuparinen, 2006; Mariani et al., 2016; Theuerkauf et al., 2013, 2016; Theuerkauf and Couwenberg, 2020). The weighting factor decreases their distance-weighted plant abundance, and therefore the pollen signal is counterbalanced with high pollen productivity. Since large basins strengthen overestimation of Abies RPP with the GPM model (Abraham et al., 2014), studies with different basin sizes lead to different RPPs for Abies. The LSM (Andersen, 1991; Kuparinen et al., 2007) is a mechanistic model that describes wind and turbulence conditions in the atmosphere, simulating the trajectory of single pollen grains over short and long distances (10¹–10⁵ m) (Theuerkauf et al., 2013), below, in and above the canopy (Kuparinen, 2006; Theuerkauf, 2025). This paper focuses on LSM-derived RPPs, since they showed lower standard deviations than the GPM-derived RPPs, as well as a better performance during the validation process for Abies and Juniperus. Nevertheless, the results of the GPM-RPPs and their validation are available in the Appendices A and B.

Pollen fall speeds (Tables S1 and S2, Supplement) were retrieved from the literature or calculated following Stoke's law for spherical particles and with Falck's assumption for ellipsoidal grains (Gregory, 1961) using the photographs contained in identification guides (Reille, 1992, 1995).

In short, we have applied a framework using the inverted REVEALS approach (Abraham et al., 2014; Kuneš et al., 2019; Theuerkauf, 2025), i.e., we estimate REVEALS before calculating RPPs, then validated the results using modern coretop samples and compared them with other pollen productivity estimates from Europe (Fig. 4).

Figure 4 Methodological workflow of this study. Pollen data (from surface and coretop samples) and vegetation data were used (1) to calculate present-day REVEALS estimates. Using both a Lagrangian Stochastic dispersal model and the Gaussian Plume model, the GenSA optimisation algorithm was applied (2) to identify the optimal set of RPPs relative to Poaceae, minimising the difference between observed and REVEALS-estimated vegetation proportions (3). The resulting RPPs were validated against present-day coretops (4) and compared with values obtained in other studies (5), following the numerical workflow described in Sect. 2.3 “RPPs comparison across European studies”.

2.2.1 Pollen and vegetation data acquisition

The taxa chosen for RPP computations include the most frequent arboreal types in both present-day STI forests and palynological sequences (Carrión et al., 2022) to which we can attribute pollen types: Abies Mill., Betula L., Corylus L., Fagus L., Pinus L., Olea L., deciduous and evergreen Quercus L. Shrub and herbaceous taxa from frequently represented families and genera in STI vegetation and fossil pollen records have also been included: Amaranthaceae/Chenopodiaceae Juss., Artemisia L., Asteraceae. SF. Asteroideae Lindl., Brassicaceae Burnett, Asteraceae. SF. Cichorioideae Chevall., Erica L., Genista-type (Genista L. and Ulex L.), Juniperus L., Plantago L., Poaceae Juss., Ranunculaceae Juss., Rosaceae Juss. and Rumex L.

Pollen data

Producing the necessary RPPs for STI requires analysis of a large number of sampling points of both modern pollen and vegetation cover. Therefore, we divided the study area into 25 grids of 150×150 km² each (Fig. S2 in the Supplement). The percentage of pollen for the main taxa was then retrieved for each grid, and compared with regional vegetation data.

We used modern pollen counts from the Eurasian Modern Pollen Database 2 (hereafter referred to as EMPD2) (Davis et al., 2020). According to the EMPD2, the region of continental Spain has the second highest number of samples, with 1110 modern pollen samples taken from terrestrial moss polsters, soils or lake sediments.

We reviewed the sample context and excluded those originating from marine and estuarine environments, resulting in the selection of 1113 samples, of which 70 were obtained through our own field surveys, conducted at various times over the past 30 years. RPPs were then obtained by excluding 51 modern coretops that have been used as a validation set (see Fig. 4 and Sect. 2.2.2 “Validation of RPPs”).

Vegetation data

Present-day arboreal vegetation cover for the 150×150 km² grids was obtained through the most recently published database of the Spanish Forestry Map (MFE) at a scale of 1:25 000 (MITECO, 2024). MFE classifies the vegetation cover into plant communities for which detailed information on the coverage of the 3 main woody, arboreal or shrub species is provided. We then estimated the forest patch size by multiplying the percentage cover of each main species by the area of the polygon in which it is located.

For shrub and herbaceous taxa we used the relevé-based database Information System on Iberian and Macaronesian Vegetation (SIVIM) (Font et al., 2017). We used all the available plant inventories from SIVIM, amounting to a total of 149 646 surveys where we performed taxa harmonisation according to pollen types (Tables S1 and S2, Supplement).

We derived information on crops and other land uses using the CORINE land cover (CLC) database (European Environment Agency, 2019) so all vegetation types could be analysed, especially in areas where human-modified landscapes are dominant. We spatially intersected CLC polygons with MFE to classify each territory unit according to both land use type and forest cover presence/absence. This dual classification differentiates between areas designated as forest by CLC that indeed retain actual forest cover and those that have undergone deforestation. The resulting landscape matrix was then combined with SIVIM data, which incorporates shrub and herb taxa, to quantify the floristic composition across different landscape contexts and forest cover conditions. Additionally, CLC provided critical information on Olea europaea crop abundance, which is absent from MFE as olive typically represents agricultural rather than natural forest systems.

2.2.2 Validation of RPPs

Validating REVEALS-based vegetation estimates using modern forest composition and raw pollen data requires an independent dataset. Accordingly, we excluded 51 coretop samples from the EMPD2 dataset when deriving RPPs (Fig. 4) to reduce potential circularity. This allowed us to reconstruct vegetation proportions for those sites using our new RPPs, and compare them with actual forest cover. From the 51 coretops validation dataset we chose 27 samples by excluding salt lakes, where the surface samples are often subject to aeolian erosion, and samples under closed canopy or from high elevations where the pollen signal might be biased. The coretops for validation were selected from a variety of landscapes and can be found as a map compilation in Fig. S3 of the Supplement.

We validated RPPs on those 27 samples using the REVEALSinR function from the “disqover' package (Theuerkauf et al., 2016) to account for the productivity and dispersal-deposition biases. Data input requires: (i) pollen counts at each lake site from which the coretop comes; (ii) estimates of pollen productivity and fall speed of pollen for all taxa; (iii) standard errors of the RPPs; (iv) distance-weighting method, for which both the LSM and GPM (Prentice-Sugita dispersal-deposition model) (Sugita, 1994) were selected; (v) basin type (peatland or lake) and diameter in meters for each and; (vi) diameter of the reconstructed region in meters. REVEALS-based estimates of the modern samples were then compared with rings of 15, 30, 45 and 100 km radius of present tree cover around each sample. Only arboreal taxa (disaggregated by taxa) were selected, since herb and shrub taxa data are devoid of surface area information. Regional plant cover for validation was obtained from the MFE, except for Olea (olive crops) which was calculated from the CORINE Land Cover dataset.

Evaluation of the validation process was conducted using a multimetric approach based on four error and bias metrics for each taxon: root mean square error, mean absolute error, mean absolute percentage error and normalised mean bias. The root mean square error and the mean absolute percentage quantify absolute deviations; the mean absolute percentage error expresses deviations relative to observed values, and the normalised mean bias captures systematic over- or underestimation. Each metric was calculated as the difference between raw pollen estimates and REVEALS-derived values, standardised using z-scores, and averaged per taxon to produce a composite improvement score – a dimensionless measure of validation performance.

Given the large volume of validation results, this paper focuses on the coretop validation of LSM-RPPs at 45 km resolution, which showed the strongest performance according to the multimetric analysis. The rest of the results regarding the LSM-RPPs validation can be found in the Supplement (Figs. S4–S6).

2.3 Comparison of RPPs across European studies

Comparing RPPs between studies has become a challenging task due to the large number of existing studies which all make partially different assumptions. We introduce a numerical workflow designed to address such a challenge, enabling more reliable comparisons and addressing methodological discrepancies between studies (e.g., different reference taxa).

Table 1 List of studies included in the comparison. Full metadata of the references are available in Table S3 in the Supplement.

We applied this workflow to compare RPPs from published studies across Europe with our LSM-RPPs, implementing the following steps:

• 1.

  Debiasing RPPs for heavy pollen grain taxa: we adopted the adjusted RPP for Abies according to Theuerkauf (2025), as the original estimate was based on the Gaussian Plume model.

• 2.

  Pairwise ratio comparison: we calculated ratios between pairs of taxa in our dataset and compared them to corresponding ratios in previous studies. A taxon pair from a previous study was considered matching if its ratio differed by no more than a factor of 1.5 (i.e., within the range of 0.67 to 1.5) from the equivalent ratio in our dataset.

• 3.

  Selecting the scaling factor taxa: we selected taxa whose pairwise ratios passed the similarity window of 1.5 and used them to compute the scaling factor. In this study, RPPs from Erica and Genista-type were compared to RPPs of Ericaceae and Fabaceae, respectively, since these represent equivalent taxonomic groups identified at different resolutions across studies.

• 4.

  Compute scaling factor and adjust RPP values. To remove the influence of the reference taxon, we calculated a scaling factor for each previously published study. Specifically, we determined the mean RPP of the scaling taxa within our dataset and compared it to the mean RPP of the same taxa in each previous study. The scaling factor was obtained as the ratio of these two means (our dataset divided by the corresponding study). All RPP values from each previous study were then multiplied by this scaling factor, assuming that the offset between studies is systematic and affects all taxa proportionally.

• 5.

  Visualisation: for selected studies, we plotted the RPPs of comparable taxa using two approaches:

  • a.

    Spatial visualisation: we produced a map comparing the similarity of our RPPs with other studies in Europe.

  • b.

    Fold-change symmetric bar plots, where values <1 were plotted below the horizontal axis as reciprocals (e.g., a RPP of 0.1 is shown as $\mathrm{1}/\mathrm{10}$ below the axis), and values >1 were plotted above. This approach allows for multiplicative inverse of RPPs.

Previous RPP studies often include multiple sets of pollen samples, analysed using different methods. For each original study conducted in Europe, we selected one representative set of RPP values corresponding to a distinct set of pollen samples. Regarding the calculation methods (e.g., the ERV models), we followed choices outlined by Githumbi et al. (2022). However, the RPPs themselves were sourced from Wieczorek and Herzschuh (2020), who provided the original, unmodified estimates without recalculation. References of the studies used for comparison are available in Table 1.

3 Results

3.1 Pollen productivity estimates

The first RPPs relative to Poaceae in STI identify low, medium and high pollen producers (Fig. 5). Low pollen producers (lower than 0.48) are Asteraceae sf. Asteroideae, Brassicaceae, Genista-type, Juniperus, Plantago, Rosaceae and Rumex; medium producers (0.49–1.42) are Betula, Asteraceae sf. Cichorioideae, Corylus, Erica, Fagus and Ranunculaceae, while high producers (>1.43) are Abies, Artemisia, Amaranthaceae, Olea, Pinus, evergreen and deciduous Quercus.

Figure 5 Relative pollen productivity estimates (RPPs) relative to Poaceae calculated using the LSM, with standard deviation (SD) and standard error (SE) for 21 characteristic taxa in continental Spain. Red bars refer to the SD. Blue dashed lines indicate the cutoff values, based on percentiles, between low, medium and high pollen producers. Colours refer to the different groups of taxa (Poaceae in grey, arboreal in green, shrubs in teal and herbs in yellow). GPM-RPPs are shown in Fig. A1. The calculated RPPs for continental Spain are publicly available at https://doi.org/10.5281/zenodo.17927544 (Jungkeit-Milla et al., 2025).

3.2 Validation of RPPs on modern coretops

We applied the arboreal RPPs (relative to Poaceae) to the 27 coretops to test their validity (Fig. 6). We compiled all available samples containing Abies pollen that had not been used to retrieve the RPPs, selecting only 4 samples in total. The results of the validation for Abies are shown in Fig. B3 in Appendix B.

The REVEALS model corrects the uneven pollen productivity across taxa, resulting in an abundance reduction of overrepresented taxa and an increase for the underrepresented taxa (Fig. 6). This correction is evident when comparing raw pollen percentages (Fig. 6, left column) with REVEALS-based vegetation estimates (Fig. 6, right column). Pinus and both evergreen and deciduous Quercus, which dominate Iberian tree communities, show a marked reduction in estimated vegetation cover. The case of Pinus is especially illustrative: the adjusted estimates align closely with the observed vegetation cover (Fig. 6), suggesting the model performs well in accounting for its high pollen productivity. Another important example is Olea, for which REVEALS estimates indicate a reduction in the extent of olive cultivation. Nonetheless, several coretops were found to have no crops within a 45 km radius and we found that pollen from Olea correlated better with olive crop coverage at 100 km distance (see Figs. S4, S5 and S6 in the Supplement).

Figure 6 Relationship between the mean tree cover at 45 km in each coretop with untransformed pollen abundances (left, in red) and with the REVEALS-based vegetation estimates (right, in blue) with standard errors. The grey diagonal line represents the 1:1 reference line, indicating perfect correspondence between pollen/vegetation estimates and actual tree cover. Coloured lines represent linear regressions with 95 % confidence intervals (shaded areas). Points above the reference line suggest pollen overrepresentation relative to actual tree cover, while points below indicate underrepresentation.

On the contrary, temperate forest taxa such as Betula, Corylus and Fagus, having relatively low RPPs (Fig. 5), tend to be underrepresented in pollen assemblages in Iberia. REVEALS corrects for this bias, resulting in higher vegetation estimates for these taxa. Juniperus, moderately underrepresented in raw pollen data likely due to low pollen productivity (Fig. 5), is partially corrected by REVEALS, resulting in slightly higher estimates of vegetation cover. The standardised composite improvement score of the multimetric analysis at 45 km revealed that temperate forest taxa perform better with raw pollen counts than with the REVEALS-based reconstructions when comparing with observed regional vegetation cover (Fig. 7). Nonetheless, REVEALS estimates still provide more ecological sense in accounting for the observed vegetation than pollen percentages for Betula at 15 and 100 km (Figs. S4 and S6 in the Supplement).

Figure 7 Multimetric analysis comparing performance of REVEALS at 45 km for the selected arboreal taxa. Blue bars refer to positive values (better performance of REVEALS); red bars refer to negative values (better performance of pollen percentages).

3.3 Comparison of RPPs across European studies

Comparison of RPPs demonstrated similar values, at least in one pair, in all or more than half of the previously published values for the following taxa: Amaranthaceae, Artemisia, Asteraceae sf. Asteroideae, Brassicaceae, Asteraceae sf. Cichorioideae, Corylus, Ericaceae (Erica), Fabaceae (Genista-type), Fagus, Juniperus, Plantago, Poaceae, evergreen and deciduous Quercus, Ranunculaceae and Rumex. In contrast, the remaining taxa – Abies, Betula, Pinus and Rosaceae – show agreement with only half or fewer of the European sites (Fig. S7a in Supplement).

The eight studies with the highest number of matching taxa are from Romania, Czech Republic-Slovakia, Norway, France, Switzerland and Sweden. The remaining studies, which show lower similarity between RPPs, come from Great Britain, Germany, Denmark, Poland, Estonia, and Finland (Fig. 8).

Figure 8 Comparison of RPP studies across Europe: teal gradient indicates the number of species appearing in matching pairs; the size of the circle indicates the number of comparable taxa per study. Numbers on the map refer to the studied regions (see Table 1). White lines indicate the 8 most similar studies. See Fig. S7 in Supplement for standard barplots with number of comparable taxa.

For several taxa – Quercus (deciduous), Asteraceae sf. Asteroideae, Corylus, Ericaceae, Plantago, Poaceae, and Pinus – the STI RPPs are near the midpoint of the observed range in similar studies (Fig. 9). Notably, Poaceae and Pinus form two distinct clusters across studies; the STI values are within the higher cluster (4 for Pinus, and 1 for Poaceae).

For the remaining taxa, STI RPPs are distributed towards the edges of the observed ranges (Fig. 9), though still close to some previously reported values. These include Amaranthaceae, Artemisia, Asteraceae sf. Cichorioideae, Fabaceae, Fagus, Juniperus, Ranunculaceae, and Rosaceae. In the case of Abies, our estimate is higher than from any other European studies, while for Betula and Rumex, the values are lower.

Only two studies have derived RPPs for Brassicaceae, Fabaceae and evergreen Quercus, while only the present study provides a RPP for Olea. Within this limited data, evergreen Quercus is consistently represented as a high pollen producer, whereas Brassicaceae appears as a low pollen producer (Fig. 9).

Figure 9 Fold-change symmetric barplot of recalculated RPPs from the eight most similar studies. Adjustments account for reference taxon differences and dispersal model correction for Abies. Values <1 plotted as reciprocals. Each bar represents one of the eight most similar studies (from study 1 to 8), coming from Southern Sweden (8), Swiss Plateau (7), Jura Mountains in Switzerland (6), Southern France (5), Central Czech Republic (4), Southern Norway (3), White Carpathians (2) and Southeastern Romania (1). Circles represent visualisation of different studies for RPP ≈ 1. See Table 1 for more details of each study and Fig. S8 in the Supplement for RPPs with confidence intervals.

4 Discussion

The results presented here are the first RPPs produced in Iberia. Our results show clear differences in pollen productivity among the 21 taxa analysed. Among the low producers, herbaceous and shrub taxa dominate, while the medium and high pollen producers include a mix of trees and shrubs. Validation using coretops indicated that RPPs of the dominant taxa in present-day landscapes are more accurately estimated than the less abundant ones. Comparison with previously published European RPPs shows that our estimates align closely with 8 different studies. Our discussion is thus framed around our three main objectives of estimating, evaluating and comparing the Iberian RPPs under the prism of the methodological challenges we found.

4.1 Validation and potential implementation of the first Iberian pollen productivities

The validation presented in this study is the first carried out in heterogeneous environments of the Western Mediterranean using vegetation survey data and disaggregated by taxa rather than by total tree cover. Validation of RPPs in modern samples remains uncommon but has been performed for the whole of Europe (Serge et al., 2023), Southern Sweden (Hellman et al., 2008), Australia (Mariani et al., 2017, 2022) and various Asian regions (Jiang et al., 2020; Wan et al., 2022; Xu et al., 2014). They generally found an excellent performance for large groups of taxa, although the fit between the modelled estimates and the vegetation data was not always as good when single taxa were considered (Hellman et al., 2008). Wan et al. (2022) presented the first evaluation of RPPs in a tropical region, suggesting that the REVEALS model performs well especially when applied at the landscape level rather than for individual species. In our dataset, we could not validate REVEALS for herb and shrub taxa because MFE only includes arboreal taxa. Nevertheless, we hypothesise that some herbaceous and shrub taxa may be overrepresented when applying the REVEALS model (as in Li et al., 2023; Marquer et al., 2020). This is because, in the Iberian context, these taxa are relatively low pollen producers, similar to certain genera such as Betula, Corylus and Fagus.

REVEALS generally performed well when reconstructing the coverage of arboreal taxa (Fig. 6), although only 5 of the 8 taxa performed better with REVEALS than with pollen abundances, according to the multimetric analysis performed at 45 km radius: Pinus, Juniperus, Olea, deciduous and evergreen Quercus (Fig. 7). Juniperus showed improved performance as the distance increased, specifically up to 45 km (see multimetric results at other distances in the Supplement Figs. S4–S6), likely because junipers occur in small copses or patchy forest stands (e.g., J. thurifera), so the greater the distance considered, the more individuals are found. In this work, and as far as we know, we calculated the first pollen productivity estimate for Olea. We found that Olea is a high pollen producer, which aligns with previous studies highlighting its high dispersal capability (Cañellas-Boltà et al., 2009; Fernández-Rodríguez et al., 2014) and its presence in modern samples where sometimes the nearest olive trees are more than 70 km away (Leunda et al., 2017). In our modern samples, most of the Olea pollen comes from olive crops in southern Spain. The MFE inventory rarely includes wild olive trees, since they do not form forests or small groves naturally. These two examples underscore the need of accounting for source area when interpreting pollen records, particularly for taxa frequently employed as indicators of environmental change. The contrasting distribution patterns and dispersal capacities of Juniperus and Olea demonstrate how spatial context and distance exert a substantial influence on pollen representation.

The standard errors of the 8 validated taxa were consistently lower than their estimates, providing a measure of precision in terms of reliability (Githumbi et al., 2022; Li et al., 2023). Nevertheless, the quality of REVEALS outputs is ultimately constrained by the quality of input datasets (EMPD2 and MFE). According to the theory, deviation of REVEALS estimates from observed vegetation may suggest unreliable RPP values. We however believe that such deviation could be explained by: (i) a poor performance of REVEALS for Betula, Corylus and Fagus, since they are not present in all the STI due to environmental heterogeneity; and (ii) inaccuracies in the vegetation dataset, especially for taxa that do not form large, continuous forests such as Juniperus, Betula or Corylus (Wan et al., 2022).

Regarding the validation of Abies (Fig. B3), only 4 coretops were available to test the performance of REVEALS, both with GPM and LSM. The GPM-REVEALS yielded estimates up to three times higher than the LSM-REVEALS. These findings are consistent with previous studies (Abraham et al., 2014; Theuerkauf et al., 2013, 2016), which showed that large basins strengthen the overestimation of this taxon. The weighting factor decreases their distance-weighted plant abundance, and therefore the pollen signal is counterbalanced with high pollen productivity, ca. 2 times higher than the LSM-RPP (Fig. A1). As larger basins strengthen the overestimation of Abies (Abraham et al., 2014), RPP estimates of Abies may vary considerably among studies with different basin sizes (Fig. B3). Moreover, according to the multimetric analyses (Figs. B2 and B4), these indicate that LSM-REVEALS produced more accurate and consistent estimates of Abies cover than GPM-REVEALS at all distances. A similar pattern was also observed for Juniperus, for which LSM-REVEALS also outperformed the GPM approach. This is particularly important because Juniperus is a key taxon, containing 5 woody species in Iberia and thoroughly present in several fossil records across the Iberian Peninsula. This is the case with the longest continuous records in the western Mediterranean from the Villarquemado paleolake (Aranbarri et al., 2014; González-Sampériz et al., 2020; Moreno et al., 2012), which shows substantial proportions of juniper during the previous interglacial (the Eemian). An accurate reconstruction of the history of Juniperus is therefore essential for assessing ecosystem resilience under extreme climate fluctuations.

In light of our results, we believe that implementing our RPPs to quantitatively reconstruct vegetation cover from fossil pollen records could be promising for both arboreal and non-arboreal taxa in Iberia. The relatively low standard errors for trees in our validated coretops (especially for Pinus, Olea and both evergreen and deciduous Quercus) indicate that they can be confidently used for regional reconstructions of past vegetation. Even for some with lower performance, such as temperate forest species, these RPP values can still provide a useful reference until region-specific estimates are available. Both GPM- and LSM-RPPs proved that they can be used for reconstructions of vegetation dynamics. However, due to the aforementioned reasons, the LSM proved to provide more robust results, as the GPM may underestimate the dispersal of heavy pollen taxa, potentially leading to unrealistic vegetation reconstructions.

As for shrub and herb taxa, while direct validation was not possible, the generally robust SEs we retrieved suggest that these values could be cautiously implemented in reconstructions of open vegetation dynamics within Mediterranean landscapes. The overrepresentation of trees in pollen records aligns with findings from other European studies and with long-lasting intuitive expectations in the palaeoecological community. Notably significant are the values among dominant tree species: Pinus (4.00), evergreen (2.31), and deciduous (1.72) Quercus, suggesting a lower prevalence of pine communities in the past than has been interpreted based on previous reconstructions that indicated a higher abundance of pines, using both pollen data and/or wood charcoal and macrofossil analyses (Aranbarri et al., 2014, 2020; Carrión et al., 2004; Ezquerra et al., 2019; Múgica et al., 2001; Rubiales et al., 2010). Conversely, some deciduous trees in Iberia, such as Betula (0.66) and Fagus (0.52), showed lower productivity, potentially leading to them being underrepresented in fossil pollen records. This underestimation could affect our understanding of Pleistocene and Holocene temperate forest dynamics unless RPPs are applied (Broström et al., 2005; Githumbi et al., 2022; Trondman et al., 2015). For instance Corylus, with a RPP close to 1, provides crucial information for linking it to, e.g., early postglacial vegetation expansions (Aranbarri et al., 2014; González-Sampériz et al., 2006; Pearce et al., 2025b; Theuerkauf et al., 2014), and post-disturbance responses (Gil-Romera et al., 2014; Leunda et al., 2020).

For herb and shrub taxa with particularly low RPPs, such as Rumex (0.23), Brassicaceae (0.28), Juniperus (0.45) or Plantago (0.48), their roles in vegetation dynamics are expected to be underestimated if only raw pollen counts in fossil records are considered. This is especially critical for taxa associated with anthropogenic activities, such as Plantago or Rumex, whose historical presence in the Mediterranean landscape may be stronger than traditionally inferred without REVEALS estimates (Grindean et al., 2019; Kuneš et al., 2019; Mazier et al., 2008; Soepboer et al., 2007).

Our RPPs demonstrate that percentage pollen assemblages may bias reconstructions of Pleistocene and Holocene vegetation by overemphasising high pollen producers such as conifers and some oaks while downplaying the ecological importance of low-productivity taxa, including several deciduous broadleaved trees and herbs. Future quantitative reconstructions of vegetation dynamics in the STI hold the potential to substantially advance debates on the contentious idea of past continuous forest canopies in Iberia (Gomes et al., 2020; Pérez-Obiol et al., 2011). Incorporating revised pollen productivity estimates may reveal an even more fragmented forest landscape in the Southwestern Mediterranean, characterised by a mix of broadleaf woodlands, coniferous patches, temperate forests and open areas, as suggested in various studies (i.e., see compilations and references in Carrión et al., 2010a and González-Sampériz et al., 2010).

Future implementations of our RPPs in Iberia expand the possibilities for reconstructing past disturbances. By integrating RPPs with fossil pollen records, we can quantify vegetation composition over time and estimate biomass dynamics affected by various disturbances. With reliable RPPs, it would become feasible to produce spatially explicit models showcasing the long-term effects of environmental drivers across time and space, including anthropogenic or natural factors as well as the long-term interactions between disturbances and ecosystem recovery (Githumbi et al., 2022; Knight et al., 2022; Theuerkauf and Couwenberg, 2017).

4.2 Pollen productivities across Europe

Dominant vegetation types of the STI – pines, evergreen and deciduous oaks – exhibited RPP values successfully validated with independent top core samples. These taxa are identified as high pollen producers, consistent with findings from previous studies, particularly after adjusting original RPP values by the average scaling factor (see Sect. 2.3 “Comparison of RPPs across European studies”). Specifically, Pinus showed a RPP of 4.0 (Fig. 9), which aligns reasonably well with values reported from Central Czech Republic (4.5) and Southern Norway (4.8) (numbers 4 and 3 in Fig. 9, respectively). Evergreen Quercus exhibited a RPP of 2.31, comparable to 3.8 reported in Southern France (number 5, Fig. 9), and deciduous Quercus had a RPP of 1.72, which closely matches the value of 2.0 from the White Carpathians (number 2), in Czech Republic and Slovakia.

We also observed that our RPPs for most taxa fall near the midpoint of the range observed in other European studies, which is of interest since the STI has different environmental conditions that could have influenced pollen productivity. High pollen producers include wind-pollinated shrubs and herbs such as Amaranthaceae (3.31) and Artemisia (1.53), which yielded similar values to those reported in the Central Czech Republic (3.1 and 2.0, respectively, number 4 in Fig. 9). Medium pollen producers comprise shrubs and herbaceous taxa, including Erica (1.18), comparable to the value from Southern France (1.5) (number 5), Asteraceae sf. Cichorioideae (1.06) and Ranunculaceae (1.45), the latter two aligning with values of 1.0 and 1.3, respectively, reported from the White Carpathians (number 2 in Fig. 9). Low pollen producers – primarily entomophilous herbs or shrubs – such as Rosaceae (0.31) and Genista-type/Fabaceae (0.42) correspond well with previously reported values of 0.3 and 0.5 from Southeastern Romania (number 1).

In contrast, other taxa did not clearly match values reported in earlier studies, either due to a limited number of existing estimates or a wide variability in published data. For instance, RPP values for Asteraceae sf. Asteroideae (0.31) and Brassicaceae (0.28) are among the first reported or there are only a few previous estimates, making direct comparison difficult. For Corylus (0.94) our result lies within the mid-range of previously published values, which vary considerably (Fig. 9). Nonetheless, these values appear consistent with biological expectations – entomophilous herbs typically show lower pollen productivity than anemophilous trees – suggesting that the estimates are reasonable despite limited comparative data. For some taxa, Iberian RPPs follow a different trend than elsewhere in Europe, especially for Abies or Betula. When values are rescaled to remove the effect of the reference taxon, Abies emerged as a high pollen producer in STI (2.93), and Betula (0.66) presents lower productivity than in other European studies (Fig. 9). As discussed above, this underrepresentation of the temperate forest taxa may point to the need for computation of specific RPPs for the Eurosiberian biogeographical region, which is where these forest taxa are mainly found in STI.

RPPs for Olea (2.60) in this study are the first estimates in Europe. Both Olea and evergreen Quercus are only present in Southern Europe, and therefore new RPPs in all other Mediterranean peninsulas are needed in future studies to unravel the palaeoecological history of these species during the Holocene. It is of critical importance to better understand the history of key pollen taxa, especially when they vary their functional role in the ecosystem and become closely linked to human presence, as with Olea, offering new insights into anthropogenic triggers in the past. Indeed, Plantago and Rumex, also key taxa related to human activities, were found to present lower productivities in STI than in the eight most comparable studies from elsewhere in Europe. We argue that two main factors explain lower productivity of these taxa in Spain: first, Mediterranean climatic conditions may constrain the growth of these herbaceous taxa, thereby reducing their reproductive performance (López-Orozco et al., 2023); and second, traditional Iberian land-use systems (e.g., extensive grazing, dehesa-like agro-silvo pastoral mosaics) tend to maintain semi-open or wooded pastures rather than the continuously open, nutrient-rich grasslands common in Northern Europe (Connor et al., 2019).

4.3 Methodological constraints

Reconstructing past vegetation in Iberia with REVEALS poses unique methodological constraints owing to the region's pronounced topographic and environmental heterogeneity. This reason directly challenges one of the core assumptions of REVEALS: the need for spatially consistent pollen-vegetation relationships. Iberia's spatial heterogeneity often leads to an uneven vegetation distribution. This is particularly evident for mixed temperate forest taxa (e.g., Betula, Corylus, Fagus) which are confined to the northern fringe and are therefore not present across the entire STI. In addition, the lack of large lake basins, otherwise common in northern Europe, often forces reliance on smaller sedimentary basins, increasing uncertainty (Sugita, 2007a). Averaging RPPs across the STI may therefore underestimate productivity in optimal habitats.

Furthermore, the surface pollen collection may have introduced spatial biases in two different ways: (1) overrepresentation of certain vegetation communities and (2) overrepresentation of certain taxa. Both are related to the use of moss polsters as pollen traps, which are typically located in humid, mountainous environments (Fig. 1), especially in the Mediterranean region, limiting the representation of open vegetation landscapes from Iberia. In addition, moss polsters are often located under the tree canopy, especially in the Mediterranean region, perhaps resulting in overrepresentation of arboreal taxa. For instance, this could have been the case with Abies (see pollen percentage per grid in Fig. S9 and REVEALS estimates per grid in Fig. S10 in the Supplement), where both species (A. pinsapo and A. alba) coexist with the moss polsters in humid, cooler regions, potentially leading to higher estimates of pollen productivity.

Vegetation survey data also introduce uncertainty. The MFE maps well the dominant tree species but underrepresents taxa forming scattered stands (J. oxycedrus, J. communis), beyond the scope of forestry interest. The SIVIM database contains almost 150 000 phytosociological plant inventories. The frequency of the taxa is estimated on a scale of 1 to 5, which we transformed into percentages of abundances (equivalence table can be found in Table S4 in the Supplement). Despite the non-continuous nature of the SIVIM data, RPPs for shrub and herb taxa yielded robust and ecologically meaningful results, except for Genista-type, whose productivity estimate was lower than the standard deviation. Regarding the main source of pollen data, the EMPD2, some of the samples were incomplete or misclassified (e.g., some salt lakes were not specified as such; some samples had inaccurate coordinates or lacked information about the sample type or the sample context).

Moreover, the choice of a loss function could influence RPP estimates. In our work, we chose a weighted sum of squared errors (WSSE) method that calculates the optimal set of taxon-specific RPP values that minimises the discrepancy between observed regional vegetation composition and vegetation proportions reconstructed by REVEALS. This loss function addresses two main challenges in pollen-vegetation modelling: heteroskedasticity and taxon-specific biases. Pollen-vegetation relationships are inherently heteroskedastic, which means that the variance of the residuals scales with vegetation abundance (Sugita et al., 2010). Dominant taxa typically exhibit high pollen production, leading to smaller relative errors in their vegetation estimates, while rare or underrepresented taxa often show disproportionate noise due to low pollen counts and localised distributions (Broström et al., 2008). In order to partially account for these uncertainties, we weight errors inversely, as the WSSE function ensures that deviations for rare taxa contribute meaningfully to the optimisation, preventing their signals from being overshadowed by dominant taxa. By doing this, we acknowledge that residual uncertainty remains higher for the less frequent taxa.

Finally, a limitation of the inverted REVEALS approach is that the optimisation procedure seeks RPPs that maximise the agreement between REVEALS-derived vegetation estimates and observed vegetation cover. Consequently, part of the observed correspondence is an expected outcome of the optimisation process rather than a fully independent validation. We acknowledge that this limitation exists, and therefore a cross-validation approach was initially implemented to reduce this potential circularity, by explicitly excluding a subset of sites during RPP estimation and subsequently using them for validation, as described in Sect. 2.2.2 “Validation of RPPs”.

5 Conclusions

We used inverted REVEALS to produce the first pollen productivity estimates (relative to Poaceae) for the Spanish Territory of Iberia. Overall, we found that conifers (Pinus and Abies), both evergreen and deciduous Quercus, and Olea are high pollen producers in continental Spain. Temperate forest taxa (Betula, Corylus and Fagus) were identified as medium pollen producers, while shrub and herb taxa generally yielded lower RPPs, except for anemophilous taxa like Amaranthaceae.

We also performed the first validation of RPPs for arboreal taxa in southern Europe, using 27 present-day coretops and forest inventory data. The most frequent arboreal taxa in present-day landscapes performed better with REVEALS-based estimates than with raw pollen counts. Additionally, we conducted a bias-free comparison of our RPPs with other European datasets, finding similar values overall, except for temperate taxa and Abies. Future studies should examine whether more accurate estimates could be achieved by producing separate RPP datasets for the Eurosiberian and Mediterranean bioclimatic regions. Region-specific RPPs, tailored to bioclimatic variability, could improve the accuracy of vegetation reconstructions and disturbance assessments across the Iberian Peninsula.

Finally, these findings suggest that previous reconstructions of past vegetation dynamics in the Iberian Peninsula may have overestimated the presence of pines and oaks, and therefore fossil records from the Southwestern Mediterranean may require reinterpretation. Future research using the new relative pollen productivity estimates for Iberian taxa to generate quantitative vegetation reconstructions could indicate a more mosaic-like pattern of broadleaf woodlands, conifers, temperate forests and open ecosystems, aligning with most recent findings.

Appendix A: Estimation of GPM-RPPs

The Gaussian Plume model (GPM) is the most established model to simulate pollen transport and dispersal. Consequently, we also computed RPPs using the GPM. The results are shown in Fig. A1, complementing the LSM-RPPs in Fig. 5.

Generally speaking, GPM-RPPs showed similar results to the LSM-RPPs. However, the standard deviations are higher (marked in red in Fig. A1) on the GPM estimates. Significant differences are visible in the GPM-RPP values of Abies and Erica, which are ca. 2 times higher than the LSM-RPP values, related to their high fall speed and low wind speed (Abraham et al., 2014). The dispersal of both species is underestimated by GPM, which results in lower species-specific distance-weighted plant abundance. Thus, the pollen signal is counterbalanced with high RPPs.

Figure A1 Relative pollen productivity estimates (RPPs) relative to Poaceae and calculated using the GPM, with standard deviation (SD) and standard error (SE) for 21 characteristic taxa in continental Spain. Red bars refer to the SD. Blue dashed lines indicate the cutoff values, based on percentiles, between low, medium and high pollen producers. Colours refer to the different groups of taxa (Poaceae in grey, arboreal in green, shrubs in teal and herbs in yellow).

Appendix B: Validation of the GPM-RPPs

GPM-RPPs were validated in 27 coretops using the GPM for unstable conditions in the REVEALSinR package (Theuerkauf et al., 2016). Vegetation rings of 15, 30, 45 and 100 km have been calculated. Here we select the 45 km ring, which showed the strongest performance according to the multimetric analysis. Figure B1 shows the relationship between pollen abundance and tree cover (left column) and between the derived estimates (REVEALS-based vegetation estimates) and tree cover (right column) at a ring of 45 km distance.

Results from the GPM-REVEALS estimates showed a slight overrepresentation of Juniperus compared with the LSM-REVEALS estimates (Fig. 6), suggesting that the GPM is less suitable than the LSM for reconstructing the arboreal cover of this particular taxon. This result is also evident in Fig. B2, where the multimetric analysis was performed to evaluate the performance of the GPM estimates.

Figure B1 Relationship between the mean tree cover at 45 km in each coretop with untransformed pollen abundances (left, in red) and with the GPM REVEALS-based vegetation estimates (right, in blue) with standard errors. The grey diagonal line represents the 1:1 reference line, indicating perfect correspondence between pollen/vegetation estimates and actual tree cover. Coloured lines represent linear regressions with 95 % confidence intervals (shaded areas). Points above the reference line suggest pollen overrepresentation relative to actual tree cover, while points below indicate underrepresentation.

As previously reported, Abies is a well-known taxon posing challenges for the Gaussian Plume models. To further assess its performance, we selected all available coretops containing Abies pollen and applied both GPM- and LSM-derived RPPs to generate REVEALS estimates (Fig. B3).

Although only 4 samples in the validation dataset contained Abies pollen, the results show that the GPM can overestimate Abies abundance by up to a factor of 3. Multimetric analysis of Abies performed at all distances also suggests that the use of REVEALS (both with GPM and with LSM) is preferable to only using pollen counts (Fig. B4), except for the GPM at 100 km. Despite the good performance of REVEALS, LSM estimates yielded higher scores than the GPM-REVEALS estimates at all distances.

Figure B2 Multimetric analysis comparing performance of GPM-REVEALS at 45 km for the selected arboreal taxa. Blue bars refer to positive values (better performance of GPM-REVEALS); red bars refer to negative values (better performance of pollen percentages).

Figure B3 Relationship between tree cover (%) at 45 km and pollen abundance (left) and with the REVEALS estimates (right) of GPM-RPPs (up) and LSM-RPPs (bottom) for Abies.

Figure B4 Multimetric analysis comparing performance of Abies at different distances, using GPM- and LSM-REVEALS estimates. Blue bars refer to positive values (better performance of GPM-REVEALS); red bars refer to negative values (better performance of pollen percentages).