Next Article in Journal
Transcriptome Analysis on the Mechanism of Ethylicin Inhibiting Pseudomonas syringae pv. actinidiae on Kiwifruit
Next Article in Special Issue
From Water into Sediment—Tracing Freshwater Cyanobacteria via DNA Analyses
Previous Article in Journal
The Microbiome-Gut-Brain Axis and Resilience to Developing Anxiety or Depression under Stress
Previous Article in Special Issue
Landscape Setting Drives the Microbial Eukaryotic Community Structure in Four Swedish Mountain Lakes over the Holocene
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Sedimentary Ancient DNA (sedaDNA) Reveals Fungal Diversity and Environmental Drivers of Community Changes throughout the Holocene in the Present Boreal Lake Lielais Svētiņu (Eastern Latvia)

1
Institute of Technology, University of Tartu, Nooruse 1, 50411 Tartu, Estonia
2
Department of Geography, Faculty of Geography and Earth Sciences, University of Latvia, Jelgavas iela 1, LV-1004 Riga, Latvia
3
Department of Geology, School of Science, Tallinn University of Technology, Ehitajate tee 5, 19086 Tallinn, Estonia
4
Institute of Ecology and Earth Sciences, University of Tartu, Ravila 14a, 50411 Tartu, Estonia
*
Author to whom correspondence should be addressed.
Microorganisms 2021, 9(4), 719; https://doi.org/10.3390/microorganisms9040719
Submission received: 17 March 2021 / Revised: 27 March 2021 / Accepted: 27 March 2021 / Published: 31 March 2021
(This article belongs to the Special Issue Understanding Ancient Microbiomes)

Abstract

:
Fungi are ecologically important in several ecosystem processes, yet their community composition, ecophysiological roles, and responses to changing environmental factors in historical sediments are rarely studied. Here we explored ancient fungal DNA from lake Lielais Svētiņu sediment throughout the Holocene (10.5 kyr) using the ITS metabarcoding approach. Our data revealed diverse fungal taxa and smooth community changes during most of the Holocene with rapid changes occurring in the last few millennia. More precisely, plankton parasitic fungi became more diverse from the Late Holocene (2–4 kyr) which could be related to a shift towards a cooler climate. The Latest Holocene (~2 kyr) showed a distinct increase in the richness of plankton parasites, mycorrhizal, and plant pathogenic fungi which can be associated with an increased transfer rate of plant material into the lake and blooms of planktonic organisms influenced by increased, yet moderate, human impact. Thus, major community shifts in plankton parasites and mycorrhizal fungi could be utilized as potential paleo-variables that accompany host-substrate dynamics. Our work demonstrates that fungal aDNA with predicted ecophysiology and host specificity can be employed to reconstruct both aquatic and surrounding terrestrial ecosystems and to estimate the influence of environmental change.

1. Introduction

In recent years, a growing number of studies have used lake sedimentary ancient DNA (sedaDNA) to reconstruct past ecosystem changes. Examples include studying the impact of invasive species [1], reconstructing a history of freshwater fisheries [2], and exploring floristic diversity changes [3,4]. However, the majority of these studies on past ecosystems and ecosystem changes have used ancient plant or animal DNA with much less attention paid to ancient fungi [5,6] which continue to be studied using traditional methods such as identification of spores or fossil remains [7,8]. To our knowledge, only a few studies from permafrost [5,6,9], sea sediments [10], lake [11,12,13], and cave sediments [14] have used molecular-based methods to assess the past fungal or whole eukaryotic diversity in historical ecosystems. Some of these studies [5,10,11,12,13] used universal 18S rDNA region primers that enable one to detect a wide variety of eukaryotes but restricts the number of identified fungal taxa. In contrast, Bellemain et al. [6] and Kochkina et al. [9] focused on the internal transcribed spacer (ITS) region that serves as the main DNA barcode marker used to assess fungal diversity [15,16]. ITS can discriminate between closely related species for a broad range of fungi compared to the nuclear ribosomal small subunit (SSU), which has poor species-level identification features in fungi.
Fungi are a highly diverse kingdom of eukaryotic organisms that play important ecological roles as saprobes, symbionts, or parasites of plants and animals. Therefore, changes in fungal communities over time indicate changes in the entire ecosystem and enable one to untangle the ecological status of the environment. Fungal richness in terrestrial environments, such as soil and plant-related habitats, is highly diverse [17]. Their spores and mycelia often transfer into lakes by soil particles, rain, wind, or together with other organic matter such as dung or plant material [18]. A recent study [11] showed that sedimentary aDNA, especially from fungi and plants, is useful for the reconstruction of the environment surrounding the lake. However, compared to terrestrial fungi, information about the biodiversity of aquatic fungi is still poorly documented. Nevertheless, fungi are common residents in aquatic habitats and represent one of the main groups of microbial eukaryote communities in deep-sea and lake sediments [10,11,19] where they play crucial roles in the turnover of organic matter. For example, parasitic chytrids act as important mediators in aquatic food webs and transfer nutrients from inedible phytoplankton to zooplankton (“mycoloop”) [20].
Here we used aDNA from lake sediments to assess the composition of past fungal communities, their functional roles, and community changes throughout the Holocene (in the last 10.5 kyr) in relation to environmental factors. We demonstrate that following a broad range of functional ecological roles of fungi offers a good proxy for interpretation of overall past biodiversity changes beyond fungal diversity. We hypothesize that the abundance and diversity of fungi with broad functionality (e.g., saprotrophs) and broad host range has been more stable over the Holocene while the diversity of fungi with narrow specificity of substrate or hosts has been more variable following the dynamics of host or specific substrate. To test our hypothesis and to achieve the above-mentioned general aims of this study, we: (a) re-constructed the past fungal community composition and their ecological diversity from lake sediments, (b) tested the predictive value of fungal diversity as a proxy for broad ecosystem changes, and (c) explored how ancient fungal communities responded to changes in climatic and environmental factors during the Holocene.

2. Materials and Methods

2.1. Sampling and Chronology of the Sediment Cores

The lake Lielais Svētiņu (LSv) is an important study site since it lies in a climatic and vegetative transition area and spans both the maritime-continentality line that runs west to east from the central Baltic area to NE Europe and transects the gradual decay of nemoral forest into boreal. Lielais Svētiņu has a long sedimentary record and well-characterized relatively late human impact, thus making it a useful site to decipher environmental changes in this climate zone (see Appendix A.1).
A sediment core was taken for aDNA analysis from the middle of the lake Lielais Svētiņu (LSv) (eastern Latvia; 56°45′ N, 27°08′ E, 4 m water depth; Figure 1) in March 2013 using a Russian type corer (1 m long, diameter of 10 cm). In total, an 11 m thick sediment was obtained. Independent samples (layers) were collected at 2.5 cm intervals from cores in the Institute of Geology in Tallinn, Estonia. Fifty-five layers (Table S1) were picked out for further analysis depending on the estimated past temperature changes [21]. Samples were collected more frequently from sediment areas with rapid temperature changes. From each individual layer, 3 biological replicate samples were taken in a row (see Figure A1 under Appendix A.6). All the subsampling was done using protective lab wear and clean equipment (Appendix A.2). The chronology of the LSv sediment core is based on radiocarbon dates from the sediment core described in Stivrins et al. [21], which were adjusted with our sediment core by the age-depth model (full description in Kisand et al. [11]).

2.2. Molecular Analysis

DNA extraction and amplification were performed in separate labs (buildings) dedicated for aDNA work under a cleaned UV-treated positive-flow hood (Kojair K-safety KR-125, Mänttä-Vippula, Finland) using clean equipment and protective lab wear (details in Appendix A.3). Total DNA was extracted from sediment samples in 3 biological replicates (total of 162 samples) using PowerSoil® DNA Isolation Kit (MoBio Laboratories, Carlsbad, CA, USA) according to the manufacturer’s protocol. Sediment samples homogenization was conducted using Fast PrepTM-24 Tissue Homogenizer (MP Biomedicals) for 55 s at 4 m/s. All extracted samples were stored at −20 °C until PCR amplification.
Metabarcoding libraries were performed in two steps (PCR1 and PCR2). PCR1 amplification targeted internal transcribed spacer 2 (ITS2) region using fungal specific multiplex of 5 forward primers (ITS3NGS1, ITS3NGS2, ITS3NGS3, ITS3NGS4, ITS3NGS5) and reverse primer ITS4NGS [17]. 0.75 μM of forward and reverse primers were added to 20 μL reactions containing 1–1.5 μL aDNA (depending on sediment depth), 2× Phusion Master Mix with HF buffer, 1 μg/μL BSA (Thermo Scientific, Vilnius, Lithuania), and ultrapure water (MoBio Laboratories, Carlsbad, CA, USA). Phusion High-Fidelity DNA polymerase (Thermo Scientific) was used for amplification program containing the following steps: denaturation at 98 °C for 30 s, followed by 30 cycles of 10 s at 98 °C, annealing at 46 °C for 30 s, and extension at 72 °C for 15 s, with a final extension at 72 °C for 10 min. The presence and quantity of amplicons were visualized using 1.4% agarose gels (1xTAE) containing ethidium bromide. All samples without visible amplicons on the gel were amplified three times, mixed, and purified/concentrated using Select-A-Size DNA Clean and ConcentratorTM kit (Zymo Research, The Epigenetics Company). Technical replicates (11 samples were amplified and sequenced twice) were used to test PCR and sequencing consistency. PCR2 was performed to tag PCR1 amplicons with Illumina TrueSeq adapters and P5/P7 tail indexing primers. Step 2 PCR reaction (20 μL) contained 2× Phusion Master Mix High-Fidelity buffer (Thermo Scientific), 25 μM Mplex primer, 2.5 μL long PCR index X, 1 μL of PCR amplicons, and nuclear-free water (Thermo Scientific). The Step 2 amplification followed the program: 2 min at 98 °C, followed by 12 cycles of 20 s at 98 °C, 30 s at 65 °C, 30 s at 72 °C; and final extension 5 min at 72 °C. In each step, negative controls (no added DNA) were used to avoid possible cross-contamination from reagents or the environment. Amplicons were sent to Illumina MiSeq 2 × 250 bp sequencing at the FIMM (Institute of Molecular Microbial Medicine Finland, Helsinki, Finland).

2.3. Bioinformatics Analysis

Raw reads were quality-trimmed using Trimmomatic (v. 0.32) (average quality score <30Q). The VSEARCH tool [22] was used for further reads pairing, dereplication, chimera removal, sorting, <4 read clusters removal, molecular operational taxonomic unit (mOTU) clustering at 97% similarity threshold (see Appendix A.4 for a detailed bioinformatics description). Additional post-clustering processing was conducted using the R-package LULU [23], which reduces erroneous mOTUs. The mOTUs taxonomic affiliation was obtained using the fungal reference database UNITE (UNITE ver. 7; [24]). We relied on 98%, 90%, 85%, 80%, and 75% sequence identity as a criterion for assigning OTUs to species, genus, family, order, or class level, respectively [17]. All non-fungal OTUs (e.g., Bacteria and Plantae) were discarded from further analyses. The sequence data has been submitted to the SRA database under BioProject PRJNA704193.

2.4. Assigning Ecological Roles of Fungi

The ecological roles of the fungi were cross-validated using various tools and composed to determine both the possible relationships between identified fungi and their hosts and identify fungal groups that are responsive to change in environmental factors. A number of tools/databases were leveraged for this purpose, including FUNGuild [25], FungalTraits [26], Plant Parasites of Europe [27], Fungal Families of the World [28], among others (Table S2). For each identified fungal mOTU, we assigned a trophic status, lifestyle, environment, mycorrhizal associations, and host-specificity.

2.5. Statistical Analysis

All statistical analyses were conducted using R version 3.5.1 [29]. Details about technical and biological replicates can be found in Appendix A.5. For downstream analysis, the technical replicates were excluded, and biological replicates were merged together (see Appendix A.5, Table S1).

2.5.1. Richness Changes between Time Periods

We used a Venn diagram (InteractiVenn; [30]) to visualize shared mOTU richness between four time periods: Latest Holocene (last 2.0 kyr), Late Holocene (2.0–4.0 kyr), Mid-Holocene (4.0–8.0 kyr), and Early Holocene (8.0–10.5 kyr). Principal component analysis (PCoA) (vegan package) was conducted to visualize the dynamics of the abundance of various mOTUs over time. In addition, we employed the PERMANOVA (adonis function in vegan R package; [31]) on mOTU normalized counts to test changes in the community between time periods.

2.5.2. Community Composition and Ecology of Fungi

We employed the codyn package [32] to measure community composition changes over the 10.5 kyr study period. For that we used the functions rate_change(), variance_ratio(), synchrony(), and turnover(). The rate change shows the rate of directional change in community composition over the study period using Euclidean distance [33]. The variance ratio and synchrony characterize patterns of species covariance and synchronous/asynchronous fluctuations. The variance ratio ranges from <1 (negative covariance) to >1 (positive covariance) and is centered at 1 (species vary independently). To detect asynchrony between species, we used “gross” synchrony metrics whose values range from −1 (complete asynchrony) to 1 (complete synchrony) and are centered at the value 0 which indicates that species fluctuate independently [34]. The turnover calculates the total species turnover metric, which shows the proportion of different species between two-time points [32].

2.5.3. Past Environmental Drivers

To test the relationship between past environmental factors and the richness of different fungal groups, we used the following paleoecological proxies: summer temperatures (Tsum), continental climate (Tcont), waterlogging (Wtol), drought (Dtol), and shade tolerance (Stol), the openness of surrounding area (OPEN), concentration of charcoal particles (Ch), pyrite (abundant in anoxic sediment, FeS2) and human-related pollen (HRP) [21]. Tsum, Tcont, OPEN, Wtol, Dtol, Stol, and HRP were reconstructed based on pollen records taken from the 2009 core [21]. Continental climate (Tcont) values were calculated as Tsum–Twin. Landscape openness (OPEN) was based on the pollen of dominant herbaceous taxa in the lake surrounding. The tolerance values of 14 trees/shrubs were used to reconstruct Wtol, Dtol, and Stol. Also, Charcoal particles (Ch) in the sediment were used as indications of fire events and pyrite (FeS2) as an indicator of anoxic conditions in the lake. HRP reconstruction was based on cereal species (Secale cereale, Hordeum vulgare, Triticum aestivum, and Avena sativa) that have been related to agricultural activities. The sediment core taken by Stivrins et al. [21] in 2009 from the same location can be compared with the sediment core used in this study. The data from two cores (2009 and 2013, this study) was parallelized using the age-depth model [11].
Raw reads data from different fungal groups were normalized (cumulative sums scaling—CSS, package metagenomeSeq, [35]), merged among biological replicates, and used for smoothing. mOTU data were combined among replicates and used for smoothing. We used the gam function from the mgcv package to fit the Generalized Additive Model (GAM) on nonlinear relationships. Smoothness estimation for the model was conducted using cubic splines (bs = “cr”), where k (knots) was chosen depending on the model fitting (the majority have 10 knots per 55 time points). GAM models were used to establish overall richness changes in fungal groups and to reduce variance over time (see Appendix A.5). We used fitted values from the model for two-sided Pearson correlation (cor.test function) with a confidence level of 0.95, which was used to study possible interactions between fungal groups and paleoecological proxies. Correlations were identified by correlation rate as very strong (>0.8 or <−0.8), strong (>0.6 or <−0.6), moderate (>0.45 or <−0.45) and weak (<0.45 or >−0.45) at p < 0.05.

3. Results

3.1. Fungal Sequences

Sequencing of the ITS2 region resulted in 53,006,690 raw read pairs of 55 triplicate samples (Table S1). After quality-filtering with Trimmomatic, 48,035,087 read pairs (90.6%) survived and were merged into 19,862,567 full-length reads. These were dereplicated into 1,539,800 identical sequence groups and all the reads <200 nt long were discarded from the dataset, leaving 19,091,590 reads. Further data denoising from chimeras and rare sequences (<4 reads) resulted in 16,361,004 reads. These sequences were clustered into 2708 molecular operational taxonomic units (mOTUs, <97% similarity). In addition, the post-clustering method LULU removed 20 mOTUs. BLASTn analysis gave no match to 1340 mOTUs and identified 115 non-fungal taxa (Plantae, Protista, and Protozoa), which were removed from downstream analyses. Out of 1233 fungal mOTUs (11,387,536 reads), 391 (32% of mOTUs) were determined as unidentified Fungi. In addition, we measured DNA fragmentation in the sediment samples using a DNA analyzer (TapeStation 2200) which showed that 300–350 bp long fragment amplification is possible even from deep sediment layers (Figure S1). For further statistical analysis, the technical replicates were removed. Unique mOTUs were combined and raw reads merged between biological replicates of each layer, forming a total of 55 samples with 1125 unique mOTUs (Table S1).

3.2. Overall Phylogenetic and Ecophysiological Richness of Fungi

Nearly 2/3 of mOTUs (68%) were classified at the Class level of taxonomy or lower (Figure 2A). In total, the fungal richness was greatest in phyla Ascomycota (30%) and Basidiomycota (27%), followed by a smaller proportion of mOTUs in Chytridiomycota (4%), Rozellomycota (6%), Zygomycota (2%), Glomeromycota (0.1%). The highest richness was recovered in the classes Agaricomycetes (96 genera), Sordariomycetes (35 genera), and Dothideomycetes (32 genera) (Figure S2). Chytridiomycota was represented by the classes Chytridiomycetes and Monoblepharidomycetes, yet most of these (67%), and all fungal mOTUs in Rozellomycota, remained unclassified at lower taxonomic levels. Most of the mOTUs (997 OTUs) were found in a small number of layers (1 to 5 layers; median = 1), while 61 mOTUs were detected in ≥10 layers (total of 55 layers) serving as a core community. The most frequent mOTUs throughout the sediment (in >40 layers) were fungi in the genera Cladosporium, Saccharomyces, and Ganoderma, which include mainly plant pathogens and saprotrophs. Further, many other aquatic and terrestrial saprotrophs (Hyaloraphidium sp., dung saprotroph Coprinellus sp., wood saprotroph Xylaria sp., yeast Candida sp., and Cryptococcus sp.) were observed throughout the sediment (in >30 layers), but not in any control samples. The full taxonomy is presented in Table S2.
In general, about 40% of the classified mOTUs belong to fungi that are expected to be terrestrial, while 23% were determined as aquatic and 12% as inhabiting both environments, leaving 25% of fungal habitats classified as unknown (Figure 3A). The most numerous guilds were terrestrial saprotrophs (294 mOTUs; 26%) and pathotrophs (245 mOTUs; 22%), in particular, wood and litter saprotrophs (128 mOTUs), as well as pathogens of plants (205 mOTUs), animals (140 mOTUs; including insects), and other fungi (61 mOTUs). Furthermore, we determined 40 mOTUs related to specific plant taxa and 23 mOTUs characteristic of certain forest types (conifer vs. deciduous trees) (see Table S2). Symbiotrophs were less numerous (41 mOTUs; 3.6%), including lichens and mycorrhizal fungi (arbuscular, ericoid, and ectomycorrhizal fungi) (Figure 3B). Among fungi with aquatic habitats, the most numerous were parasites of plankton (phytoplankton and other planktonic organisms) (223 mOTUs). In addition, we were able to detect a small number of fish pathogenic fungi (12 mOTUs). However, multiple ecological roles were assigned to 25.4% of mOTUs with the same identity and no ecological traits to 23% of mOTUs due to insufficient similarity to their described species or genera (see Appendix A.7).

3.3. Dynamics and Development of the Fungal Community

We conducted PCoA analysis and linear fitting of age and the concentration of human-related pollen (HRP) to visualize the dynamic changes in the abundance of mOTUs over time (Figure 2B and Figure S3). We see a smooth change in the community during the Holocene with more discontinuous and rapid changes during the Latest Holocene which distinctively contrasts from the rest of the samples. These results are also supported by a PERMANOVA analysis that shows significant differences between time periods (Table S3). We observed a three-fold increase in mOTUs during the Latest Holocene period (last 2 kyr, 37% of mOTUs, 2 new classes, Figure S4). Analysis of the shared mOTUs between time periods supports a model of gradual change in the community from the Early to Late Holocene (11 kyr to 2 kyr), while the Latest Holocene distinctively differed (Figure S4A). An increasing number of mOTUs were shared between consecutive time periods: 10% between Early and Mid-Holocene; 13% between Mid-Holocene and Late Holocene; 18% between Late and Latest Holocene).
To measure changes in the community composition over time, we used temporal diversity indices (rate change and turnover) and community stability metrics (variance ratio and synchrony) ([32]; Table 1, Figure 4). We observed a high total species turnover in the majority of fungal groups (mean ~80%) between adjacent time points (Figure 4) due to many appearing and disappearing species. Yet, the composition of fungal mOTUs that originate from terrestrial habitats was rather stable all over the 10.5 kyr (Table 1, Figure S5), which shows the importance of reappearing species. The highest rate of directional change in a community (rate change; [33]) appeared in the groups of aquatic fungi i.e., plankton parasites and in mycorrhizal fungi (Figure S5). In contrast, (especially in plankton parasitic fungi) the proportion of species turnover was gradually lower (Figure 4), illustrating a growing number of appearing species. Likewise, the variance ratio was high in these groups, which indicates species covariance over the entire timespan (Table 1). In addition, “gross” synchrony [34] showed that ~25% (median = 0.24; min–max: 0.14–0.26) of these species fluctuations depend on each other, thus illustrating that a quarter of species in the community responded in a similar way to environmental changes. In contrast, the highest asynchrony (32%) appeared in lichen symbionts, where species fluctuated opposite to each other. All the other larger ecological groups (terrestrial fungi, plant pathogens, saprotrophs, symbiotrophs) showed a relatively low rate change in community composition and “gross” synchrony values. Also, the species turnover was lower for saprotrophic fungi (i.e., wood, litter, and dung saprotrophs) (Figure 4).

3.4. Past Environmental Drivers

We used Generalized Additive Model (GAM) to denoise the variance in abundance data of fungal ecophysiological groups and the fitted values were related (correlation analysis) with available environmental paleoproxy variables (Figure 5 and Figure S6). Two-sided Pearson correlations between abundance in fungal ecophysiological groups and environmental paleoproxies illustrate significant positive correlations with continental climate (Tcont), openness of surrounding area (OPEN), concentration of charcoal particles (Ch), and human-related pollen (HRP), and significant negative correlation with summer temperatures (Tsum), shade tolerance (Stol), and pyrite (FeS2) (Figure 6, Table S4). In addition, plant pathogens also hold strong positive correlations with drought tolerance (Dtol) (r = 0.62) and strong negative correlations with Stol (r = −0.66) and Tsum (r = −0.61). Fish pathogenic fungi with low richness had a strong positive relationship to waterlogging tolerance (Wtol) (r = 0.68), whereas saprotrophs of litter, wood, and dung showed only a moderate correlation with different environmental paleoproxies.

4. Discussion

We demonstrated that fungal aDNA is well preserved in lake sediment shown to have a good record of various paleo-environmental and ecological materials including pollen, micro- and macrofossils [7,36]. As with the archived material analyzed by traditional paleoecological methods, fungal aDNA originates from both terrestrial and aquatic habitats and mirrors the community composition in the lake ecosystem and watershed. Therefore, when the ecophysiology, including the host specificity of fungal species archived in the lake sediment as aDNA, can be predicted, this information can be employed to reconstruct both the aquatic ecosystem and the ecosystem in the surrounding terrestrial landscape. Earlier, non-aDNA-based studies have demonstrated gradual ecosystem changes in this lake and its surroundings [21,36], driven by gradual climate changes during the Holocene. At the same time, an abundance of fungi with broad ecophysiology (e.g., saprotrophs) were relatively stable over the ~10.5 kyr study period and thus, more robust to environmental changes than narrow-substrate and host-specific fungi, such as pathotrophs (e.g., plankton parasites) or symbiotrophs [26,37,38], which follow the host and/or substrate dynamics. Contrastingly, more rapid ecosystem changes occurred during the last two millennia with a slow but steady increase in the anthropogenic influence in the area. Characteristic to this phenomenon is higher fungal richness among pathogenic fungi (plankton parasites, plant pathogens) and mycorrhizal fungi. Increased human activity is associated with an increased level of forest fires, floristic diversity, and landscape openness [36,39]. These environmental drivers have a different nature when compared to natural fluctuations in the structure of vegetation described by paleoproxies such as waterlogging, drought, and shade tolerance. In addition, although terrestrial and aquatic fungal communities were affected by the same environmental drivers, the more prevalent community changes appeared in aquatic fungi, possibly indicating a susceptibility of aquatic organisms to adjust to modest environmental disturbances.

4.1. Fungal Community Composition and Their Ecological Diversity

The general proportions in the fungal taxonomic affiliations at the phylum level (Figure 2A) were compared to previous studies from permafrost [6] and aquatic environments [18,19,40]. In contrast with sea sediments, we observed members of Glomeromycota and a higher proportion of Basidiomycota, which is more similar to terrestrial environments in this region [41]. Agaricales and Polyporales were the prevalent Basidiomycota orders (Figure S2). The species in these orders are mostly litter or wood saprotrophs, plant pathogens, or ectomycorrhizal symbionts (Agaricales) and likely represent transient species of terrestrial origin that have been carried to the lake with organic matter or spores. Thus, the high richness of terrestrial fungi preserved in lake sediments in all layers studied (Figure 3A), confirms the value of these archives to reconstruct the surrounding past terrestrial biodiversity [11], especially through interactions with their host species [7]. The high total fungal richness and successful species- and genus-level identification enabled the detection of host-specific fungi (46 mOTUs) (e.g., the plant pathogens of Alnus, Betula, Salix, Frangula, Picea, Poaceae, or mycorrhizal fungi of Pinus and Salix roots). The host-specific fungi offer valuable information about host populations in the past because fungi likely leave a stronger signal (higher concentration of DNA) within the sediment compared to host DNA (e.g., mammal DNA; [42]) due to their DNA being both widespread and protected by spore persistence in the sediment. Therefore, fungal host-specificity or associations to organism groups would enable us to use them as signatures, that could be used to study the appearance and dynamics of plant diseases in paleoenvironments.
In addition to terrestrial fungi, we also observed many Ascomycota, Chytridiomycota, and Rozellomycota characteristics of aquatic environments [19,43,44]. These fungi are known to play important ecological roles in aquatic habitats as saprobes on a variety of substrates or as parasites on phytoplankton, zooplankton, other fungi, or even aquatic invertebrates [45,46]. As a result, the presence of species from all 4 orders of Chytridiomycota (phytoplankton parasitic Rhizophydiales and Lobulomycetales, saprotrophic Monoblepharidales and Chytridiales) were observed in our samples. Also, multiple studies have shown Rozellomycota, which infects various organisms including large diatoms [44,47] and smaller planktonic eukaryotes [46,48].

4.2. Fungal Community Changes and Paleoecological Drivers

In conjunction with changes in climate and the dynamics of ecosystems (e.g., average temperatures, vegetation proxies in the region) we observed a gradual change in the fungal community during most of the Holocene (Figure 2B and Figure S4). Possibly due to the low number of mOTUs in the core community and relatively long time intervals obtained by subsampling the cores with a few centimeter intervals, the turnover at the species level was high (~80%) in most groups between consecutive time points (Figure 4). However, a low rate of community change was observed in most ecophysiological groups over the ~10.5 kyr (Figure S5). For example, we observed that terrestrial fungi (including about 40% of all mOTUs) had a relatively stable community over the entire timespan (Table 1, Figure 3A and Figure S5), exhibiting 3 times lower rate of community change than in aquatic fungi. Therefore, the aquatic fungi could be more susceptible to changes in environmental factors than terrestrial groups. Similarly, in the ecological groups of saprotrophs, animal, and fungal pathogens, we did not observe any major community changes over the Holocene (Table 1, Figure S5) and the proportion of species turnover was lowest for saprotrophic fungi (e.g., dung, wood and litter saprotrophs). All of these changes in richness can be considered natural fluctuations in diversity (Figure 3B). Thus, terrestrial groups of litter, wood, and dung saprotrophs demonstrated only moderate correlations with different environmental paleoproxies that do not induce community change. In addition, we show that fish pathogen fluctuations are correlated to waterlogging (Figure 6A) even though their communities did not drastically change. This may indicate possible variations in fish populations.
We observed that fungal richness accelerated in the last two or three millennia with an increase in the total number of mOTUs by ~37% (Figure 2B, Figure 3 and Figure S4). Most of these mOTUs belong to aquatic fungi (i.e., plankton parasites) with some exceptions. In contrast with the rise in richness, we observe a decrease in species turnover (Figure 4). This corroborates the formation of a community with newly added species and therefore more ecological niches in the lake. Our observation agrees with studies by Stivrins [21] and Kisand [11], who documented an increase in phytoplankton, especially Chlorophyta (green algae), over the last 2.0 kyr in LSv sediments. Phytoplankton species are the primary producers in lakes and zoosporic fungal parasites help to control the extent of their blooms [37]. Phytoplankton parasitic fungi have many possible roles in food web dynamics [45]. One role is contributing to nutrient transfer from inedible phytoplankton to zooplankton (“mycoloop”, [37]). Therefore, any shift in phytoplankton also impacts other aquatic organisms. In our study, the community change in plankton parasitic fungi is likely connected to an increase of phytoplankton or other planktonic organism richness.
In contrast with other terrestrial fungi, the species number (as mOTUs) of plant pathogens but also in groups of conifer-related fungi and mycorrhizal fungi increased during the last two millennia. Interestingly, the higher rate of community change appeared only in plankton parasites and in mycorrhizal fungi groups, where ~25% of plankton parasitic species and ~13% of mycorrhizal species fluctuated dependent on each other and respond similarly to environmental factors (Table 1, Figure S5). Thus, regardless of the rise in richness during the last ~2 kyr, and highly fluctuating species in adjacent timepoints (~60–80%), the communities of plant pathogens and conifer related fungi were rather stable over the entire timespan (Figure 3 and Figure 4, Table 1, Figure S5). Therefore, we conclude that some fungi with narrower host group or substrate specificity (e.g., plankton parasites, pathotrophs, and ectomycorrhizal fungi) could be more responsive to ecosystem changes than wide spectrum species such as saprotrophs or terrestrial fungi in general [37,38].

4.3. Increased Human Impact and Change in Richness of Fungi

The immediate surrounding of the lake under study has remained relatively pristine from direct human activity which can be associated with the appearance of HRP which gradually rose over the last 2 kyrs with a distinct increase during the last 700 yrs. Human activities have been shown to positively influence the vegetation richness of LSv [21] that possibly enriched the lake with additional nutrients (e.g., agriculture, grazing) which could explain the increase in richness of plant pathogenic and mycorrhizal fungi over the last 2 kyrs. Also, strong associations with drought indicate stress-related effects where vegetation is more prone to pathogens. In addition, the strong correlations to Ch and OPEN (Figure 6A) could also be results of increasing and widening human activity in the region, especially because Ch and OPEN are not the main factors impacting terrestrial communities in the late to early Holocene (2–10.5 kyr) (Table S5, Figure 6A,C). Ch is used as evidence of fire events in the past and is suggested to increase with drought [49] and human activities [36] such as slash-and-burn agriculture which alters the landscape. Forest fires are key factors that explain disturbances in forest and landscape structure [39,49] and can at least partly explain the increasing openness of the area (Figure 5A). Open areas enable more widespread transport of plant material to the lake and may reflect the larger-scale events. Therefore, a rapid increase in OPEN and Ch could be the result of wider human impact over the region (>2.0 kyr) whereas temperature, Wtol, Stol, and Dtol were mainly driven by natural factors. When excluding the last ~2 kyr the fungal communities were affected by waterlogging, shade, and drought tolerance as natural environmental drivers.
An increase in nutrients likely brought along the rise of phytoplankton and other planktonic organisms that further lead to an increase in the richness of plankton parasites in the lake. However, we also detected an increase in the richness of plankton parasitic fungi in the last ~4 kyr which is earlier than the time we see eukaryotic algae dominance (>2.0 kyr). This is possibly due to a late Holocene climatic-environmental change either driven by a shift in cyanobacteria-eukaryotic algae (eukaryotic algae dominance only appears at lower temperatures, [21]) or possibly due to cooler temperatures favored by parasitic fungi. Chytrid epidemics can appear even on optimal growth conditions for the host [37]. Therefore, the drivers of the observed community changes in plankton parasitic fungi are likely interactions with the shift to a cooler climate during Late Holocene and human activity and landscape openness in the Latest Holocene.

4.4. Challenges with Using Fungal Richness and Diversity as a Proxy for Total Paleo-Diversity

One of the first challenges of broad fungal community detection is efficient DNA extraction. Obtaining high yield and purity of aDNA preserved in the sediment is a demanding task due to the effects of humic-acid binding and the presence of clay components [50,51]. This can be achieved via taking extreme care with DNA extraction, PCR bias, possible contamination, and DNA degradation, among other concerns (see Appendix A.8). Another key challenge is to assign the habitat and ecophysiology to the resulting mOTUs based on a short ITS sequence because the origin of acquired fungal taxa is very broad in our sediment archives. The sequences can represent fungal diversity in sediment biomes, originate from water columns, or terrestrial environments [18]. One of the challenges of determining the fungal origin or their ecophysiology lies in the inability to assign a fungal identity to all detected fungi and unresolved taxonomy, especially for aquatic fungi. Therefore, species detection may be the key to acquiring reliable fungal signatures/proxies for paleo-environmental research. In addition, the habitat generalists may be less usable as a proxy than fungi with narrow ecological roles. Nevertheless, we were able to detect a wide variety of fungal taxa and their ecological roles and were able to reconstruct the past biodiversity. Despite this success, our work would benefit from additional research focused on generating reference sequences for databases and studies on fungal ecology. This additional information would enable us to detect an even more detailed understanding of fungal signatures in the sediment and increase their predictive value in ecosystem changes. Although the multiplex primer approach enables a reliable estimate of fungal richness, moving to a quantitative assessment of fungal biodiversity dynamics is an additional challenge due to high variation in ITS region repeats for different fungal species [52]. Still, we were able to determine the major community shifts in ecological groups and relate these to environmental changes. In addition, our data can be a reflection of changes in host diversity, for example, an increase of plankton parasites through plankton blooms. Molecular approaches that target host-specific fungal taxa and their quantification using qPCR would enable a more detailed understanding of the ecophysiology of fungi in specific habitats that are critical for paleoenvironmental research. Also, studies on the ITS area repeats and the status of DNA degradation would enable future researchers to make stronger conclusions about taxon abundances and therefore community shifts in the context of ecosystem change.

5. Conclusions

In summary, our study demonstrates that the ecophysiology of fungi archived into historical sediment records can be used to study the environmental drivers of past communities and processes in both lakes and their surrounding landscapes. We revealed various community changes and related these to paleoproxies and possibly to the changes in host dynamics of host-associated fungi. The observed major community shifts of plankton parasites and mycorrhizal fungi emerging together with increased human activity and climate shifts allow us to infer their susceptibility to environmental factors when most ubiquitous terrestrial groups (e.g., saprotrophs) with wide habitat ranges did not induce extreme community shifts. The effects of temperature, human activity, and landscape openness explained the community changes of fungal groups through plant material transfer and blooms of the planktonic organisms in the lakes. In addition, we highlight that fungi with narrow specificity of substrate or hosts could be more useful as proxies to follow the host and/or substrate dynamics compared with fungi that have broad functionality. A clearer understanding of the links between fungal species presence, abundance, and their host/habitat range over longer time periods needs to be established before we can broadly use ancient fungal communities as independent paleo-variables.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/microorganisms9040719/s1, Table S1: Samples taken from the sediment core and used for analysis. Table S2: Supplemental Data. Figure S1: Change of the DNA size distribution in various sediment depth: (a) 90 cm; (b) 200 cm; (c) 600 cm, (d) 900 cm measured with TapeStation 2200 (D5000 DNA chip). Figure S2: Detected classes and their richness changes over time in fungal phyla groups. Figure S3: Ordination space of PCoA analysis for read abundance (CSS normalized read counts) of mOTUs in divided time periods. Colored points mark groups of time periods and grey points mark the single mOTUs. Arrows indicate linear regression (envfit in vegan package) of different environmental factors. Table S3: PERMANOVA test analysis variation within time periods. Figure S4: Venn diagram demonstrates (A) mOTU variation and overlap between time periods (B) at class level variation and overlap between time periods. Figure S5: Showing the rate of community change in fungal ecological groups. rate_change_interval function generates the distances by time lag intervals. All the mOTUs detected in <2 samples were removed. Figure S6: mgcv plots of fitted values from GAM models of fungal groups. Table S4: Two-sided Pearson correlations between abundance in fungal groups and environmental paleoproxies using (A) normalized reads and (B) richness (mOTUs). All significant correlations (p < 0.05) are color-coded as a gradient from green (positive correlations) to red (negative correlations), not significant correlations are marked as n.s. Table S5: Two-sided Pearson correlations between abundance in fungal groups and environmental paleoproxies using richness (mOTUs), covering the Holocene period (2.0 kyr–10.5 kyr). All significant correlations (p < 0.05) are color-coded as a gradient from green (positive correlations) to red (negative correlations), not significant correlations are marked as n.s.

Author Contributions

Conceptualization, V.K. and L.T. (Liisi Talas); methodology, V.K., L.T. (Liisi Talas), and L.T. (Leho Tedersoo); formal analysis, L.T. (Liisi Talas), V.K., and N.S.; investigation, L.T. (Liisi Talas), V.K., S.V., L.T. (Leho Tedersoo), and N.S.; data curation, L.T. (Liisi Talas), V.K., and L.T. (Leho Tedersoo); writing—original draft preparation, L.T. (Liisi Talas) and V.K; writing—review and editing, L.T. (Liisi Talas), V.K., N.S., L.T. (Leho Tedersoo), and S.V.; visualization, L.T. (Liisi Talas) and V.K.; supervision, V.K. All authors have read and agreed to the published version of the manuscript.

Funding

Veljo Kisand was supported by grants PUT134 and PUT1389 (ETAg) and Institute of Technology, University of Tartu basic funding grant. Siim Veski was supported by IUT1-8 and PRG323 (ETAg). Normunds Stivrins was supported by the Latvian Council of Science project no. LZP-2020/2-.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The sequence data can be found from the SRA database under BioProject accession number PRJNA704193. The bioinformatic filtering scripts are available on Github at https://github.com/liisital/bioinformatic-filtering/blob/main/filtering-script.sh (accessed on 25 March 2021).

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Appendix A.1. Study Site Description

Sediment samples were collected from the lake Lielais Svētiņu (LSv) which is located in the Rezekne district of eastern Latvia. LSv is a lake with an area of 18.8 ha, a water depth of 4 m, and an elevation of 96.2 m above sea level. It was once part of a larger waterbody and is now a relic of an ancient lake [53]. The majority of the catchment area (~12 km2) of LSv is covered with forest, including, for example, species of Betula, Picea, Pinus, Ulmus, Tilia, Alnus and Quercus [36]. The current topography of LSv was created as a result of the Weichselian glaciation/deglaciation 11,700 years ago [54]. LSv is mostly characterized using common paleological analysis (pollen and non-pollen polymorphs) that show relatively late human impact in sediments [36,55].

Appendix A.2. Subsampling

Subsampling of the cores was finished within a one-month period at the Institute of Geology in Tallinn, Estonia. All the sediment surfaces that were in contact with oxygen were scraped off with a sterile spatula to expose the fresh sediment. Subsamples were collected with sterile 1 mL syringes and only half of that sample was placed into sterile 2 mL Eppendorf tubes to avoid possible contamination by sediment that entered the syringe first. All the subsample collecting was done using clean equipment under a positive-flow hood (Kojair K-safety KR-125), which was decontaminated with 98% ethanol and UV light prior to use. Subsamples were stored at −80 °C.

Appendix A.3. DNA Extraction and PCR

The DNA extraction and amplification were done in separate laboratories at the Institute of Technology, Tartu, Estonia. DNA extraction and PCR set-up were executed under separate positive-flow hoods (Kojair K-safety KR-125), which were beforehand swept with DNA AWAY™ Surface Decontaminant (Thermo Scientific, Vilnius, Lithuania) and exposed to UV light. All the surfaces and equipment were cleaned before each sampling process. All the pipetting, sample, or DNA handling were conducted only under the clean positive-flow hood (Kojair K-safety KR-125). The extraction controls (no added sediment) and negative PCR controls (no added DNA) were used to eliminate possible contamination from commercial extraction kits, used buffers, and from the environment.

Appendix A.4. Bioinformatic Analysis

Raw fastq file reads were quality-trimmed using Trimmomatic (V0.32) to remove Illumina-specific sequences and regions with low-sequence quality (<30Q cutoff). For further processing, VSERACH tool [22] was used. Quality trimmed reads were paired together (16bp overlap, --fastq_mergepairs), dereplicated and all the reads with length <200 bp were discarded. Thereafter the sequences were screened for chimeric fragments (--uchime_denovo) and detected chimeras were removed from the dataset. Unique identical sequence clusters were sorted by abundance and all the clusters containing <4 reads were discarded. Rare sequences have been shown to represent sequencing errors and technical artifacts [56]. Reads were clustered together into molecular operational taxonomic units (mOTU) using a VSEARCH command --cluster_smallmem, based on a 97% similarity threshold (vsearch --cluster_smallmem). Further post-clustering processing was carried out using R-package LULU [23]. LULU algorithm enables to accomplish ’more realistic’ biodiversity measures compared with only traditional denoising (removal of rare mOTUs) by reducing erroneous mOTUs and keeping real OTUs [23]. All mOTUs with BLASTn match under the <75% sequence similarity was determined as unidentifiable mOTUs.

Appendix A.5. Statistical Analyses: GAM Model

For GAM models we used summarized mOTU data and normalized reads data. Some of the groups (insect parasitic fungi, conifers-related fungi, mycorrhizal fungi, litter saprotrophs) showed no trends nor detectable dynamics and were therefore removed from the following correlation analysis. Two-sided Pearson correlations conducted with mOTU data and normalized reads showed similar results between correlations (only ± 0.1 to 0.3 correlation change between datasets) (Table S4).

Appendix A.6. Reproducibility of Amplicon-Based Community Descriptors

For general evaluation, we used the PERMANOVA test (adonis function in vegan R package; [31]) on mOTU richness and abundance to examine the variation between biological and technical replicates (Table A1). Results showed no significant variation (p = 1.0) between technical replicates in OTU richness or abundance, indicating no bias during PCR and sequencing steps. However, for analyzing which biological replicates, including samples depths as co-variate, may differ from each other, Linear Mixed Models (LMM; lme function in nlme R package; [57]) and likelihood ratio test were used to determine possible variation source and complexity by using raw read counts of mOTUs and a maximum number of mOTUs (Table A2). LMM showed variation in mOTU richness between biological replicates (p < 0.01), but no significant variation in read abundances (p = 0.76) or impact of sample depth (p = 0.1) in case of read abundances (Table A2). Thereafter, we used the Multcomp Tukey test (glht function in multcomp package; [58]) to show how raw read counts and max number of mOTU values differ between biological replicates (Table A3). We saw a shift from 1st replicate to 2nd and 3rd biological replicate in mOTU richness (Table A3, Figure A1), which could be due to vertical displacement of sediment layers where samples were collected or by natural differences between samples. Further, technical replicates were removed and unique mOTUs between biological replicates were combined, forming 55 samples (1125 unique mOTUs; 9,323,817 reads) taken from different depths.
Table A1. PERMANOVA test analysis variation within technical and biological replicates.
Table A1. PERMANOVA test analysis variation within technical and biological replicates.
Count of mOTUs
Adonis (formula = binary data ~ Tech_repl + Bio_repl,
data = metadata 3, permutations = 99, method = “jaccard”)
-Source of variationDfSumsOfSqsMeanSqsF.ModelR2Pr (>F)
Tech_repl 151.64450.328900.811510.075001.00
Bio_repl 222.04501.022522.522910.093260.01 **
Residuals4518.23820.40529-0.83174-
Total5221.9278--1.00000-
Raw Read Counts Data
Adonis (formula = read data ~ Tech_repl + Bio_repl,
data = metadata 3, permutations = 99, method = “bray”)
-Source of variationDfSumsOfSqsMeanSqsF.ModelR2Pr (>F)
Tech_repl 151.31180.262370.67430.057931.00
Bio_repl 223.82551.912744.91580.168920.01 **
Residuals4517.50950.38910-0.77315-
Total5222.6468--1.00000-
1 Tech_repl- Technical replicates; 2 Bio_repl- Biological replicates; 3 metadata- data frame for independent variables. Significance level p ≤ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1.
Table A2. LMM analysis and likelihood ratio test to determine the impact of depth to variation within biological replicates using (A) count of mOTUs and (B) raw reads.
Table A2. LMM analysis and likelihood ratio test to determine the impact of depth to variation within biological replicates using (A) count of mOTUs and (B) raw reads.
(A) Count of mOTUsm1rs = lme(chao ~ Biol 1 + Depth, random = ~1| Tech_repl 2, data = d, method = “ML”)
m1r_s = lme(chao ~ Biol * Depth, random = ~1|Tech_repl, data = d, method = “ML”)
anova (m1rs, m1r_s)
ModeldfAICBIClogLikTestL.Ratiop-Value
Biol + Depthm1rs171757.7041780.897−871.8521--
Biol * Depthm1r_s2101742.9461776.078−861.47311 vs. 2 20.758010.0001
summary (m1r_s)
Random effects:Formula: ~1 | Tech_repl
(Intercept)Residual
StdDev:4.03505216.69432
Fixed effects:chao ~ Repl * Depth2
ValueStd.ErrorDFt-Valuep-Value
(Intercept)81.723758.335641909.8041330.0000
Repl2−62.7425411.09800190−5.653500.0000
Repl3−23.7964111.14714190−2.1347540.0341
Depth2−0.065930.00994190−6.6306700.0000
Repl2:Depth20.058900.013711904.2958390.0000
Repl3:Depth20.012390.013801900.8973950.3706
(B) Raw Read Counts Datam3rs = lme(read_sum ~ Biol 1 + Depth, random = ~1|Tech_repl2, data = d2, method = “ML”)
m3r_s = lme(read_sum ~ Biol * Depth, random = ~1|Tech_repl, data = d2, method = “ML”)
anova (m3rs, m3r_s)
ModeldfAICBIClogLikTestL.Ratiop-Value
Biol + Depthm3rs175253.5775276.769−2619.789
Biol * Depthm3r_s2105253.3145286.446−2616.6571 vs. 2 6.2631130.0995
summary(m3r_s)
Random effects:Formula: ~1 | Tech_repl
-(Intercept)Residual-
StdDev:5.54790695,893.84
Fixed effects:read_sum ~ Repl * Depth2
-ValueStd.ErrorDFt-Valuep-Value-
(Intercept)108,192.1745,577.481902.37380790.0186-
Repl2−69,843.4263,663.73190−1.09706760.2740-
Repl396,198.0863,972.611901.50373860.1343-
Depth2−71.3457.09190−1.24961750.2130-
Repl2:Depth284.4878.641901.07424890.2841-
Repl3:Depth2−105.6379.22190−1.33332520.1840-
1 Tech_repl- Technical replicates; 2 Biol- Biological replicates.
Table A3. Multcomp Tukey test showing a correlation between biological replicates taken from a sediment core in the shown order.
Table A3. Multcomp Tukey test showing a correlation between biological replicates taken from a sediment core in the shown order.
-CountBiological
Replicates
EstimateStd. Errorz ValuePr (>|z|)
mOTUs2–1−17.6843.445−5.133<1 × 10−4 ***
3–1−15.2303.471−4.388<1 × 10−4 ***
3–22.4543.4320.7150.881
Raw Read2–1−453617,140−0.2650.993
3–111,78217,2660.6820.895
3–216,31717,0730.9560.757
Significance level p ≤ 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1.
Figure A1. The order of biological replicate sample collection from sediment layers.
Figure A1. The order of biological replicate sample collection from sediment layers.
Microorganisms 09 00719 g0a1

Appendix A.7. Ecology of Fungi

Of 25.4% of mOTUs with multiple ecological roles assigned, 11.7% of fungi were determined as pathotrophs-saprotrophs, which have parasitic to saprophytic lifestyles depending on environmental conditions. Likewise, 12% (135 mOTUs) of fungi were determined as pathotroph-saprotroph-symbiotroph, small amount as saprotroph-symbiotroph (1.4%; 16 mOTUs) and as pathotroph-symbiotroph (0.3%; 4 mOTUs) groups possibly having many adaptations in the environment. However, many of these fungi were not detectable at lower taxonomy than the order level.

Appendix A.8. Methodological Considerations

Possible contamination, richness bias, and DNA degradation are all very serious problems when trying to reconstruct realistic past community changes. The measures we took to prevent possible contamination were described in Appendix A.2 and Appendix A.3. For DNA extraction and amplification, the negative controls (no added sediment/DNA) were used to ensure no cross-contamination occurred from commercial kits or laboratory. In addition, only samples with no false-positive controls went to sequencing. After the bioinformatic workflow, only two fungal mOTUs (unidentified fungal mOTU, Ascomycota mOTU with 4 reads) were detected in PCR negative. Both mOTUs were removed from the dataset for further analysis.
We are discussing possible richness bias that can occur through (i) primer mismatches and taxonomic biases or (ii) bioinformatic workflow (e.g., quality trimming). ITS region has been proposed as the primary DNA marker for fungal barcoding from environmental samples [15,16]. However, commonly used ITS primers can introduce different types of primer biases (specificity to fungal group, taxonomy, mismatches, and length), therefore, it is suggested to use primer mixture of ITS region or the mixture of the nuclear ribosomal large subunit (LSU) with ITS region [16,59]. Thus, we used a combination of primers (five forward and one reverse primer that are analogous to ITS3 and ITS4) which were matching >99.5% of all fungi, including phylum Chytridiomycota and Glomeromycota [17]. We were able to detect a variety of fungi assigned to six phyla (Ascomycota, Basidiomycota, Chytridiomycota, Rozellomycota, Zygomycota, and Glomeromycota), however, a relatively high proportion of unmatched/unidentified clusters (32%) can reflect poor coverage of representatives in the sequence database. In addition, many mOTUs of Chytridiomycota and Rozellomycota remained unidentified at higher taxonomic levels due to the limited sequence availability of described species in UNITE database (v. 8). Still, the combined ITS region primer approach showed to be useful when targeting a wide range of fungal taxa in lake sediments. Final richness can be also affected by bioinformatic workflow through quality trimming, the efficiency of chimera removal, and clustering of mOTUs. Quality trimming (threshold > Q30) has been suggested for Illumina paired-end reads to avoid an increase in artificial “richness” and poorly overlapping reads that have been shown to appear using only modest quality trimming (threshold > Q20) [11]. Likewise, we used quality-filtering (threshold > Q30) to obtain successfully paired reads and to avoid potential sources of errors. From several tools available, we removed chimeras effectively using VSEARCH de novo tool [22]. Our comparative analysis for clustering showed a 10% decrease in obtained mOTUs when using a 97% similarity threshold compared to a 98% similarity threshold. Clustering with a higher similarity threshold (98%) resulted in separate mOTUs with the same taxonomic affiliation and GeneBank ID when the lower 97% threshold clustered these sequences together. Additional post-clustering processing (LULU) [23] reduced further 0.7% of clustered mOTUs achieving a “more realistic” survey of biodiversity by eliminating erroneous mOTUs.
Degradation of DNA in sediment is another issue for reconstructing paleoenvironment and community changes over time. We conducted DNA fragmentation measurements for sedimentary samples from different layers that showed targeting 300–350 bp DNA fragments will not bias the analysis severely even in old (and deep) sediment samples (Figure S1). When studying the changes in richness and community structures of fungal ecological groups, some large groups e.g., terrestrial fungi, showed relatively stable community over the whole timespan (Table 1, Figure S5). Therefore, regardless of the impact of depth to richness (Table A2) which seems to occur mainly due to high richness in the last ~2.5 kyr, we suggest that DNA degradation is not the main force creating richness fluctuations, but community changes are rather due environmental factors.

References

  1. Ficetola, G.F.; Poulenard, J.; Sabatier, P.; Messager, E.; Gielly, L.; Leloup, A.; Etienne, D.; Bakke, J.; Malet, E.; Fanget, B.; et al. DNA from Lake Sediments Reveals Long-Term Ecosystem Changes after a Biological Invasion. Sci. Adv. 2018, 4. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Nelson-Chorney, H.T.; Davis, C.S.; Poesch, M.S.; Vinebrooke, R.D.; Carli, C.M.; Taylor, M.K. Environmental DNA in lake sediment reveals biogeography of native genetic diversity. Front. Ecol. Environ. 2019, 17, 313–318. [Google Scholar] [CrossRef]
  3. Clarke, C.L.; Edwards, M.E.; Brown, A.G.; Gielly, L.; Lammers, Y.; Heintzman, P.D.; Ancin-Murguzur, F.J.; Bråthen, K.A.; Goslar, T.; Alsos, I.G. Holocene floristic diversity and richness in northeast Norway revealed by sedimentary ancient DNA (sedaDNA) and pollen. Boreas 2019, 48, 299–316. [Google Scholar] [CrossRef] [Green Version]
  4. Alsos, I.G.; Sjögren, P.; Edwards, M.E.; Landvik, J.Y.; Gielly, L.; Forwick, M.; Coissac, E.; Brown, A.G.; Jakobsen, L.V.; Føreid, M.K.; et al. Sedimentary ancient DNA from Lake Skartjørna, Svalbard: Assessing the resilience of arctic flora to Holocene climate change. Holocene 2015, 26, 627–642. [Google Scholar] [CrossRef] [Green Version]
  5. Lydolph, M.C.; Jacobsen, J.; Arctander, P.; Gilbert, M.T.P.; Gilichinsky, D.A.; Hansen, A.J.; Willerslev, E.; Lange, L. Beringian Paleoecology inferred from permafrost-preserved fungal DNA. Appl. Environ. Microbiol. 2005, 71, 1012–1017. [Google Scholar] [CrossRef] [Green Version]
  6. Bellemain, E.; Davey, M.L.; Kauserud, H.; Epp, L.S.; Boessenkool, S.; Coissac, E.; Geml, J.; Edwards, M.; Willerslev, E.; Gussarova, G.; et al. Fungal palaeodiversity revealed using high-throughput metabarcoding of ancient DNA from arctic permafrost. Environ. Microbiol. 2013, 15, 1176–1189. [Google Scholar] [CrossRef]
  7. Doyen, E.; Etienne, D. Ecological and human land-use indicator value of fungal spore morphotypes and assemblages. Veg. Hist. Archaeob. 2017, 26, 357–367. [Google Scholar] [CrossRef]
  8. Stivrins, N.; Aakala, T.; Ilvonen, L.; Pasanen, L.; Kuuluvainen, T.; Vasander, H.; Gałka, M.; Disbrey, H.R.; Liepins, J.; Holmström, L.; et al. Integrating fire-scar, charcoal and fungal spore data to study fire events in the boreal forest of northern Europe. Holocene 2019, 29, 1480–1490. [Google Scholar] [CrossRef]
  9. Kochkina, G.; Ivanushkina, N.; Ozerskaya, S.; Chigineva, N.; Vasilenko, O.; Firsov, S.; Spirina, E.; Gilichinsky, D. Ancient fungi in Antarctic permafrost environments. FEMS Microbiol. Ecol. 2012, 82, 501–509. [Google Scholar] [CrossRef] [Green Version]
  10. De Schepper, S.; Ray, J.L.; Skaar, K.S.; Sadatzki, H.; Ijaz, U.Z.; Stein, R.; Larsen, A. The potential of sedimentary ancient DNA for reconstructing past sea ice evolution. ISME J. 2019, 13, 2566–2577. [Google Scholar] [CrossRef] [Green Version]
  11. Kisand, V.; Talas, L.; Kisand, A.; Stivrins, N.; Reitalu, T.; Alliksaar, T.; Vassiljev, J.; Liiv, M.; Heinsalu, A.; Seppä, H.; et al. From microbial eukaryotes to metazoan vertebrates: Wide spectrum paleo-diversity in sedimentary ancient DNA over the last ~14,500 years. Geobiology 2018, 16, 628–639. [Google Scholar] [CrossRef]
  12. Capo, E.; Debroas, D.; Arnaud, F.; Guillemot, T.; Bichet, V.; Millet, L.; Gauthier, E.; Massa, C.; Develle, A.L.; Pignol, C.; et al. Long-term dynamics in microbial eukaryotes communities: A palaeolimnological view based on sedimentary DNA. Mol. Ecol. 2016, 25, 5925–5943. [Google Scholar] [CrossRef]
  13. Capo, E.; Ninnes, S.; Domaizon, I.; Bertilsson, S.; Bigler, C.; Wang, X.R.; Bindler, R.; Rydberg, J. Landscape Setting Drives the Microbial Eukaryotic Community Structure in Four Swedish Mountain Lakes over the Holocene. Microorganisms 2021, 9, 355. [Google Scholar] [CrossRef]
  14. Epure, L.; Meleg, I.N.; Munteanu, C.-M.; Dumitru, R.; Oana, R.; Moldovan, T. Bacterial and Fungal Diversity of Quaternary Cave Sediment Deposits. Geomicrobiol. J. 2014, 31, 116–127. [Google Scholar] [CrossRef]
  15. Schoch, C.L.; Seifert, K.A.; Huhndorf, S.; Robert, V.; Spouge, J.L.; Levesque, C.A. Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi. Proc. Natl. Acad. Sci. USA 2012, 109, 6241–6246. [Google Scholar] [CrossRef] [Green Version]
  16. Bellemain, E.; Carlsen, T.; Brochmann, C.; Coissac, E.; Taberlet, P.; Kauserud, H. ITS as an environmental DNA barcode for fungi: An in silico approach reveals potential PCR biases. BMC Microbiol. 2010, 10, 189. [Google Scholar] [CrossRef] [Green Version]
  17. Tedersoo, L.; Bahram, M.; Põlme, S.; Kõljalg, U.; Yorou, N.S.; Wijesundera, R.; Ruiz, L.V.; Vasco-Palacios, A.M.; Thu, P.Q.; Suija, A.; et al. Global diversity and geography of soil fungi. Science 2014, 346, 1256688. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Khomich, M.; Davey, M.L.; Kauserud, H.; Rasconi, S.; Andersen, T. Fungal communities in Scandinavian lakes along a longitudinal gradient. Fungal Ecol. 2017, 27, 36–46. [Google Scholar] [CrossRef]
  19. Xu, W.; Gao, Y.; Gong, L.; Li, M.; Pang, K.L.; Luo, Z.H. Fungal diversity in the deep-sea hadal sediments of the Yap Trench by cultivation and high throughput sequencing methods based on ITS rRNA gene. Deep Sea Res. Part I Oceanogr. Res. Pap. 2019, 145, 125–136. [Google Scholar] [CrossRef]
  20. Kagami, M.; Miki, T.; Takimoto, G. Mycoloop: Chytrids in aquatic food webs. Front. Microbiol. 2014, 5, PMC4001071. [Google Scholar] [CrossRef]
  21. Stivrins, N.; Kołaczek, P.; Reitalu, T.; Seppä, H.; Veski, S. Phytoplankton response to the environmental and climatic variability in a temperate lake over the last 14,500 years in eastern Latvia. J. Paleolimnol. 2015, 103–119. [Google Scholar] [CrossRef]
  22. Rognes, T.; Flouri, T.; Nichols, B.; Quince, C.; Mahé, F. VSEARCH: A versatile open source tool for metagenomics. PeerJ 2016, 2016, e2584. [Google Scholar] [CrossRef] [PubMed]
  23. Frøslev, T.G.; Kjøller, R.; Bruun, H.H.; Ejrnæs, R.; Brunbjerg, A.K.; Pietroni, C.; Hansen, A.J. Algorithm for post-clustering curation of DNA amplicon data yields reliable biodiversity estimates. Nat. Commun. 2017, 8, 1188. [Google Scholar] [CrossRef] [PubMed]
  24. Kõljalg, U.; Nilsson, R.H.; Abarenkov, K.; Tedersoo, L.; Taylor, A.F.S.; Bahram, M.; Bates, S.T.; Bruns, T.D.; Bengtsson-Palme, J.; Callaghan, T.M.; et al. Towards a unified paradigm for sequence-based identification of fungi. Mol. Ecol. 2013, 22, 5271–5277. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  25. Nguyen, N.H.; Song, Z.; Bates, S.T.; Branco, S.; Tedersoo, L.; Menke, J.; Schilling, J.S.; Kennedy, P.G. FUNGuild: An open annotation tool for parsing fungal community datasets by ecological guild. Fungal Ecol. 2016, 20, 241–248. [Google Scholar] [CrossRef]
  26. Põlme, S.; Bahram, M.; Jacquemyn, H.; Kennedy, P.; Kohout, P.; Moora, M.; Oja, J.; Öpik, M.; Pecoraro, L.; Tedersoo, L. Host preference and network properties in biotrophic plant–fungal associations. New Phytol. 2018, 217, 1230–1239. [Google Scholar] [CrossRef] [Green Version]
  27. Ellis, W. Plant Parasites of Europe. Available online: https://bladmineerders.nl/ (accessed on 10 September 2019).
  28. Cannon, P.; Kirk, P. Fungal Families of the World; CABI: Wallingford, UK, 2007; pp. 1–456. [Google Scholar]
  29. R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, 2018. Available online: http://www.r-project.org/ (accessed on 29 March 2021).
  30. Heberle, H.; Meirelles, V.G.; da Silva, F.R.; Telles, G.P.; Minghim, R. InteractiVenn: A web-based tool for the analysis of sets through Venn diagrams. BMC Bioinform. 2015, 16, 1–7. [Google Scholar] [CrossRef] [PubMed]
  31. Oksanen, J.; Blanchet, F.; Friendly, M.; Kindt, R.; Legendre, P.; McGlinn, D.; Minchin, P.; O’Hara, R.; Simpson, G.; Solymos, P.; et al. Vegan: Community Ecology Package. R package version 2.5–4. 2019. Available online: https://cran.r-project.org/package=vegan (accessed on 29 March 2021).
  32. Hallett, L.M.; Jones, S.K.; MacDonald, A.A.M.; Jones, M.B.; Flynn, D.F.B.; Ripplinger, J.; Slaughter, P.; Gries, C.; Collins, S.L. Codyn: An R package of community dynamics metrics. Methods Ecol. Evol. 2016, 7, 1146–1151. [Google Scholar] [CrossRef]
  33. Collins, S.L.; Micheli, F.; Hartt, L.; Collins, S.L. A Method to Determine Rates and Patterns of Variability in Ecological Communities. Oikos 2000, 91, 285–293. Available online: https://www.jstor.org/stable/3547549 (accessed on 29 March 2021). [CrossRef]
  34. Gross, K.; Cardinale, B.J.; Fox, J.W.; Gonzalez, A.; Loreau, M.; Polley, H.W.; Reich, P.B.; Van Ruijven, J. Species richness and the temporal stability of biomass production: A new analysis of recent biodiversity experiments. Am. Nat. 2014, 183, 1–12. [Google Scholar] [CrossRef] [Green Version]
  35. Paulson, J.N.; Colin Stine, O.; Bravo, H.C.; Pop, M. Differential abundance analysis for microbial marker-gene surveys. Nat. Methods 2013, 10, 1200–1202. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Stivrins, N.; Kalnina, L.; Veski, S.; Zeimule, S. Local and regional Holocene vegetation dynamics at two sites in eastern Latvia. Boreal Environ. Res. 2014, 19, 310–322. [Google Scholar]
  37. Kagami, M.; De Bruin, A.; Ibelings, B.W.; Van Donk, E. Parasitic chytrids: Their effects on phytoplankton communities and food-web dynamics. Hydrobiologia 2007, 578, 113–129. [Google Scholar] [CrossRef] [Green Version]
  38. Větrovský, T.; Kohout, P.; Kopecký, M.; Machac, A.; Man, M.; Bahnmann, B.D.; Brabcová, V.; Choi, J.; Meszárošová, L.; Human, Z.R.; et al. A meta-analysis of global fungal distribution reveals climate-driven patterns. Nat. Commun. 2019, 10, 5142. [Google Scholar] [CrossRef] [Green Version]
  39. Feurdean, A.; Veski, S.; Florescu, G.; Vannière, B.; Pfeiffer, M.; O’Hara, R.B.; Stivrins, N.; Amon, L.; Heinsalu, A.; Vassiljev, J.; et al. Broadleaf deciduous forest counterbalanced the direct effect of climate on Holocene fire regime in hemiboreal/boreal region (NE Europe). Quat. Sci. Rev. 2017, 169, 378–390. [Google Scholar] [CrossRef]
  40. Röhl, O.; Peršoh, D.; Mittelbach, M.; Elbrecht, V.; Brachmann, A.; Nuy, J.; Boenigk, J.; Leese, F.; Begerow, D. Distinct sensitivity of fungal freshwater guilds to water quality. Mycol. Prog. 2017, 16, 155–169. [Google Scholar] [CrossRef]
  41. Tedersoo, L.; Anslan, S.; Bahram, M.; Drenkhan, R.; Pritsch, K.; Buegger, F.; Padari, A.; Hagh-Doust, N.; Mikryukov, V.; Gohar, D.; et al. Regional-Scale In-Depth Analysis of Soil Fungal Diversity Reveals Strong pH and Plant Species Effects in Northern Europe. Front. Microbiol. 2020, 11, 1953. [Google Scholar] [CrossRef]
  42. Giguet-Covex, C.; Ficetola, G.F.; Walsh, K.; Poulenard, J.; Bajard, M.; Fouinat, L.; Sabatier, P.; Gielly, L.; Messager, E.; Develle, A.L.; et al. New insights on lake sediment DNA from the catchment: Importance of taphonomic and analytical issues on the record quality. Sci. Rep. 2019, 9, 14676. [Google Scholar] [CrossRef]
  43. Wurzbacher, C.; Rösel, S.; Rychła, A.; Grossart, H.P. Importance of Saprotrophic Freshwater Fungi for Pollen Degradation. PLoS ONE 2014, 9, e94643. [Google Scholar] [CrossRef] [Green Version]
  44. Grossart, H.P.; Wurzbacher, C.; James, T.Y.; Kagami, M. Discovery of dark matter fungi in aquatic ecosystems demands a reappraisal of the phylogeny and ecology of zoosporic fungi. Fungal Ecol. 2016, 19, 28–38. [Google Scholar] [CrossRef] [Green Version]
  45. Gleason, F.H.; Kagami, M.; Lefevre, E.; Sime-Ngando, T. The ecology of chytrids in aquatic ecosystems: Roles in food web dynamics. Fungal Biol. Rev. 2008, 22, 17–25. [Google Scholar] [CrossRef]
  46. Wurzbacher, C.; Warthmann, N.; Bourne, E.C.; Attermeyer, K.; Allgaier, M.; Powell, J.R.; Detering, H.; Mbedi, S.; Grossart, H.P.; Monaghan, M.T. High habitat-specificity in fungal communities in oligo-mesotrophic, temperate Lake Stechlin (North-East Germany). MycoKeys 2016, 16, 17–44. [Google Scholar] [CrossRef] [Green Version]
  47. Ishida, S.; Nozaki, D.; Grossart, H.P.; Kagami, M. Novel basal, fungal lineages from freshwater phytoplankton and lake samples. Environ. Microbiol. Rep. 2015, 7, 435–441. [Google Scholar] [CrossRef]
  48. Jones, M.D.M.; Forn, I.; Gadelha, C.; Egan, M.J.; Bass, D.; Massana, R.; Richards, T.A. Discovery of novel intermediate forms redefines the fungal tree of life. Nature 2011, 474, 200–204. [Google Scholar] [CrossRef]
  49. Aakala, T. Forest fire histories and tree age structures in Värriö and Maltio Strict Nature Reserves, northern Finland. Boreal Environ. Res. 2018, 23, 209–219. [Google Scholar]
  50. Hurt, R.A.; Qiu, X.; Wu, L.; Roh, Y.; Palumbo, A.V.; Tiedje, J.M.; Zhou, J. Simultaneous Recovery of RNA and DNA from Soils and Sediments. Appl. Environ. Microbiol. 2001, 67, 4495–4503. [Google Scholar] [CrossRef] [Green Version]
  51. Lakay, F.M.; Botha, A.; Prior, B.A. Comparative analysis of environmental DNA extraction and purification methods from different humic acid-rich soils. J. Appl. Microbiol. 2007, 102, 265–273. [Google Scholar] [CrossRef]
  52. Stadler, M.; Lambert, C.; Wibberg, D.; Kalinowski, J.; Cox, R.J.; Kolařík, M.; Kuhnert, E. Intragenomic polymorphisms in the ITS region of high-quality genomes of the Hypoxylaceae (Xylariales, Ascomycota). Mycol. Progress 2020, 19, 235–245. [Google Scholar] [CrossRef] [Green Version]
  53. Seglins, V.; Kalnina, L.; Lacis, A. The Lubans Plain, Latvia, as a reference area for long term studies of human impact on the environment. PACT 1999, 57, 105–129. [Google Scholar]
  54. Zelčs, V.; Markots, A. Deglaciation history of Latvia. In Quaternary Glaciations—Extent and Chronology of Glaciations, 1st ed.; Ehlers, J., Gibbard, P.L., Eds.; Elsiever: Amsterdam, The Netherlands, 2004; Volume 2, pp. 225–243. [Google Scholar]
  55. Veski, S.; Amon, L.; Heinsalu, A.; Reitalu, T.; Saarse, L.; Stivrins, N.; Vassiljev, J. Lateglacial vegetation dynamics in the eastern Baltic region between 14,500 and 11 400 cal yr BP: A complete record since the Bølling (GI-1e) to the Holocene. Quat. Sci. Rev. 2012, 40, 39–53. [Google Scholar] [CrossRef]
  56. Unterseher, M.; Jumpponen, A.; Öpik, M.; Tedersoo, L.; Moora, M.; Dormann, C.F.; Schnittler, M. Species abundance distributions and richness estimations in fungal metagenomics—Lessons learned from community ecology. Mol. Ecol. 2011, 20, 275–285. [Google Scholar] [CrossRef]
  57. Pinheiro, J.; Bates, D.; DebRoy, S.; Sarkar, D.; R Core Team. Nlme: Linear and nonlinear mixed effects models. R package version 3. 2019, pp. 1–140. Available online: https://CRAN.R-project.org/package=nlme (accessed on 15 October 2019).
  58. Hothorn, T.; Bretz, F.; Westfall, P.; Heiberger, R.M.; Schuetzenmeister, A.; Scheibe, S. Multcomp: Simultaneous inference in general parametric models. R package version 1.4–5. 2016. Available online: http://cran.stat.sfu.ca/web/packages/multcomp/multcomp.pdf (accessed on 29 March 2021).
  59. Frenken, T.; Alacid, E.; Berger, S.A.; Bourne, E.C.; Gerphagnon, M.; Grossart, H.-P.; Gsell, A.S.; Ibelings, B.W.; Kagami, M.; Küpper, F.C. Integrating chytrid fungal parasites into plankton ecology: Research gaps and needs. Environ. Microbiol. 2017, 19, 3802–3822. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Map showing the location of the lake Lielais Svētiņu and coring site of analyzed sediment cores taken March 2013 in Latvia.
Figure 1. Map showing the location of the lake Lielais Svētiņu and coring site of analyzed sediment cores taken March 2013 in Latvia.
Microorganisms 09 00719 g001
Figure 2. Dynamics of the fungi in the sediment core: (A) Richness of observed fungi grouped by phylum; (B) Ordination space of PCoA analysis using read abundance (normalized read counts) of mOTUs. Colored points mark various time periods in Holocene and grey points mark the single mOTUs. Arrows indicate linear regression (envfit in vegan package) of sample age (Age, R2 = 0.57, p < 0.05) and concentration of human-related pollen (HRP, R2 = 0.43, p < 0.05) with PCoA space.
Figure 2. Dynamics of the fungi in the sediment core: (A) Richness of observed fungi grouped by phylum; (B) Ordination space of PCoA analysis using read abundance (normalized read counts) of mOTUs. Colored points mark various time periods in Holocene and grey points mark the single mOTUs. Arrows indicate linear regression (envfit in vegan package) of sample age (Age, R2 = 0.57, p < 0.05) and concentration of human-related pollen (HRP, R2 = 0.43, p < 0.05) with PCoA space.
Microorganisms 09 00719 g002
Figure 3. Changes in (A) richness of fungi and (B) in distributed ecophysiological groups over 10.5 kyr. Dotted lines with colored backgrounds illustrate time periods.
Figure 3. Changes in (A) richness of fungi and (B) in distributed ecophysiological groups over 10.5 kyr. Dotted lines with colored backgrounds illustrate time periods.
Microorganisms 09 00719 g003
Figure 4. The total species turnover between consecutive sampling (time) points of fungal groups over the last 10.5 kyr. Total species turnover shows the proportion of different species between two time points.
Figure 4. The total species turnover between consecutive sampling (time) points of fungal groups over the last 10.5 kyr. Total species turnover shows the proportion of different species between two time points.
Microorganisms 09 00719 g004
Figure 5. The dynamics of (A) environmental paleoproxies and (B) ecophysiological groups of fungi over the last 10.5 kyr, modeled using GAM with formula = y ~ splines::bs(x, 5). Assigned ecophysiological roles of fungi abbreviations: plankton parasites (Plankton), animal pathogens (Animal), fish pathogens (Fish), insect parasites (Insect), fungal parasites (Fungal), plant pathogens (Plant), dung saprotrophs (Dung), litter saprotrophs (Litter), wood saprotrophs (Wood), mycorrhizal fungi (Myco), deciduous trees related fungi (Dec), conifers related fungi (Con).
Figure 5. The dynamics of (A) environmental paleoproxies and (B) ecophysiological groups of fungi over the last 10.5 kyr, modeled using GAM with formula = y ~ splines::bs(x, 5). Assigned ecophysiological roles of fungi abbreviations: plankton parasites (Plankton), animal pathogens (Animal), fish pathogens (Fish), insect parasites (Insect), fungal parasites (Fungal), plant pathogens (Plant), dung saprotrophs (Dung), litter saprotrophs (Litter), wood saprotrophs (Wood), mycorrhizal fungi (Myco), deciduous trees related fungi (Dec), conifers related fungi (Con).
Microorganisms 09 00719 g005
Figure 6. Associations between environmental paleoproxies and abundance in fungal ecophysiological groups over ~10.5 kyr. The two-way Pearson correlations are illustrated using Chord diagrams, where (A) shows strong positive correlations (r > 0.6) and (B) strong negative correlations (r < −0.6) over the whole timespan (~10.5 kyr); (C) shows strong positive correlations (r > 0.6) and (D) strong negative correlations (r < −0.6) over the period (2.0–10.5 kyr) excluding possible human impact (last ~2.0 kyr). Assigned abbreviations: water and terrestrial environment fungi (Aq_terr), terrestrial fungi (Terr), aquatic fungi (Aqua), pathotrophs (Patho), saprotrophs (Sapro), symbiotrophs (Symb), plankton parasites (Plankton), animal pathogens (Animal), fish pathogens (Fish), insect parasites (Insect), fungal parasites (Fungal), plant pathogens (Plant), mycorrhizal fungi (Myco), deciduous trees related fungi (Dec), conifers related fungi (Con).
Figure 6. Associations between environmental paleoproxies and abundance in fungal ecophysiological groups over ~10.5 kyr. The two-way Pearson correlations are illustrated using Chord diagrams, where (A) shows strong positive correlations (r > 0.6) and (B) strong negative correlations (r < −0.6) over the whole timespan (~10.5 kyr); (C) shows strong positive correlations (r > 0.6) and (D) strong negative correlations (r < −0.6) over the period (2.0–10.5 kyr) excluding possible human impact (last ~2.0 kyr). Assigned abbreviations: water and terrestrial environment fungi (Aq_terr), terrestrial fungi (Terr), aquatic fungi (Aqua), pathotrophs (Patho), saprotrophs (Sapro), symbiotrophs (Symb), plankton parasites (Plankton), animal pathogens (Animal), fish pathogens (Fish), insect parasites (Insect), fungal parasites (Fungal), plant pathogens (Plant), mycorrhizal fungi (Myco), deciduous trees related fungi (Dec), conifers related fungi (Con).
Microorganisms 09 00719 g006
Table 1. Community change indexes of fungal groups over time.
Table 1. Community change indexes of fungal groups over time.
Ecological GroupRate Change 1Synchrony 2Variance Ratio 3
Water environment0.040.24916.7
Terrestrial environment0.0130.0583.2
Both environments0.0170.1213.4
Pathotroph0.0380.16615
Animal pathogen0.0170.1754.9
Insects parasite0.006−0.1660.5
Fish pathogen0.003−0.2760.3
Plankton parasite0.0450.26917.8
Fungal parasite0.0090.0621.7
Plant pathogen0.0150.0772.7
Saprotroph0.0250.096.3
Litter saprotroph0.0040.031.1
Wood saprotroph0.0090.0461.8
Dung saprotroph0.004−0.0370.7
Symbiotroph0.0180.1345.2
Lichen symbiont0.0007−0.3260.1
Mycorrhizal fungi0.0470.14123.7
Conifer-related fungi0.001−0.0280.7
Deciduous tree-related fungi0.00020.0021.0
1 The rate change shows the rate of directional change in community composition over the timespan of the study. 2 The variance ratio is a positive number ranging from <1 (negative covariance) to >1 (positive covariance) centered at 1 (species vary independently). 3 The “gross” synchrony values range from −1 (complete asynchrony) to 1 (complete synchrony) with a center value of 0 (species fluctuate independently).
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Talas, L.; Stivrins, N.; Veski, S.; Tedersoo, L.; Kisand, V. Sedimentary Ancient DNA (sedaDNA) Reveals Fungal Diversity and Environmental Drivers of Community Changes throughout the Holocene in the Present Boreal Lake Lielais Svētiņu (Eastern Latvia). Microorganisms 2021, 9, 719. https://doi.org/10.3390/microorganisms9040719

AMA Style

Talas L, Stivrins N, Veski S, Tedersoo L, Kisand V. Sedimentary Ancient DNA (sedaDNA) Reveals Fungal Diversity and Environmental Drivers of Community Changes throughout the Holocene in the Present Boreal Lake Lielais Svētiņu (Eastern Latvia). Microorganisms. 2021; 9(4):719. https://doi.org/10.3390/microorganisms9040719

Chicago/Turabian Style

Talas, Liisi, Normunds Stivrins, Siim Veski, Leho Tedersoo, and Veljo Kisand. 2021. "Sedimentary Ancient DNA (sedaDNA) Reveals Fungal Diversity and Environmental Drivers of Community Changes throughout the Holocene in the Present Boreal Lake Lielais Svētiņu (Eastern Latvia)" Microorganisms 9, no. 4: 719. https://doi.org/10.3390/microorganisms9040719

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop