Modern Pollen–Plant Diversity Relationships Inform Palaeoecological Reconstructions of Functional and Phylogenetic Diversity in Calcareous Fens

Blaus, Ansis; Reitalu, Triin; Gerhold, Pille; Hiiesalu, Inga; Massante, Jhonny Capichoni; Veski, Siim

doi:10.3389/fevo.2020.00207

ORIGINAL RESEARCH article

Front. Ecol. Evol., 03 July 2020

Sec. Paleoecology

Volume 8 - 2020 | https://doi.org/10.3389/fevo.2020.00207

This article is part of the Research Topic Functional Traits as Indicators of Past Environmental Changes View all 11 articles

Modern Pollen–Plant Diversity Relationships Inform Palaeoecological Reconstructions of Functional and Phylogenetic Diversity in Calcareous Fens

$\r\nAnsis Blaus*$ Ansis Blaus^1*

Triin Reitalu¹

Pille Gerhold^2,3

Inga Hiiesalu²

Jhonny Capichoni Massante²

Siim Veski¹

¹Department of Geology, Tallinn University of Technology, Tallinn, Estonia
²Institute of Ecology and Earth Sciences, University of Tartu, Tartu, Estonia
³Institute of Agricultural and Environmental Sciences, Estonian University of Life Sciences, Tartu, Estonia

Predicting the trajectory of ongoing diversity loss requires knowledge of historical development of community assemblages. Long-term data from palaeoecological investigations combined with key biodiversity measures in ecology such as taxonomic richness, functional diversity (FD), phylogenetic diversity (PD) and environmental factors expressed as Ellenberg indicator values (EIVs) could provide that knowledge. We explored the modern pollen–plant (moss polster pollen vs. surrounding vegetation) diversity relationships for herbaceous and woody taxa in calcareous fens from two different regions in Estonia, NE Europe. Associations of taxonomic richness, vegetation composition, FD (including functional alpha diversity and trait composition), PD and EIVs in modern pollen vs. plant data were studied with correlation analysis, Procrustes analysis and linear regression models. To test their potential use in palaeoreconstructions, diversity measures were applied on pollen data from Kanna spring fen reflecting fen vegetation development over the last nine millennia and diversity changes through time were studied using generalized additive models. Results showed significant pollen–plant richness correlations for herbaceous taxa at vegetation estimate scales up to 6 m radius and Procrustes analysis showed significant compositional associations at all plant estimate scales (up to 100 m). Woody taxa had no significant pollen–plant richness correlations but composition relationships were significant at plant estimate scales of 6–100 m. Traits that were best reflected by pollen data (both in terms of trait composition and functional alpha diversity) among woody and herbaceous taxa were seed number, clonality, SLA and LDMC. PD of herbaceous species was reflected by pollen data. Among the EIVs, Ellenberg L and T were significantly reflected by pollen data for both woody and herbaceous communities. Palaeoreconstruction from Kanna fen indicates that trends of woody taxa are mostly related to long-term changes in climate while diversity variables of herbaceous taxa closely follow autogenic processes within the fen. We suggest that pollen-based diversity estimates should be calculated separately for woody and herbaceous taxa as they clearly represent different spatial scales. Present study suggests that linking sedimentary pollen data with FD, PD and EIVs provides possibilities to examine long-term trends in community assembly and ecosystem processes that would be undetectable from traditional pollen diagrams.

Introduction

Identifying the drivers underlying biodiversity changes is among the key research questions both in ecology (Vellend et al., 2017) and palaeoecology (Birks et al., 2016b). Biodiversity can be measured at different organizational levels – from genetic diversity to taxonomic (species) diversity, to landscape diversity. Coupling these measures with community functional diversity (FD), which reflects the type, scope, and abundance of functional traits in communities provides a link between species diversity and ecosystem functioning (Tilman, 2001; Lavorel and Garnier, 2002; Díaz et al., 2007). In addition, phylogenetic diversity, which estimates the phylogenetic relatedness among the coexisting species (Faith, 1992) connects local community assembly and evolutionary processes (Gerhold et al., 2015, 2018). The concepts of functional and phylogenetic diversity are of central interest in a large number of contemporary ecological studies of land plants (e.g., Cadotte et al., 2009; Cavender-Bares et al., 2009; Purschke et al., 2013; Díaz et al., 2016). In palaeoecology, plant functional types have been used for several decades to enable the transformation of pollen records to biomes or land-cover types for comparison with climate models (e.g., Prentice et al., 1996; Fyfe et al., 2010; Davis et al., 2015). The more detailed studies of functional and phylogenetic aspects of palaeo communities are relatively rare (e.g., Lacourse, 2009; Reitalu et al., 2015; Brussel et al., 2018; Carvalho et al., 2019; Jabłońska et al., 2019) and have demonstrated the great potential of functional and phylogenetic palaeodiversity to contribute to a better understanding of processes underlying long-term patterns of community assembly. For example, Brussel et al. (2018) showed that long-term fire frequency variations may drive directional selection for fire–adapted plant community attributes.

In palaeoecology, sedimentary pollen analysis is one of the main tools for reconstructing changes in past vegetation and plant diversity in the late Quaternary (Berglund and Ralska-Jasiewiczowa, 1986; Birks and Birks, 2006; Gaillard et al., 2008; Reitalu et al., 2015). However, interpreting pollen data in terms of past plant diversity is not straightforward because of the low taxonomic resolution of pollen data and interspecific variation in pollen production and dispersal (Odgaard, 1999; Liu et al., 2014; Birks et al., 2016a; Carvalho et al., 2019). To clarify the pollen–plant relationships and to aid the pollen-based biodiversity reconstructions, modern pollen assemblages acquired from moss polsters, pollen traps, and lake surface sediments have been compared with detailed information of the surrounding vegetation (Bunting and Hjelle, 2010; Meltsov et al., 2011, 2012; Matthias et al., 2015; Felde et al., 2016; Julier et al., 2018; Reitalu et al., 2019). Several studies have demonstrated a positive association between modern pollen richness and plant richness (Meltsov et al., 2011; Matthias et al., 2015; Felde et al., 2016; Reitalu et al., 2019). Pollen of woody taxa is expected to be produced in greater quantities, to travel larger distances and to reflect broader spatial scales than the pollen of herbaceous and/or insect-pollinated taxa (Meltsov et al., 2011; Felde et al., 2016; Reitalu et al., 2019). In the present study, we will, therefore, examine the woody and herbaceous species/pollen taxa separately.

Carvalho et al. (2019) have highlighted the importance of using modern pollen–plant relationships also for functional palaeodiversity. They showed in modern pollen–plant study from fen peatland that some leaf traits (C and N content) were better reconstructed from pollen data than others (specific leaf area and leaf dry matter content). The choice of traits for analyzing functional diversity is often debated and depends on the research system, hypotheses and data availability (Lavorel and Garnier, 2002). In the current study, we have focused on the commonly used traits in plant functional ecology that could be assumed to be relevant for shedding light to the past plant communities. The chosen traits reflect different aspects of plant life history related to vegetative growth [clonality, plant height, specific leaf area (SLA), leaf dry matter content (LDMC)], reproduction (seed size and number) and biotic interactions (mycorrhizal type). In palaeoecological context, plant height, SLA, LDMC, seed size and mycorrhizal type can be assumed to be related to long-term climate variability and indicative of niche differentiation during favorable climatic conditions (Lacourse, 2009; Reitalu et al., 2015; Knight et al., 2020). Clonality, plant height, seed number and mycorrhizal type can be assumed to be associated with plant establishment in the harsh conditions of novel habitats (Ye et al., 2014; Pither et al., 2018) and with waterlogging in wetlands (Jabłońska et al., 2019).

There is a wide array of measures and several conceptual frameworks for estimating the functional and phylogenetic aspects of community assembly (see overviews in Vellend et al., 2011; Garnier et al., 2016). The functional structure of the community (functional diversity sensu lato, FD) can be characterized by average trait values and by measures estimating the variation in trait values (functional diversity sensu stricto or functional alpha diversity) (Garnier et al., 2016). Community weighted mean (CWM) trait values are widely used for evaluating links between community dynamics and ecosystem processes (Lavorel et al., 2008; Muscarella and Uriarte, 2016; Duarte et al., 2017). Functional alpha diversity (FDa) reflects the variability of traits among the coexisting species and allows to differentiate between functional divergence and convergence (Mason et al., 2005; de Bello et al., 2010; Laliberté and Legendre, 2010). While stochastic effects also influence patterns of functional trait variation, selection and adaptation to changing environment will lead to non-random patterns. Environmental changes are assumed to select for or against species with certain traits and result in functional convergence (species being more similar in their traits than expected by random) (Pillar et al., 2009). Competition can be assumed to lead to functional divergence due to niche differentiation (species being less similar in their traits than expected by random) (Purschke et al., 2013). Phylogenetic diversity (PD) reflects the variability of phylogenetic distances among the coexisting species and helps to differentiate between phylogenetic divergence (coexisting species are phylogenetically more distant than expected by random) and phylogenetic convergence (coexisting species are phylogenetically more similar than expected by random) (Grime, 2006; Vellend et al., 2011). In the case of phylogenetically conserved traits, FDa and PD should show similar patterns but trait variation could also be large within phylogenetic lineages (Prinzing et al., 2008). In the current study, we follow the example of Reitalu et al. (2015) and use three measures to characterize the functional and phylogenetic structure of community assembly: (1) functional composition representing the community weighted mean trait values (CWMs), (2) functional alpha diversity (FDa) reflecting the variability of traits among the coexisting species, (3) phylogenetic diversity (PD) reflecting the variability in distances between species/taxa in a phylogenetic tree. We test the applicability of these measures for pollen-based biodiversity reconstructions.

In addition to CWM, FDa, and PD, we will use Ellenberg indicator values (EIVs) (Ellenberg et al., 1991) to characterize plant/pollen assemblages in terms of their environmental tolerances. EIVs provide an understanding of environmental differences between sites where environmental conditions cannot be directly measured, such as in palaeoecological reconstructions (Kuneš et al., 2011; Reitalu et al., 2015; Felde and Birks, 2019).

To study palaeoecological changes in biodiversity, calcareous spring fens are suitable environments – they are recognized as hotspots of biodiversity (Bedford and Godwin, 2003; Hájková et al., 2015) and provide continuous sediment archives that preserve proxy records for biodiversity and environmental reconstructions. Calcareous fens are of ecological interest to both palaeo- and contemporary ecologists as they provide data of recent and long-term vegetation changes (Almendinger and Leete, 1998; Hájková and Hájek, 2003; Grootjans et al., 2006; Blaus et al., 2019; Carvalho et al., 2019; Jabłońska et al., 2019). Human activities such as artificial drainage, agricultural practices, and eutrophication have negatively influenced the directional development of calcareous spring fens leading to a gradual disappearance of characteristic species or habitats (Stammel et al., 2003; Topiæ and Stanèiæ, 2006; Osadowski et al., 2018). There are several palaeoecological studies of calcareous fens in Europe (e.g., Pidek et al., 2012; Hájková et al., 2012, 2015; Gałka et al., 2018; Jamrichová et al., 2018; Blaus et al., 2019; Jabłońska et al., 2019) but none of them have investigated modern pollen-plant diversity relationships. The current study of modern pollen-plant diversity from spring fens provides valuable information for palaeo reconstructions of fen vegetation by examining the functional and phylogenetic aspects of fen biodiversity offering a new angle on the development of these unique environments.

We focus on modern pollen–plant diversity relationships in spring fens from two different regions in Estonia: southern Estonia and Saaremaa Island in western Estonia. Pollen data is derived from moss polsters and plant cover data is estimated from the surrounding vegetation. We compare taxonomic richness, vegetation composition, FDa, trait CWM, PD and EIVs in modern pollen and plant data (Supplementary Figure S1). We aim to test the following hypotheses:

(1) Pollen-based richness and composition estimates of woody and herbaceous taxa are positively associated with corresponding plant-based estimates but at different spatial scales;

(2) Pollen- and plant-based CWM, FDa, PD and EIV estimates are positively correlated both in woody and herbaceous taxa;

(3) Pollen–plant diversity relationships do not differ between study areas (southern Estonia and western Estonia).

Further, we will use the information from the modern pollen–plant diversity relationships and apply the CWM, FDa, PD, and EIVs that have significant modern pollen–plant correlations to sedimentary pollen data (Supplementary Figure S1). Using existing well-dated pollen dataset from Kanna spring fen in western Estonia (Blaus et al., 2019), we will reconstruct the changes in different diversity components over the 9000-year history of the Kanna fen. Kanna fen developed first as a small bog and a shift to minerotrophic fen started ca 7000 years ago (Blaus et al., 2019). We expect to find high values of CWM of traits related to survival in nutrient-poor habitats (clonality, LDMC, mycorrhizal type) and tendency toward functional convergence due to environmental filtering in the initial bog phase of the mire development. In the minerotrophic fen phase during the mid-Holocene warm period, we expect to find high values of CWM of traits related to competitive abilities (SLA, plant height, seed size) and a tendency toward functional divergence due to niche differentiation. Among the Ellenberg indicator values, we expect the reaction (R), nutrient availability (N) and temperature (T) values to show the highest rate of change through the fen development.

Materials and Methods

Study Area

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.

FIGURE 1

Figure 1. Study area and sampling scheme. (A) Red dots indicate the sampled spring fens in southern Estonia and Saaremaa Island. Blue star marks the location of a Kanna spring fen used for sediment pollen based reconstructions. (B) The schematic overview of vegetation survey strategy modified from Bunting et al. (2013).

Pollen Data

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 Data

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.

Pollen−Plant Richness and Compositional Relationships in Different Radii

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.

Functional Diversity and Ellenberg Indicators

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).

TABLE 1

Table 1. Plant traits used to examine functional composition (CWM) and alpha diversity (FDa) of pollen and plant communities.

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

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.

Pollen-Plant Functional and Phylogenetic Diversity Comparisons

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.

TABLE 2

Table 2. Correlation test results for pollen- and plant-based community weighted mean (CWM) trait values (functional trait composition).

TABLE 3

Table 3. Correlation test results for pollen- and plant-based functional alpha diversity FDa estimates.

TABLE 4

Table 4. Correlation test results between pollen and plant based phylogenetic diversity.

TABLE 5

Table 5. Results of correlation tests between pollen-plant based Ellenberg indicators (EIVs) of light (L); temperature (T); continentality (K); moisture (F); soil pH reaction (R); nutrient availability (N) for woody and herbaceous plants.

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).

Sediment Pollen-Based Reconstructions

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).

Results

Pollen–Plant Richness and Compositional Relationships

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.

FIGURE 2

Figure 2. Pollen–plant richness and composition correlations. Y-axis: the values of Pearson’s correlation coefficients (for richness on the left) and Protest correlations (for vegetation composition on the right). X-axis: the radius of plant estimate. Vertical lines indicate bootstrapped 95% confidence intervals for the correlations.

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).

FIGURE 3

Figure 3. Association between pollen-based and plant-based herbaceous richness showing the influence of a geographical region [South (orange) vs. West (blue)]. Herbaceous plant richness at 6 m scale is used. Linear model with backward selection of variables indicates that plant richness (t = 4.26, p < 0.001) and region (t = –5.02, p < 0.001) are significant predictors of pollen richness (determination coefficient R² = 0.5).

Pollen-Plant Functional and Phylogenetic Diversity Comparisons

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).

Reconstructions in Time

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.

FIGURE 4

Figure 4. Sedimentary pollen based richness in Kanna spring fen. Generalized additive model of herbaceous taxa (green line) and woody taxa (red line) in time. Dashed lines enclose the 95% confidence intervals of trend lines.

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).

FIGURE 5

Figure 5. Sedimentary pollen based reconstructions of functional trait composition (CWM) in Kanna spring fen. Generalized additive models in relation to time for woody (red line) and herbaceous taxa (green line). Red and green dashed lines enclose the 95% confidence intervals of trend lines. Traits of herbaceous taxa is plotted on the secondary y-axes. Abbreviations of trait names are explained in Table 1.

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.

FIGURE 6

Figure 6. Sedimentary pollen based reconstructions of functional alpha diversity (FDa) in Kanna spring fen. Generalized additive models in relation to time for woody (red line) and herbaceous taxa (green line). FDa is expressed as a standardized effect size of mean pairwise distance (SESmpd). Red and green dashed lines enclose the 95% confidence intervals of trend lines. Abbreviations of the trait names are explained in Table 1.

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.

FIGURE 7

Figure 7. Sedimentary pollen based reconstruction of phylogenetic diversity (PD) in Kanna spring fen. Generalized additive models in relation to time. Dashed lines enclose the 95% confidence intervals of trend lines.

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.

Discussion

The use of sedimentary pollen to reconstruct functional and phylogenetic diversity is a relatively novel approach with previous studies focusing mainly on functional trait composition (community weighted mean trait values) (e.g., Reitalu et al., 2015; Brussel et al., 2018; Carvalho et al., 2019; van der Sande et al., 2019). In addition, Reitalu et al. (2015) calculated sedimentary pollen based functional alpha diversity and phylogenetic diversity that allow to study convergence and divergence patterns in community assembly. Our study adds to the evidence from modern pollen-plant studies of functional aspects (Carvalho et al., 2019) that both functional composition and alpha diversity can be studied using sedimentary pollen data and have the potential to provide valuable knowledge on palaeo-environments in addition to the conventional methods. However, not all studied variables had significant pollen-plant associations highlighting that uncritical use of different diversity metrics in pollen-based reconstructions may be misleading.

Taxonomic Diversity

Modern Pollen–Plant Richness and Composition

We hypothesized that modern pollen richness and composition reflect vegetation richness and composition, and that this relationship is dependent on plant group (woody vs. herbaceous) and spatial scale. The results partly supported our predictions showing that pollen richness of herbaceous taxa in moss polster samples best reflects herbaceous vegetation richness within 0.5–6 m radius and pollen-plant compositional relationship for herbaceous taxa was significant at all scales (Figure 2). The positive richness association existed in both study areas (Figure 3). This adds to the growing evidence that pollen taxonomic diversity is positively associated with the surrounding plant diversity (Urrego et al., 2011; Meltsov et al., 2012; Jantz et al., 2014; Matthias et al., 2015; Birks et al., 2016a; Mourelle and Prieto, 2016; Reitalu et al., 2019). In pollen analytical studies, the scale of the vegetation reflected in pollen sites is expressed as the relevant source area of pollen (RSAP). The RSAP is the area beyond which the pollen–plant correlation does not improve (Sugita, 1994). Our results for herbaceous taxa agree with existing studies from moss polsters where the RSAP of non-arboreal (non-woody) taxa is in the range of 0.5 – 10 m (Hjelle, 1998; Bunting, 2003; Bunting and Hjelle, 2010). When the RSAP calculation involves woody taxa, the RSAP estimations from the moss polsters are in the range of 300 – 400 m (Broström et al., 2005; Mazier et al., 2008) exceeding the scale of our vegetation survey. The lack of significant richness correlation for woody taxa in our study may be caused by the fact that our study radius was too small. However, the forests in the study regions are relatively uniform and increasing the radius would not add many new woody species. Pollen richness of woody taxa has been shown to be strongly positively associated with woody plant richness in a broad-scale study across entire N Europe (Reitalu et al., 2019) where the number of tree species varied from 3 to 80. In the current study, the richness of woody plant species varied from 11 to 22 and we can surmise that the variation in woody plant richness is too small to be reflected by the pollen data. However, the compositional relationship between pollen-plant data of woody taxa was significant at scales of 6 – 100 m radius indicating that vegetation composition of woody taxa is reflected by the pollen data.

While palaeo-diversity is traditionally calculated for all taxa combined (see overview in Birks et al., 2016a), our results support the studies suggesting to look at the diversity of wind-pollinated and insect-pollinated taxa (Meltsov et al., 2011) or at woody and herbaceous taxa separately (Carvalho et al., 2019; Reitalu et al., 2019). To conclude, dividing woody and herbaceous taxa in pollen-based diversity estimates improves our understanding of spatial scales reflected in palaeo-diversity reconstructions. Woody and herbaceous plants clearly form different functional groups (Díaz et al., 2016), and analyzing their diversity separately allows for better comparisons with contemporary plant ecological studies.

Sediment Pollen Based Richness Reconstruction

The richness of both woody and herbaceous taxa was relatively low in the beginning of the fen development from 9.2 to 7.2 ka and increased to the highest levels by ca 6 ka. Similar trend of increasing palynological richness in the beginning of the Holocene can be observed in other studies from Saaremaa (Poska and Saarse, 2002) and in the entire region (Reitalu et al., 2015). Palaeoecological studies from N Europe often show an increase in pollen richness during the last 2 or even 4 ka, associated with human-induced landscape opening (Birks and Line, 1992; Berglund et al., 2008; Reitalu et al., 2015). In the current study, herbaceous richness was relatively stable and woody richness decreased during the last millennia. Based on the modern pollen-plant correlations, herbaceous pollen richness from the fen samples is most likely to reflect local within-fen richness changes and taxa associated with agricultural activities and landscape opening coming from outside the fen are less likely to contribute. The decline in the richness of woody species is most likely a reflection of climate-driven decrease in broadleaved trees and increase in pine dominance that today has established a relatively homogeneous forest cover in the landscape.

Functional Diversity

Modern Pollen−Plant Relationships

We documented significant positive pollen-plant correlations for six CWM traits for woody taxa and for four CWM traits for herbaceous taxa (Table 2). Functional alpha diversity (FDa) had significant positive associations for five traits in woody taxa and for three traits in herbaceous taxa (Table 3). Woody taxa might show stronger relationships because fewer species are included in each pollen type, and when trait values for pollen types are calculated by trait averaging more information is lost for herbaceous species. In our dataset, the only woody pollen type including large number of species is Salix but for herbaceous taxa many pollen types have comparatively large numbers of species (e.g., Asteraceae n = 27; Cyperaceae n = 51; Poaceae n = 41; Orchidaceae n = 17). Similarly to our study, the modern pollen-plant study from fens in England showed that trait averaging for herbaceous pollen types might reduce the pollen-derived trait variability (Carvalho et al., 2019). Trait averaging is likely to have a stronger effect on the FDa (measured as mean trait distance between species/taxa in a community) than on the CWM that is already a result of averaging. FDa is particularly difficult to assess for taxon-trait combinations where the trait variation within pollen taxon is large (for example AM and plant height in Cyperaceae) leading to significant effects of trait averaging in pollen data. Even though mycorrhizal type, plant height and seed weight can be expected to be associated with waterlogging and nutrient stress in fen systems (Jabłońska et al., 2019), our modern pollen-plant analyses indicate that pollen might not be the best proxy to reconstruct the long-term changes in these traits. Plant macro remains (Felde and Birks, 2019; Jabłońska et al., 2019) or sedaDNA (Parducci et al., 2014) may give better taxonomic resolution. Another possible reason for the lack of significant correlation in modern pollen-plant data is related to short gradients within the study area, for example, EM associations are characteristic to the majority of present-day Estonian woody taxa and the CWM of EM varied between 0.8 and 1. Similarly to the missing woody taxon richness association, we can surmise that the variation is too small to be reflected by the pollen data and modern pollen-plant studies across larger geographic extent and longer diversity gradients would improve our understanding of the pollen plant diversity associations (Davis et al., 2013).

The diversity estimates that showed significant pollen-plant associations for woody taxa can be assumed to be associated with climate change (mycorrhizal type, clonality, height, LDMC, SLA), disturbance (height, clonality, seed number) and competitive strength (mycorrhizal type, height, clonality, SLA) (Cornelissen et al., 2003; Brundrett and Tedersoo, 2018). The woody taxon functional diversity – reflecting the landscape-scale changes – can therefore be used to reconstruct trait responses to climate change and long-term human disturbances and to quantify the role of competition for woody community assembly (Reitalu et al., 2015). As herbaceous pollen diversity mostly reflects change in local fen community, we can expect traits to respond to waterlogging and nutrient stress (SLA, LDMC, clonality) (Jabłońska et al., 2019), changes in soil reaction (ERM), disturbance by drainage (seed number) and to indicate competition within the fen community (clonality, SLA) (Cornelissen et al., 2003).

We hypothesized that the pollen–plant diversity relationships do not differ between the study areas (southern Estonia and western Estonia). This was true for most of the functional diversity estimates – out of the traits that had significant pollen-plant correlations, the interaction with region was only significant for herbaceous SLA (Tables 2 −3). The region effect on SLA is most likely the result of small variation in SLA in the fens of western region while in the southern region, there is much longer gradient in the SLA (Supplementary Figures S3k,ab). For the majority of traits, pollen-plant association was similar in both regions giving assurance that the diversity estimate can be used in palaeoreconstructions.

Sediment Pollen Based Functional Diversity Reconstruction

Woody taxa

The trait composition (CWM) of woody taxa showed low clonality, low proportion of AM associations, low SLA, and high LDMC during the first millennia of Kanna site development (9.2-7 ka) (Figure 5). Both SLA and LDMC have been widely used to reflect species’ response to soil fertility and climate change (Cornelissen et al., 2003; Laughlin et al., 2010) with SLA being positively and LDMC negatively associated with soil fertility and warmer climate. After the ice retreated, Saaremaa was under the waters of the Baltic Ice Lake and the highest parts of the island emerged after the drainage of the lake at about 10.3 ka (Poska and Saarse, 2002). The peat in Kanna site began to accumulate right after the site emerged from the water due to land uplift (Blaus et al., 2019) and the beginning of vegetation development represents establishment of the pioneer communities. Early site development was dominated by Pinus sylvestris (Blaus et al., 2019) that is a good example of species with resource conservation strategy by reducing water loss as an adaptation to nutrient-poor environments through low LDMC (Qin and Shangguan, 2019). Clonality, which is associated with abiotic stress (Klimeš et al., 1997; Ye et al., 2014), reached maximum around 8 ka. However, the variation in clonality of woody plants is relatively small and is likely to be mainly associated with encroachment by clonal shrubs like Salix and Myrica in the wet depressions.

Most of the trait CWMs showed relatively rapid turnover during 8-6 ka and the period between 6 and 4 ka is characterized by the highest SLA and AM values and the lowest LDMC for the entire period. These responses are likely to be associated with fundamental trade-offs related to nutrient conservation, nutrient acquisition and turnover (Wright et al., 2004), that in turn reflect warmer and drier climatic conditions during the mid-Holocene (Hammarlund et al., 2002) when favorable conditions for expansion of thermophilous tree species existed in Estonia (Saarse and Veski, 2001; Sillasoo et al., 2007). High CWM of AM coincides with immigration of some AM associated species e.g., Ulmus gabra (Thomas et al., 2018) and Fraxinus excelsior (Seven and Polle, 2014) and the trends are in line with findings that show correlation between warmer temperature and general AM fungal activity (Compant et al., 2010; Hempel et al., 2013). Around 2 ka, CWM of SLA and AM declined rapidly coinciding with climate cooling in the area (Hammarlund et al., 2002) and spruce (Picea abies) expansion in Saaremaa (Saarse et al., 1999). Overall, the trait CWM changes of woody taxa in our study area are compatible with reconstructed biome dynamics for Northern Europe based on plant functional types (Davis et al., 2015).

The FDa values of woody plants were above zero for most of the traits except clonality (Figure 6). Although most of the values were below significant trait overdispersion, there was a strong overall tendency toward niche differentiation in AM, seed number and SLA, especially in the beginning of the study period. This pattern is most likely associated with the many fundamental functional differences between gymnosperms and angiosperms (Brodribb et al., 2012) allowing for their long-term coexistence.

Herbaceous taxa

Between 9 and 7 ka Kanna fen developed as a small bog with a high abundance of Ericaceae species, presence of Poaceae and Cyperaceae (Blaus et al., 2019) with a shift to minerotrophic fen starting around 7 ka. For the “bog phase” of the mire development, we expected to find trait composition indicative of low productivity and adaptations to acidic conditions. In accordance with the expectations, our results indicated high values for ERM, seed number, LDMC, clonality and low values for SLA and plant height (Figure 5). Several of the traits have been hypothesized to mitigate nutrient scarcity in wetlands (LDMC, mycorrhizal type, clonality) (Moor et al., 2017; Jabłońska et al., 2019) and to be associated with waterlogging gradient (LDMC, SLA) (Violle et al., 2011; Baastrup-Spohr et al., 2015; Jabłońska et al., 2019). The “bog phase” was characterized by high proportion of Ericaceae taxa (that have low SLA, high LDMC and ERM) that are typical for nutrient poor and acidic soils (Cairney and Meharg, 2003).

Substantial changes in fen plant CWM traits during 7-5 ka coincide with the site transition to minerotrophic fen environment caused by the unique combination of climatic, topographic and hydrological conditions (Blaus et al., 2019). The “fen phase” of the mire development is characterized by increased water, calcium carbonate and organic matter content and occurrence of fen habitat specialists (e.g., Potentilla erecta, Parnassia palustris) as well as by high palynological richness (Blaus et al., 2019). High SLA and low LDMC during the fen phase are in good accordance with earlier studies showing that high SLA and low LDMC are associated with increase in soil reaction (Bartelheimer and Poschlod, 2015), in water and nutrient content (e.g., Díaz et al., 2004; Ordońez et al., 2010; Moor et al., 2015), anoxia (Klimkowska et al., 2019), organic matter accumulation (Nieder and Benbi, 2008), and species richness (Violle et al., 2011). High values for AM during the fen phase suggest similarly to Ramírez-Viga et al. (2018) that AM associations have an important role in wetlands.

After 2 ka, herbaceous CWM decreased for SLA and seed number but increased for LDMC and clonality. This functional change is related to the increasing proportions of Cyperaceae and the decline in the abundance of forbs that characterizes the “recent phase” of Kanna fen development (Blaus et al., 2019). Decreasing SLA and increasing LDMC have been associated with drier conditions in wetlands (Baastrup-Spohr et al., 2015). Jabłońska et al. (2019) showed that the importance of tussock-forming species with small clonal spread increased in fens prior to the termination of fen phase (and development of bog or mire woodland) and the reported functional changes might be an early “warning” for nature conservation. The actual reasons behind the recent functional changes and increase in Cyperaceae are rather unclear and it is difficult to draw a border between local and larger scale processes. Even if the changes are autogenic or successional these processes are usually following to allogenic disturbances (Bodini and Klotz, 2009). Besides, some studies have shown that climate induced hydrological changes in wetlands might have greater impact on plant community changes compared to autogenic processes (Wilcox, 2004; Dieleman et al., 2014).

Functional alpha diversity of ERM and seed number showed significant functional divergence (Figure 6). In grassland communities, functional divergence has been associated with niche differentiation at fine scales (within one square meter and less) where herbaceous species interact (de Bello et al., 2013). At larger spatial scales (within habitat patches), the functional divergence is likely a result of fine-scale within-habitat environmental heterogeneity where species with different traits occur side by side (de Bello et al., 2013; Bergholz et al., 2017). In addition to Ericaceae and other bog indicators, several typical fen species (Potentilla, Parnassia, Cladium mariscus) were present in low numbers in Kanna fen initial phase pollen samples (Blaus et al., 2019) and most probably grew on the edges of the small bog leading to the functional divergence pattern in the pollen signal.

Phylogenetic Diversity

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.

Ellenberg Indicator Values

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).

FIGURE 8

Figure 8. Sedimentary pollen based reconstructions of Ellenberg indicator values (EIVs) in Kanna spring fen. Generalized additive models in relation to time. Dashed lines enclose the 95% confidence intervals of trend lines.

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.

Conclusion

In a modern pollen-plant study from Estonian calcareous fens, our results suggest that while pollen of herbaceous taxa in fens reflects vegetation at local fen scale, the pollen of woody taxa is likely to reflect larger landscape scale forest vegetation. Dividing woody and herbaceous taxa in pollen-based diversity estimates improves our understanding of spatial scales reflected in palaeo-diversity reconstructions. Woody and herbaceous plants clearly form different functional groups (Díaz et al., 2016), and analyzing their diversity separately allows for better comparisons with contemporary plant ecological studies.

Correlations between pollen- and plant-based estimates of functional and phylogenetic diversity indicated that pollen data reflected the functional and phylogenetic aspects of plant communities reasonably well. For most of the indicators pollen-plant association did not differ between the two study regions. However, not all tested variables exhibited positive modern pollen-plant associations. The CWM values were better reflected in pollen data than the FDa values based on mean pairwise distances between taxa in a sample. The interpretation of functional variables that have large variation within pollen types may be difficult because information is lost by averaging trait values for pollen types. Testing for various pollen–plant diversity associations in our dataset may have been partly hindered by relatively short gradients in diversity variables and we, therefore, encourage further studies of modern pollen–plant diversity relationships not only from fens but across various environments to help to increase the reliability of interpreting historical processes from sedimentary pollen data.

Reconstructions of different diversity aspects in Kanna spring fen through 9.2 ka of fen development showed that diversity of woody taxa was closely related to the main climate changes with the mid-Holocene warm period having contrasting functional and phylogenetic diversity values compared to both early- and late-Holocene. The largest changes in the FD and PD of herbaceous taxa were related to the bog-fen transition around 7 ka. Our results show that climate and other abiotic processes have considerably influenced community FD and PD through the 9.2 ka. FDa of both woody and herbaceous taxa indicated functional divergence for several periods of community development where coexisting taxa were less similar in their traits than expected by random.

To conclude, pollen-based functional and phylogenetic diversity estimates and EIVs provide valuable knowledge in addition to the conventional pollen analysis. The use of different diversity metrics helps to achieve a better understanding of environmental effects on different strategy mechanisms of plants and to better interpret long-term processes affecting community assembly in various palaeoenvironments.

Data Availability Statement

The datasets generated for this study can be found in the Dryad data repository the Dryad Digital Repository (10.5061/dryad.wstqjq2hh).

Author Contributions

AB and TR were the principal conceivers of the manuscript with the study design and conceptual idea by TR and the leading of writing by AB, with the input from all authors. AB and TR conducted the fieldworks of modern vegetation and pollen sample data collection and performed the data analysis. JM provided the script for phylogenetic analysis. IH compiled mycorrhizal type data. PG, IH, JM, and SV reviewed the manuscript and contributed significantly with the comments and valuable insights. All authors contributed to the article and approved the submitted version.

Funding

The authors would like to acknowledge the funding from the Eesti Teadusagentuur (Estonian Research Council) for the research grants PUT1173 to TR, IUT1-8 and PRG323 to SV, PUT1170 to IH, and PUT1006 to PG and JM.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We are grateful to Mari Reitalu and Mari Lepik for assistance during the fieldwork and Lilian Freiberg for her local knowledge of the southern Estonian sites. We thank Jennifer Clear for the help on language editorial work.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2020.00207/full#supplementary-material

References

Almendinger, J. E., and Leete, J. H. (1998). Regional and local hydrogeology of calcareous fens in the Minnesota river basin, USA. Wetlands 18, 184–202. doi: 10.1007/BF03161655