Pollen and plant data were collected from 34 calcareous fens in Estonia. Of the 34 sites, seven open and seven forested sites are located in southern Estonia (South region), with ten open and eight forested sites located in Saaremaa Island (West region) (Figure 1). The sites distributed between the two regions were located at least 2.5 km apart but the open and forested fens were sampled relatively close to each other (60 to 300 m apart). The terrain in the South region is predominantly dome-shaped consisting of undulating hillocky dead-ice landforms covered with forest and grasslands. Between hillocks, there are wet meadows, lakes and patches of mires (transitional mires, bogs, fens, alder fens). About 70% of the region is covered with boreonemoral forest dominated by Pinus sylvestris, Picea abies, Betula pubescens, and Alnus glutinosa. Saaremaa is the largest island in Estonia, located in the Baltic Sea. The island is composed of Silurian carbonate rocks, primarily limestone, that are covered by glacial sediments reworked by fluvial geomorphology during the development of the Baltic Sea (Saarse et al., 2009). Over 40% of the island is covered with forests, predominantly pine forests (Pinus sylvestris), mixed forests and broad-leaved (deciduous) trees (e.g., Quercus robur) are found. Sandy soils (Podzols), peatland soils (Histosols), waterlogged soils (Gleysols) and lessive soils (Luvisols) comprise most of the soil cover in southern Estonia, whereas Limestone rendzinas (Rendzic Leptosols) and calcareous soils (Calcaric Regosols) are typical in western Estonia (Reintam et al., 2005). Both study areas are regions with the densest calcareous spring fen network in Estonia. The fens in West region in Saaremaa, have relatively uniform vegetation, dominated by Schoenus ferruginaeus, Carex davalliana and Molinia caerulea, while the fens in South region are more variable in species composition but usually dominated by Carex spp., Molinia caerulea, Phragmites australis.

Moss polster samples (size of 5 cm in radius) were collected to determine the modern pollen assemblage from each of the fen sites (n = 34). The moss samples consisted of a relatively wide range of species, most abundant being Calliergonella cuspidata, Plagiomnium ellipticum, Scorpidium cossonii, Campylium stellatum, Sphagnum subnitens. In the majority of sites, the mosses did not form tense tussocks but the structure was relatively loose. For pollen analysis, only the green (living) upper part of the mosses was collected. Moss sample collection was synchronized with vegetation inventories. All the sampling and preservation of moss polsters followed the Crackles Bequest Project protocol (Bunting et al., 2013). For diversity reconstructions, we used sediment pollen assemblages from 53 sedimentary samples of Kanna spring fen (Blaus et al., 2019).

Both modern (polster) and historic (sediment) pollen samples were treated with HCl and 10% KOH followed by standard acetolysis method (Berglund and Ralska-Jasiewiczowa, 1986; Faegri and Iversen, 1989). Both moss and sediment samples were examined under a light microscope at magnifications of 250, 400, and 1000x. Approximately 1000 terrestrial and aquatic pollen grains per moss polster and sediment pollen sample were counted (min = 920; max = 1115), spores of vascular plants were also counted.

The potential of pollen data to reconstruct vegetation proportions (and vegetation diversity) is hampered by interspecific differences in pollen production and dispersal causing over-representation of some taxa and under-representation of other taxa in pollen assemblages (Birks et al., 2016a). To reduce the influence of high pollen producers, we used Andersen’s correction factors (Andersen, 1970) to downweight the dominant pollen taxa – a method that has been shown to improve the pollen-plant diversity relationships in earlier studies (Felde et al., 2016; Reitalu et al., 2019). All the pollen-based diversity estimates were, therefore, calculated on Andersen-transformed pollen counts.

Pollen data was divided into woody taxa (including trees and shrubs) and herbaceous taxa (including grasses and dwarf shrubs from the Ericaceae family) prior to the diversity calculations. Spores of pteridophytes are included in diversity calculations but for simplicity, we refer to the dataset as “pollen data” in the article.

Vegetation survey around each moss sample was carried out at the end of the flowering season during the second half of July and August (2017 and 2018). Survey timing and vegetation recording methodology followed Crackles Bequest Project protocol (Bunting et al., 2013), which is designed to produce vegetation abundance estimates in different distance classes in concentric areas around the pollen sample. Although vegetation mapping was conducted relatively late, only a few spring-flowering ephemerals (e.g., Anemone) might have been overlooked. At each site, vegetation was recorded within a 100 m radius around the moss sample (Figure 1). Within a 10 m radius, the vegetation was described in detail in five concentric areas at radii of 0.5, 1.5, 3, 6, and 10 m. In each 10 m area, twenty-one 1 × 1 m quadrats were systematically placed: one quadrat on top of the moss sample, four quadrats positioned at each cardinal direction at 1, 2.25, and 4.5 m from the central moss polster, and eight quadrats at 8 m distance from the center (Figure 1). In each 1 × 1 m quadrat, the percentage cover of all vascular plants was recorded. Additionally, species not occurring in the quadrats but within the circles were recorded. Between 10 and 100 m, vegetation types were mapped by using orthophoto maps (incl. infrared orthophoto maps) and field-work observations. The vegetation types were delimited keeping in mind the pollen perspective – for example, single trees and shrubs in open fen areas were mapped as separate entities because their role in pollen signal is potentially high (Bunting et al., 2013). Species composition of each vegetation type in 100 m radius was characterized in the field. For woody taxa, the percentage cover of each species was estimated. For herbaceous taxa, the Braun-Blanquet cover-abundance scale was used (Braun-Blanquet, 1964). The abundances were then translated to cover estimates within each vegetation type as follows: ±0.01%; 1 – 5%; 2 – 10%; 3 – 25%; 4 – 50%; 5 – 75%.

Altogether, plant cover estimates were calculated for nine different radii; 0.5, 1.5, 3, 6, 10, 25, 50, 75, and 100 m. For the radiuses up to 10 m, the plant cover estimates were averaged from the quadrats belonging to the particular radius. For the radiuses 25, 50, 75, and 100 m, the cover of each species was calculated from species cover estimates in different vegetation types weighted by the area of the vegetation type within the radius (Figure 1).

Similarly to pollen data, plant data was divided into woody taxa (including trees and shrubs) and herbaceous taxa (including grasses and dwarf shrubs from the Ericaceae family) prior to the diversity calculations.

In order to clarify which spatial scale of vegetation is best reflected by modern pollen, we used plant cover estimates of nine different radii (0.5, 1.5, 3, 6, 10, 25, 50, 75, and 100 m) and plant richness of woody and herbaceous species was calculated for each radius. In addition to richness estimates, we used Principal Component Analysis (PCA) with Hellinger-transformed plant cover data (Legendre and Gallagher, 2001) to characterize the vegetation composition. PCA was performed separately for woody and herbaceous plant species in each radius. PCA ordinations of pollen data and plant data for two radii (6 m for herbaceous species and 100 m for woody species) are given in Supplementary Figure S2. Pollen richness estimates and PCA ordinations of Hellinger-transformed woody and herbaceous pollen percentages were used for comparison with corresponding plant estimates in different radii. Pollen-plant richness relationships in different radii were studied with the help of Pearson’s correlation coefficients. Pollen–plant compositional relationships were analyzed with the help of Procrustes analysis (Gower, 2001) that introduces uniform scaling and rotation to minimize the squared differences between pollen-based and vegetation-based ordinations. The Procrustes symmetric sum of squares was used as “Procrustes correlation” to quantify the pollen–plant compositional relationship. We used the bootstrap method to provide confidence intervals for Pearson’s correlation coefficients and Procrustes correlations (Davison and Hinkley, 1997).

The R environment (version 3.4.4) (R Core Team, 2018) was used for the analyses with packages “vegan” (Oksanen et al., 2017) and “boot” (Canty and Ripley, 2019) for ordination analyses and bootstrapping, respectively.

To calculate the functional trait composition (CWM) and functional alpha diversity (FDa), we used nine traits that reflect different aspects of plant life history: clonality, height, specific leaf area (SLA), leaf dry matter content (LDMC), seed size, seed number, mycorrhizal type. For mycorrhizal type, we used the affiliation of plant species in three plant mycorrhizal types: arbuscular mycorrhiza (AM), ectomycorrhiza (EM) and ericoid mycorrhiza (ERM). Traits, the sources and availability of trait data and the functional significance of different traits are summarized in Table 1. Trait availability was the lowest for LDMC (only 40% of the species) and for SLA (62%), however, trait values were available for the most dominant taxa and we chose to include them in functional diversity calculations. Pteridophytes were included in trait assignments (except for seed weight and number).

Assigning trait values to pollen types that have only one species (for example, Pinus pollen type only includes P. sylvestris in Estonia) was straightforward. To calculate trait values for pollen types that included several different species the average trait values were calculated from the species recorded in our vegetation survey (for example, Potentilla pollen type includes P. anserina, P. erecta, P. reptans). Averages were also used for mycorrhizal types. In case one species can form two types of mycorrhizal associations, it is used in both estimates. Averaged trait values for pollen types are given in Supplementary Table S1.

Functional alpha diversity and CWM trait values were used to express FD. Species cover estimates and Andersen-transformed pollen percentages were used as weights for CWM calculations in plant and pollen data, respectively. “Functcomp” function in the “FD” R package (Laliberté et al., 2014) was used to calculate CWMs. FDa was calculated as an abundance weighted standardized effect size of mean pairwise distance (SESmpd) according to de Bello et al. (2016) for individual traits as well as all traits combined. Negative SESmpd indicates functional convergence and positive SESmpd indicates functional divergence; values ≤−2 and ≥2 indicate significant (P < 0.05) convergence and divergence, respectively (Gotelli and Graves, 1996). All diversity measures were calculated separately for herbaceous and woody taxa.

In addition to the functional diversity, community weighted EIVs of light (L), temperature (T), continentality (K), moisture (F), pH reaction (R) and nitrogen or nutrient availability (N) were tested with modern pollen-plant relationships. EIVs were calculated as community weighted mean values similarly to CWM trait values described above.

Phylogenetic diversity (PD) is based on the phylogenetic relationship amongst the coexisting species in the community and is represented by a phylogenetic tree (Brocchieri, 2016). To estimate PD, we calculated the phylogenetic tree from the latest megaphylogeny for vascular plants in the “V.PhyloMaker” R package (Jin and Qian, 2019). Similarly to FDa, PD was expressed as the standardized effect size of mean pairwise distance (SESmpd). The mean pairwise distance was estimated across the whole phylogenetic tree by averaging all species pairwise distances (Pavoine and Bonsall, 2010). SESmpd was calculated by comparing the observed community mpd to the null distribution of randomized communities with equal species richness using the null model “taxa.labels,” implemented in the “picante” R package (Kembel et al., 2010). SESmpd was calculated separately for woody and herbaceous species, and additionally for only angiosperms in both woody and herbaceous groups. Negative values of SESmpd indicate phylogenetic clustering (i.e., the coexistence of closely related species), whereas positive values indicate over-dispersion (i.e., the coexistence of distantly related species) (Webb et al., 2002). When pollen taxa included more than one species, we randomly chose one representative species from the list of species in our vegetation survey, for example, Primulaceae was represented by Primula farinosa.

Based on our results of richness and vegetation composition analyses in different radii (Figure 1) and published literature on the source area of pollen for the herbaceous (Hjelle, 1998; Bunting, 2003; Bunting and Hjelle, 2010) and woody species (Broström et al., 2005; Mazier et al., 2008), we used two radii – 6 m for herbaceous species and 100 m for woody species – to test the modern pollen–plant functional and phylogenetic diversity and EIV relationships. To test for significant pollen-plant associations of CWM, FDa, PD, and EIV measures we used correlation analyses. Prior to the correlation tests, we tested for the normal distribution of the variables and several variables were log-transformed to achieve normality (Tables 2−5). Pearson’s correlation coefficient was used when both variables were normally distributed. Alternatively, Kendall’s Tau coefficient was used. To control for the effect of multiple tests (52 correlation tests altogether), we used Benjamini and Hochberg (1995) correction approach.

To test the influence of the region (West and South) on the pollen–plant diversity relationships, linear regression models were fitted for all pollen–plant diversity associations (CWM, FDa, PD, and EIV) with pollen diversity (pollenD) as the response variable and corresponding plant diversity (plantD), region (Reg) and their interaction term (plantD:Reg) as explanatory variables. The full model was given as: pollenD ∼ plantD + plantD:Reg. Backward selection with the function “stepAIC” [and with k = 4 (the multiple of the number of degrees of freedom used for the penalty)] in the “MASS” R package (Venables and Ripley, 2002) was used to clarify which combination of predictors significantly explained pollen richness. Assumptions for using linear models were tested. In some cases, outlier(s) were excluded in order to meet the assumptions for linear modeling (Supplementary Figure S3).

Finally, we applied the different diversity metrics (richness, CWM, FDa, PD) and EIVs tested in modern pollen–plant data to Kanna spring fen pollen sequence (53 pollen samples). Kanna spring fen is located in Viidumäe Nature Reserve (22.096721 °E; 58.325031 °N) in the western part of Saaremaa Island (Figure 1). The peat accumulation in Kanna began 9.2 ka (calibrated kilo annum before present). The age-depth model is based on 11 AMS radiocarbon dates and dates from spherical fly ash particles, peat accumulation was generally uniform and continuous (Blaus et al., 2019). Blaus et al. (2019) separated three main stages of the Kanna fen development: bog phase (9.2–7.2 ka) characterized by ombrotrophic conditions; fen phase (7.2–0.4 ka) characterized by minerotrophic conditions; and “recent phase” characterized by shifts in dominant fen taxa, possibly caused by human impact. These three development phases were also used in the current study to examine the functional (FDa, CWM), PD and EIV shifts associated with major changes in Kanna fen development. We used the diversity metrics which showed significant correlations between modern pollen and plant data (Tables 2–5) for reconstructions. Generalized additive models (GAM) with the 95% confidence intervals were used to fit the relationships of each diversity variable with time using the “mgcv” R package (Wood, 2011).

The results of Pearson’s correlations with bootstrapped 95% confidence intervals showed significant positive correlations for herbaceous taxa at the scales of vegetation estimate up to 6 m (r = 0.42–0.47, Figure 2). For woody taxa, the pollen-plant richness correlation was not significant at any of the vegetation estimate scales. Compositional relationships assessed by Procrustes’ analysis showed positive associations for herbaceous taxa at all vegetation estimate scales. For woody taxa, the Procrustes correlations were significantly positive at the scales of 6 to 100 m.

Linear regression model showed that for herbaceous taxa both plant richness (t = 4.26, p < 0.001) and region (t = −5.02, p < 0.001) were significant predictors of pollen richness (determination coefficient R² = 0.5). Pollen richness increased with increasing plant richness in both regions and region South had significantly higher pollen richness than region West (Figure 3). In contrast, for woody taxa, the regression model did not indicate any significant variables describing pollen richness (Supplementary Figure S1b).

Correlations between pollen–plant CWM traits of woody taxa showed significant positive pollen–plant associations for SLA, seed number, AM, LDMC, clonality and plant height (Table 2). Linear regressions for CWMs of woody taxa indicated that both plant-based CWM estimate and region (South vs. West) were significant predictors of pollen CWM (Table 2 and Supplementary Figure S3). However, the interaction term of the region and plant CWM was not significant. For herbaceous taxa, pollen and plant CWM traits were significantly positively correlated for four traits: ERM, SLA, LDMC and seed number. For the herbaceous CWM of LDMC and seed number, the plant-based estimate was the only significant predictor of pollen-based CWM estimate (Table 2). For SLA, the pollen-plant CWM relationship differed significantly between the regions.

Pollen and plant functional alpha diversity (FDa) of woody taxa was significantly positively correlated for four traits: seed number, SLA, clonality and AM (Table 3). Linear models for FDa of woody taxa showed that for seed number and clonality, plant estimate was the only significant predictor of pollen diversity, for SLA and clonality, region was also significant (Table 3 and Supplementary Figure S3). The correlation between pollen and plant FDa calculated across all traits was significantly positive for woody taxa and the association did not depend on the region. Herbaceous taxa showed positive significant FDa correlation for three traits: ERM, clonality and SLA. Linear regressions showed that for ERM and clonality, corresponding plant-based diversity estimate was the only significant predictor, for SLA the pollen-plant relationship differed between the regions (Table 3).

Phylogenetic diversity (PD) was significantly positively correlated between pollen and plant data only for herbaceous taxa, the significant correlation disappeared when pterophytes were excluded from the calculation (Table 4). Linear regression indicated that the slope of herbaceous pollen-plant PD relationship differed between the regions (Table 4). Woody taxa did not show significant PD correlations.

The Ellenberg indicator values (EIVs) based on woody taxa were significantly positively correlated between pollen- and plant-based estimates for the indicators of light (L), temperature (T), soil reaction (R) and moisture (F) (Table 5). Linear regressions showed that for woody taxa pollen-based T estimate, plant-based T was the only significant predictor, pollen-based L and F were significantly associated with corresponding plant-based estimate and region (Table 5). Ellenberg R estimate based on woody pollen taxa was significantly associated only with the region. EIVs based on herbaceous taxa were significantly positively correlated for L, T and nutrient content (N). Linear regressions showed that for herbaceous pollen-based N estimate, plant-based N was the only significant predictor, for L and T region was also significant (Table 5 and Supplementary Figure S3).

Pollen richness of both woody and herbaceous taxa was the lowest in the beginning of the mire development (9-8 ka) and increased between 8 and 6 ka (Figure 4). Richness of herbaceous taxa was more or less stable until present day but richness of woody taxa decreased during the last 2 ka.

Community weighted mean values of AM, SLA and LDMC showed relatively similar trends for woody and herbaceous taxa through the fen development – AM and SLA had their maximum values between 6 and 2 ka and LDMC was the lowest during the same period (Figure 5). Plant height, seed number and clonality had contrasting trends through time for woody and herbaceous taxon groups. For woody taxa, clonality was the highest and height was the lowest between 7 and 2 ka. For herbaceous taxa, clonality was the lowest around 7 ka and the plant height had its maximum between 6 and 2 ka. Herbaceous seed number and ERM were the highest in the beginning of the study period (9-7 ka) (Figure 5).

Functional alpha diversity was the highest in the beginning of the study period for several traits: AM and clonality for both woody and herbaceous taxa, seed number and ERM for herbaceous taxa and SLA for woody taxa (Figure 6). During the middle of the study period, the FDa values of several traits decreased notably indicating a shift from functional divergence toward functional clustering. FDa values were close to or above 2 indicating significant functional divergence for seed number, AM and ERM in herbaceous taxa and for SLA and AM in woody taxa. The tendency toward functional convergence was detected only for woody taxon clonality.

Phylogenetic diversity showed relatively little change through time (Figure 7). The highest PD based on woody taxa was shown for the beginning of the study period (9-7 ka) when PD was close to or above 2 indicating phylogenetic divergence. When gymnosperms were excluded from the PD calculation, the values dropped drastically and indicated phylogenetic convergence throughout the study period (SESmpd < −2). PD of herbaceous angiosperms increased around 7 to 6 ka.

Ellenberg indicator valuess based on woody taxa indicated that moisture (F), light (L), and continentality (K) were higher during the beginning of the study period and close to the present-day and had the lowest values between 6-5 ka. For herbaceous taxa, the same indicators show little change or slight increase toward the present-day. The R for herbaceous taxa showed a strong increase from 7 to 5 ka, R for woody taxa increased during the same period but the change was less pronounced. Nutrient availability (N) increased according to both woody and herbaceous taxa between 7 and 5 ka but the N-value based on woody plants decreased and N-value based on herbaceous taxa stayed more-or-less stable from 5 ka until the present.

Pteridophytes and gymnosperms, being evolutionarily distant from angiosperms, have a strong impact on the phylogenetic structure of the vegetation (e.g., Massante et al., 2019). Modern pollen-plant PD associations in our study were also strongly influenced by pteridophytes (especially ferns) and gymnosperms. Pollen-plant PD correlation was significant for herbaceous taxa but only when pteridophytes were included in the calculation (Table 4). Ferns and their spores were closely related between vegetation and pollen samples, meaning that if ferns were present in vegetation their spores occurred frequently in the moss polster.

Results of PD reconstruction in Kanna fen exhibited significant clustering (phylogenetic convergence) of woody taxa throughout the 9.2 ka when gymnosperms were excluded from the PD calculation but overdispersion (phylogenetic divergence) was evident when gymnosperms were included in the reconstruction (Figure 7). This pattern is consistent with the study of Reitalu et al. (2015) who showed that gymnosperms significantly increased the phylogenetic overdispersion. Massante et al. (2019) found that woody taxa are phylogenetically overdispersed at high latitudes compared to low latitudes, because after originating in the tropics, specific and distant lineages of woody taxa, i.e. gymnosperms, were able to adapt and survive extinctions in cold temperate conditions. The slight decrease in overdispersion in woody taxa in the mid-Holocene is probably caused by climate changes that shaped PD by colonization of broad-leaved species (Eiserhardt et al., 2015). Weakening of the clustering of woody angiosperms at the same time supports the studies that identify overdispersion as related to warm and nutrient-rich conditions (Cavender-Bares et al., 2004; Spasojevic and Suding, 2012). These favorable conditions for plant growth most probably weaken the effect of environmental filtering for specific phylogenetic clades with adaptations for extreme environmental conditions, and allow species from distant lineages to coexist.

The PD of angiosperm herbaceous taxa followed the main phases of Kanna fen development (Figure 7). During the ombrotrophic stage (ca. 9-7 ka), the PD of herbaceous angiosperms tended toward phylogenetic clustering, most probably influenced by the high abundances of different Ericaceae species (Blaus et al., 2019). Clustering has been observed for species pools of relatively young habitats (Lososová et al., 2015), and could be related to site emergence from the sea and extreme environmental conditions. Clustering has commonly been explained as a result of environmental filtering, since closely related species are usually expected to be ecologically similar (Webb et al., 2002). In particular, low nitrogen and soil reaction 9.2 ka might have led to environmental stress and filtering the lineages not adapted to particular conditions (González-Caro et al., 2014). However, the functional diversity reconstructions showed functional over-dispersion for several traits during the same period, indicating that the functional-phylogenetic diversity relationships are more complex.

The fen phase (ca. 7-2 ka) was characterized by relatively high PD as a result of phylogenetically distinct species co-occurring together, and may depend mechanistically on climate as shown by Svenning et al. (2015). This stage coincides with the increase in summer temperatures (Renssen et al., 2009) and the maximum of palynological richness (Figure 4). Empirical studies have justified that PD typically covaries with species richness at different scales (Mace and Purvis, 2008; Mooers et al., 2008; Kluge and Kessler, 2011), reflecting patterns of species migration and diversification (Forest et al., 2007). PDs of both herbaceous and woody taxa showed a slight decrease toward the present day, likely in response to the increase in Cyperaceae within the fen, and decline of broad-leaved taxa in the surrounding landscape, respectively. The decline also coincides with decreasing pollen richness. In contrast, Reitalu et al. (2015) showed a decrease in PD in spite of the increase in pollen-derived richness associated with human impact, particularly, the facilitated increase in ruderal communities. To conclude, PD can be reconstructed from pollen data at least to some extent, it has a different dynamic in time among the plant growth forms and our example shows that period of high taxonomic richness in the natural fen system is also characterized by high PD.

The applicability of Ellenberg indicator values (EIVs) for different habitat types and in different regions in Europe has evoked considerable discussion in plant ecological literature (e.g., Schaffers and Sýkora, 2000; Diekmann, 2003; Smart and Scott, 2004; Williams et al., 2011; Bartelheimer and Poschlod, 2015). EIVs have been widely used in different mire systems to study a wide range of hypothesis (e.g., Cornwell and Grubb, 2003; Williams et al., 2011; Andersen et al., 2013; Gawenda-Kempczyńska, 2016; Klimkowska et al., 2019) and our findings add to those studies by showing the indicators that could be applied on sediment pollen and therefore could be used to interpret historical environmental gradients.

In current study, we found significant positive pollen-plant correlations for woody taxon EIVs of light (L), temperature (T), humidity (F), and soil reaction (R) (Table 5). Diekmann (2003) points out that short environmental gradients and within-sample heterogeneity may lead to unreliable EIV results. Because pollen samples of woody taxa are “sampling” vegetation from a relatively large area, results from heterogeneous study areas may be misleading as they can only provide average values for the entire area.

For herbaceous taxa significant pollen-plant associations were found for EIVs of light (L), temperature (T) and nutrient availability (N) (Table 5). Humidity (F) and soil reaction (R) that are recognized as main drivers of plant species diversity and composition of spring fen ecosystems (Hájková and Hájek, 2003; Sekulová et al., 2011) did not show significant pollen-plant associations. Our modern pollen-plant sampling scheme was designed to include open and overgrown fens but did not include large gradients of pH and humidity. The mean pH was 6.5 ± 0.7 and only three sites had pH < 6. Spring fens have relatively constant water availability, therefore moisture levels are high and do not vary much between the fens. The lack of positive pollen-plant associations for F and R may, therefore, be related to the relatively short gradient lengths in these indicators (cf. Diekmann, 2003). Ellenberg R has been shown to be relatively weakly associated with pH in wetlands (Williams et al., 2011) and might be better associated with calcium content (Schaffers and Sýkora, 2000). Our results showing significant pollen-plant correlations for L and N are promising as both indicators have been shown to be significantly associated with major gradients in fen species composition (Kotowski and van Diggelen, 2004; Andersen et al., 2013) and have great potential for pollen-based palaeoreconstructions from fens.

The reconstruction of EIVs from Kanna sediment pollen shows that the beginning of site development is characterized by low Ellenberg values of N and R (Figure 8) – which is in line with functional characteristics of both woody and herbaceous taxa and with earlier interpretations of site development history (Blaus et al., 2019). Indicators of L, F and K revealed by woody taxa were relatively low from 7 to 3 ka coinciding with mid-Holocene thermal maximum (Davis et al., 2003) agreeing with the studies reporting dryer climatic conditions during the mid-Holocene (Hammarlund, 2003; Seppä and Poska, 2004). The opposite trends in R and L in the mid-Holocene support the results showing negative R and L correlations in forest, mostly because acido-tolerant but light-demanding species are outshaded by closed forest canopy (Diekmann, 2003).

From the herbaceous EIV reconstructions, R and N showed an increase between 7 and 5 ka (Figure 8) corresponding to the period when minerotrophic conditions developed in the fen. Ellenberg R has been shown to be well correlated with fen calcium carbonates (Schaffers and Sýkora, 2000) while Ellenberg N reflects biomass or fertility (Hill et al., 2000; Schaffers and Sýkora, 2000; Diekmann, 2003). Our results agree well with Bartelheimer and Poschlod (2015) who showed that N is positively associated with SLA and plant height among herbaceous species. High N coincides with high pollen richness and with high phylogenetic diversity of herbaceous taxa.

Temperature reconstructions from pollen data based on modern pollen calibration datasets are widely used in Quaternary palaeoecology (e.g., Seppä and Poska, 2004; Salonen et al., 2012). Our results together with earlier studies (e.g., Kuneš et al., 2011; Reitalu et al., 2015; Felde and Birks, 2019) indicate that there is also potential for reconstructing long-term changes in soil reaction (calcium content), fertility, moisture and light.