Marine eukaryote community responses to the climate and oceanographic changes in Storfjordrenna (southern Svalbard) over the past  ∼  13.3 kyr BP: insights from sedimentary ancient DNA analysis

Sedimentary ancient DNA (sedaDNA) metabarcoding is an emerging method for reconstructing the responses of marine organisms to past climate and oceanographic changes, including rare and non-fossilized taxa. Marine sedaDNA records from the Arctic are scarce, especially those focusing on the impact of environmental shifts on the biodiversity and functional composition of marine eukaryote communities. Here, we present a sedaDNA eukaryotic record from a sediment core retrieved in Storfjordrenna, southern Svalbard, spanning the termination of the Bølling–Allerød, the Younger Dryas, and the Holocene (13.3–1.3 kyr BP). We successfully recovered the eukaryotic communities and identified them by their ecological roles. Our study showed that the eukaryotic biodiversity in Storfjordrenna remained relatively stable, except during transitions between major climatic intervals. These shifts were characterized by changes in richness and relative abundance, driven by factors such as perennial ice cover, surface water cooling, and subsurface Atlantic water influx. Cercozoans and Marine Stramenopiles (MAST) emerged as dominant heterotrophs, characterized by high ecological flexibility and broad tolerance. Primary productivity was primarily driven by Arctic water (ArW) associated phytoplankton, including diatoms (Thalassiosira and Chaetoceros), green algae (Micromonas), and autotrophic dinoflagellates (Polarella glacialis) as well as the mixoplanktonic silicoflagellate Pseudopedinella elastica. Amplicon sequence variant (ASV)-based indicator analysis revealed that uncultured cercozoan lineages and MAST taxa were primarily associated with Atlantic water (AW) proxies, whereas parasitic dinoflagellates (Dino-group I) and choanoflagellates were more closely aligned with ArW proxies. Analysis of indicator responses shows the complex interactions within eukaryotic communities, and reveals a strong association among functional ecological groups that impact ecosystem productivity and regulation. This complexity highlights the limitations of traditional single-proxy approaches to accurately reconstructing paleoenvironmental conditions. Our study demonstrates the potential of high-resolution marine sedaDNA metabarcoding in elucidating responses to past climate changes, improving our understanding of the intricate interactions within eukaryotic communities, and enhancing our knowledge of marine ecosystems.

1 Introduction

The Arctic marine ecosystem is undergoing rapid and profound changes, primarily driven by climate warming (IPCC, 2023; Polyakov et al., 2017, 2020). A prominent feature of these changes is the increased influx of AW into the region, a phenomenon known as Atlantification. This process is associated with warming and rise of sea surface temperatures (SST), reduced sea-ice cover, altered salinity patterns, and changes in nutrient dynamics (Årthun et al., 2012; Polyakov et al., 2017).

These transformations in the marine environment are altering the biodiversity of the Arctic region and impacting the function and resilience of its ecosystems (Benner et al., 2019; Hallegraeff, 2010; Ribeiro et al., 2024). The Storfjordrenna region in southern Svalbard is an ideal location to study these shifts, having experienced significant climate-driven changes over the last ∼ 14 000 years, driven by meltwater discharge and the interaction between cold ArW and warmer AW inflows (Łącka et al., 2019, 2015; Pawłowska et al., 2020; Telesiński et al., 2024). The region's biodiversity has been shaped by, and remains sensitive to, these fluctuating and dynamic environmental conditions (Bensi et al., 2024; Deb and Bailey, 2023; Górska et al., 2022; Hop et al., 2019). Understanding how biodiversity adapts to such changes is essential for reconstructing past ecological responses to climate change and for predicting future trends. Although the impact of climate change on Arctic marine ecosystems is well documented (Deb and Bailey, 2023; Wassmann et al., 2010), relatively few studies have examined marine ecosystems using sedaDNA to analyze long-term biodiversity patterns (Grant et al., 2024; Pawłowska et al., 2020; Zimmermann et al., 2023). Recent developments in sedaDNA techniques have increased our ability to extract and analyze DNA from marine environments, providing valuable insights into eukaryotic communities and their responses to environmental changes over geological time scales (Grant et al., 2024; Harðardóttir et al., 2024; Zimmermann et al., 2021). Moreover, several studies have suggested that the sedaDNA or environmental DNA (eDNA) approach has the potential to use specific ASVs as proxies/bioindicators, even when their taxonomy is unknown thereby strengthening the connection between ASVs and ecological functions and environmental changes (Grant et al., 2024; Harðardóttir et al., 2024; Li et al., 2023; Lin et al., 2022; Pawłowska et al., 2020; Perret-Gentil et al., 2021, 2017; Zimmermann et al., 2021). Studies of marine eukaryotic sedaDNA have demonstrated that even low-resolution records can provide significant data on the shifts in marine communities over time, offering insight into past ecosystem dynamics (Grant et al., 2024). For example, recent studies have demonstrated that the sedaDNA approach can be used to reconstruct interactions between sea ice cover, ocean temperatures, and eukaryotic community composition (Armbrecht, 2020; Grant et al., 2024; Harðardóttir et al., 2024; Zimmermann et al., 2023, 2021). However, there is still a significant lack of suitable-resolution marine eukaryotic sedaDNA records from the Arctic, especially those focusing on the, ecosystem-oriented approaches that integrate ecological roles and biotic interactions to better link past biodiversity changes with ecosystem functioning and environmental changes.

This study addresses this issue by reconstructing the long-term history of marine eukaryotic communities from Storfjordrenna, southern Svalbard, using sedaDNA metabarcoding using V1V2 primers that capture both planktonic and benthic taxa with adequate taxonomic resolution for ecological interpretation (Fonseca et al., 2010; Lindeque et al., 2013; Sinniger et al., 2016). The sedaDNA record is supported by previously published sedimentological, micropaleontological, and geochemical records (Łącka et al., 2019, 2020, 2015; Telesiński et al., 2024). By focusing on eukaryotic communities associated with ArW and AW masses, we aim to assess their structure, ecological roles, and potential as indicators of past environmental conditions. Our approach seeks to identify how marine ecosystem have responded to significant climate-driven changes in this region and how these responses can improve our understanding of the future trajectories of Arctic marine biodiversity in the context of ongoing climate warming.

Figure 1 (a) Map of the study area and (b) location of core JM09-020GC (yellow dot). Red arrows indicate warm currents, and blue arrows indicate cold currents. Abbreviations: NAC: North Atlantic Current, WSC: West Spitsbergen Current, RAC: Return Atlantic Current, ESC: East Spitsbergen Current, BIC: Bear Island Current, EGC: East Greenland Current.

2 Study area

Storfjorden is an open fjord system located in the Svalbard archipelago, between the islands of Spitsbergen, Barentsøya, and Edgeøya (Fig. 1a). The cross-shelf through Storfjordrenna is located south of Storfjorden's mouth. The hydrography of Storfjorden and Storfjordrenna is primarily governed by the interplay of two major water masses: AW and ArW. AW is characterized by relatively high temperatures (> 3 °C) and high salinity (> 34.95 PSU), whereas ArW exhibits lower temperatures (< 0.5 °C) and salinity ranging from 34.3 to 34.8 PSU (Bensi et al., 2024; Skogseth et al., 2020; Sundfjord et al., 2017). AW is transported northwards by the Norwegian Atlantic Current, which bifurcates upon entering the Barents Sea into the West Spitsbergen Current and the North Cape Current (Blindheim and Østerhus, 2005). In contrast, ArW enters the region via the East Spitsbergen Current and the Bear Island Current, bringing cold, less saline waters into the Barents Sea (Hopkins, 1991). AW enters Storfjordrenna in a cyclonic manner, following the bathymetry. The Polar Front, which separates AW and ArW water masses, is located along the slope of Storfjordrenna (Bensi et al., 2024). The biological and geochemical signals preserved in the sediments represent an integrated response to both gradual and abrupt climatic and oceanographic changes (Łącka et al., 2019, 2015). Therefore, the site's proximity to the Arctic Front facilitates the detection of subtle ecosystem responses to long-term warming and changing oceanographic regimes.

3 Materials and Methods

3.1 Sediment core and age model

Gravity core JM09-020-GC was collected in 2009 at a depth of 253 m in Storfjordrenna, northwestern Barents Sea during the expedition of R/V Jan Mayen (Fig. 1b). The core was stored and processed according to the methods described by (Łącka et al., 2019, 2020, 2015). The chronology of the core was established based on AMS¹⁴C radiocarbon dating and one additional tie point defined by the appearance of vivianite micro-concretions in a sediment layer set to 12.8 kyr BP, coinciding with the onset of the Younger Dryas (Table S1 in the Supplement) (Łącka et al., 2020). The dates published first by Łącka et al. (2015), and Łącka et al. (2020) were recalibrated using the Marine20 calibration curve (Fig. S1 in the Supplement, Table S1) (Heaton et al., 2020). The palaeoceanographic history of Storfjordrenna over the past ∼ 14 000 years is well documented through detailed, multi-proxy reconstructions. These include analyses of fossil foraminifera assemblages, isotopic composition of foraminiferal tests, grainsize and elemental composition of sediments, alkenones (Łącka et al., 2019, 2015), and dinoflagellate cysts (Telesiński et al., 2024).

3.2 sedaDNA workflow

3.2.1 DNA extraction, amplification, and sequencing

Approximately 10 g of sediment was collected from 55 sediment layers using sterile spoons and transferred to sterile containers. DNA extractions were performed using the DNeasy PowerMax Soil Kit (Qiagen), following the manufacturer's instructions. All DNA extracts were stored at −20 °C until PCR amplification.

The V1V2 region of the 18S rDNA (with a length of ∼ 340 bp) was amplified by PCR using the forward primer SSU_FO4mod (5^′-GCT TGW CTC AAA GAT TAA GCC-3^′) and the reverse primer SSU_R22 (3^′-CCT GCT GCC TTC CTT RGA-5^′) (Fonseca et al., 2010; Lindeque et al., 2013), which were tagged with a unique 8-nucleotide sequence at their 5^′ ends (Esling et al., 2015). We followed the protocols established and positively tested by Pawłowska et al. (2020, 2014) and Lejzerowicz et al. (2013), with minor adjustments to the number of PCR cycles based on sample-specific amplification performance. Each sample was amplified in triplicate and each PCR reaction was performed in a total volume of 25 µL, which included 1.5 µL of 1.5 mM MgCl₂ (Applied Biosystems, USA), 2.5 µL of 10× PCR buffer II (Applied Biosystems), 0.5 µL of 0.2 mM deoxynucleotide triphosphates (Promega, USA), 0.5 µL of 20 mg mL⁻¹ bovine serum albumin (Invitrogen Ultrapure, USA), 1 µL of 10 µM of each primer, 0.2 µL of AmpliTaq Gold DNA polymerase (Applied Biosystems) and 2 µL of template DNA. The amplification conditions consisted of a pre-denaturation step at 95 °C for 5 min, followed by 50 cycles of denaturation at 95 °C for 30 s, annealing at 57 °C for 30 s and extension at 72 °C for 1 min, followed by a final extension step at 72 °C for 5 min. PCR products, including negative control for each unique combination of tag-encoded primers, were verified by agarose gel electrophoresis. PCR products were purified using the High Pure PCR Cleanup Micro Kit (Roche) and quantified using a Qubit 2.0 fluorometer. Libraries were pooled in equimolar quantities and the sequence library was prepared using a TruSeq library-preparation kit (Illumina). Samples were then loaded into a MiSeq instrument for a paired-end run of 2 × 250 cycles. The sequencing was performed at the University of Geneva.

3.2.2 Data quality control and processing

The raw sequencing reads for each sample were processed using the SLIM web application (Dufresne et al., 2019). In brief, the module demultiplexer was used to demultiplex the raw reads according to their unique tags in the forward and reverse reads. Quality filtering, chimera removal and ASVs table generation were performed using DADA2 v.1.16 with pseudo-pool parameters (Callahan et al., 2015).

The ASVs were then curated using the LULU package v.0.1.0 (Froslev et al., 2017) with the default parameters. The taxa assignment of the ASVs was performed using VSEARCH against the taxonomically curated PR2 database v.4.14.1 (Guillou et al., 2012), which contains functional annotations. We used a Last Common Ancestor approach, assigning to the consensual taxonomic rank to up to reference sequences with at least 80 % similarity as a threshold for the dataset. The ASVs were also assigned to functional groups with at least 95 % similarity; with the functional attributes of Ibarbalz et al. (2019). The ASVs assigned to prokaryotes (bacteria and archaea) were removed in order to analyze only eukaryotic ASVs. Additionally, fungi and gymnamoebae were removed due to the high risk of contamination (Armbrecht, 2020). Unique ASVs (occurring in only one sample), short sequences (< 200 bp), rare ASVs (having < 100 reads), and low read count samples (< 1000 reads) were removed from the dataset. Additionally, the unassigned sequences were blasted with NCBI to further clarify the taxonomic composition. The Cumulative Sum Scaling (CSS) technique (scale factor – 0.9) was used to transform the read counts in the dataset and using `cssNorm' function in the genomeSeq package (Paulson et al., 2013), and used for the downstream statistical analysis in the study.

3.3 Statistical analysis

Data analysis was performed in R v.4.5.2 (Team, 2025) using several R packages. The relative abundances of reads and ASVs for each eukaryote group were calculated and plotted using ggplot2 (Wickham, 2016) and Grapher 24.1.213.

Three alpha diversity indices were calculated for all samples (ASV richness q=0, Shannon index q=1, and Simpson index q=2) based on the functions in the vegan package (Oksanen et al., 2025). Shannon diversity was compared between the main groups using the Kruskal–Walli's rank test in the stats package (Kruskal and Wallis, 1952). Significance between the groups was determined using a pairwise Wilcoxon rank sum test with an adjusted p value (Benjamini–Hochberg) in the ggpubr package (Kassambara, 2026). A Principal Coordinate Analysis (PCoA) ordination was generated using the Bray–Curtis dissimilarity matrix calculated using the ape package (Paradis and Schliep, 2019) and stringer package (Wickham, 2025) to assess beta diversity and visualize dissimilarities in eukaryotic community composition among the samples. Permutational Multivariate Analysis of Variance (PERMANOVA) with 999 permutations via the adonis2 function in the vegan package (Oksanen et al., 2025).

A co-occurrence heatmap representing most of the families in the study was generated. The corrplot and pheatmap packages (Kolde, 2025; Wei and Simko, 2024) were used to analyze the correlation between eukaryote families based on spearman method (a cutoff mark in correlation as > 0.5, and p value Benjamini–Hochberg (BH) adjusted < 0.05). Environmental parameters represented by palaeoceanographic proxies were used to identify the response of eukaryote species using three analytical methods. Seven proxies were used, including indicators of: (i) sea surface temperature (SST ${\mathrm{U}}_{\mathrm{37}}^{\mathrm{K}\ast }$) Łącka et al. (2019), (ii) AW (foraminifera Nonionellina labradorica and Buccella frigida Łącka et al. (2015), and dinocyst Operculodinium centrocarpum (Telesiński et al., 2024), (iii) ArW/meltwater (% C37:4 Łącka et al. (2019), (iv) glaciomarine condition (foraminifera Elphidium excavatum and Cassidulina reniforme Łącka et al. (2015), (v) sea ice (dinocyst Echinidinium karaense Telesiński et al. (2024), and (vi) bottom current dynamics (mean grain size 0–63 µm) Łącka et al. (2015). Multidimensional Fuzzy Set Ordination (MFSO) correlation plot, and Fuzzy set ordination (FSO) plots were generated for each environmental variable to assess their influence on eukaryote communities and identify key proxies for downstream analysis (Roberts, 2008). Firstly, a heatmap of sparse partial least squares (sPLS) regression between ASVs, and proxies was generated using the spls and cim function in mixOmics package (Froslev et al., 2017; Rohart et al., 2017). The potential ASV based indicators were selected based on a correlation coefficient threshold of > 0.3 and BH adjusted p value of < 0.05. Secondly, a Spearman correlation heatmap of the top 100 most significant ASVs (rho > 0.3 and p-adjust < 0.05) as generated using the pheatmap package. Finally, the dataset was analyzed using DEseq2 analysis, and the dataset was curated based on log2FoldC ≥ 1 (BH-adjusted p value < 0.05), and lower base mean (Love et al., 2014). Potential indicator ASVs were categorized based on their correlation strength and consistent detection across at least two methods or strong association with multiple paleo-proxies.

Figure 2 PCoA based on the Bray–Curtis dissimilarity matrix method with the eukaryote dataset (raw data converted into CSS formation).

4 Results

4.1 Metabarcoding data

A total of 2 620 808 raw sequence reads were generated from 55 samples. After initial quality filtering and denoising with DADA2 in SLIM, eight samples were excluded due to low number of reads (< 300). An additional five samples were removed during the downstream eukaryote-specific quality control process described in the Methods. This reduced our dataset to 1 609 500 sequence reads and 273 ASVs in 42 samples (Tables S2 and S3 in the Supplement).

4.2 Alpha diversity

Alpha diversity indices varied across time interval; however, Shannon diversity showed no significant differences among intervals (Kruskal–Wallis, p = 0.48; Fig. S2 in the Supplement, Table S2). This apparent discrepancy reflects that species richness is sensitive to rare taxa, whereas Shannon diversity accounts for both richness and evenness, making it less responsive to occasional low-abundance ASVs. Overall, the number of observed ASVs ranged from 4 to 144 (Fig. S2, Table S2). The highest values were observed in the Younger Dryas, the Early and Late Holocene, particularly at 12.3, 11.3, 9.5, 4.0, and 2.8–2.3 kyr BP. In contrast, a significant decrease in richness was observed around 13.3, 11.8, 2.17, and 1.8 kyr BP (Fig. S2, Table S2). Similar trends were revealed by both the Shannon and Simpson indices, with minimal diversity observed around 12.8, 11.7, 9.2 kyr, and 3.4 kyr BP. Due to limited data resolution, no clear trends in alpha diversity could be discerned between 9.2 kyr BP and 3.4 kyr BP. However, a decline was evident after 3.4 kyr BP, continuing towards 1.3 kyr BP (Fig. S2).

4.3 Beta Diversity

Beta diversity analyses revealed minor changes in community composition during the transitions from the Bølling–Allerød to the Younger Dryas, and from the Younger Dryas to the Holocene (Fig. 2). The PCoA plot revealed partial overlap across different time intervals, with Early Holocene and Younger Dryas samples were widely dispersed, while Late Holocene samples clustered separately along Axis 2, respectively. In contrast, Bølling–Allerød and Mid Holocene samples were largely scattered along the Axis 2 (Fig. 2). PERMANOVA results confirmed a significant effect of the time intervals on eukaryote community (R² = 0.16, p < 0.05), with significant pairwise differences detected between the Late and Early Holocene, Late Holocene and Younger Dryas, Early Holocene and Younger Dryas, and Younger Dryas and Bølling–Allerød (Table S4 in the Supplement). FSO plots revealed a significant relationship between the samples and several paleo-environmental proxies, including the dinocyst Operculodinium centrocarpum from (Telesiński et al., 2024), ArW/meltwater indicator (% C37:4) (Łącka et al., 2019), and sea surface temperature (SST ${\mathrm{U}}_{\mathrm{37}}^{\mathrm{K}\ast }$) (Łącka et al., 2019) (Fig. S3 in the Supplement). In the MFSO framework, these proxies together explained a moderate proportion of the variation in community composition (cumulative fuzzy set correlation r = 0.454) (Fig. S3).

Figure 3 (a) Bar plot showing the downcore distribution of eukaryotic sub-divisions based on their relative abundance. (b) Proportional richness of distinct ecological groups across selected time periods (Bølling–Allerød, Younger Dryas, Early Holocene, Mid Holocene, and Late Holocene), expressed as the percentage of ASVs. The accompanying table provides the number of samples, the total number of observed ASVs and unique ASVs within each climate time interval.

4.4 Community composition

Within the dataset, a total of 236 ASVs were assigned to 14 Sub-Divisions, while 37 ASVs remained unassigned (Table S3 in the Supplement). The Cercozoa was the most abundant sub-division, comprising 67 ASVs, accounting for 24.54 % of the total ASVs (Fig. 3, Table S3). Overall, the taxonomic structure of eukaryotes based on read abundance fluctuated significantly between samples (Fig. 3a). In contrast, ASV richness remained stable across different time periods except during the Bølling–Allerød period (Fig. 3b). The number of unique ASVs was highest during the Late Holocene, with five unique ASVs identified. The Younger Dryas and Early Holocene each exhibited three unique ASVs. Conversely, no distinctive ASVs were identified during the Mid Holocene and Bølling–Allerød periods. Across the entire Holocene, a total of 52 unique ASVs were recorded (Fig. 3b).

The ASVs were categorized based on their ecological roles, such as phytoplankton (31 ASVs), mixoplankton (7 ASVs), mesozooplankton (8 ASVs), planktonic heterotrophic protists (23 ASVs), parasites (24 ASVs), zoobenthos (36 ASVs), nekton (2 ASVs), and other heterotrophs (104 ASVs) (Fig. 3b, Table S3). The latter category comprises multiple taxonomic groups characterized by complex habitats and feeding behaviors, many of which have poorly understood ecological roles. This group includes Cercozoa (67 ASVs), Labyrinthulea (11 ASVs), Choanoflagellatea (5 ASVs), MAST (15 ASVs), benthic ciliophora (2 ASVs), Apusomonada (3 ASVs), and one ASV from the Telonemea flagellate group. The unassigned taxa (37 ASVs) also remained ecologically uncategorized (Fig. 3b).

Figure 4 Relative abundance and richness (expressed as observed ASVs percentage) of major ecological groups, along with Shannon index, richness, and sea surface temperature (SST UK'37) from Łącka et al. (2019). Lines represented ASVs abundance (%), and area represented the read abundance (%).

The phytoplankton community consisted of diatoms, green algae, haptophytes, and autotrophic dinoflagellates, most of which were associated with ArW (Fig. 4, Table S3). In terms of read abundance, Thalassiosira spp. and Chaetoceros sp. dominated among diatoms, while Micromonas polaris was the dominant species within the green algae. The haptophytes group was primarily represented by Phaeocytis sp., whereas the Gymnodinium spp. and sea-ice-associated species Polarella glacialis were dominant within the autotrophic dinoflagellate group (Fig. S4 in the Supplement). The mixoplankton community was primarily composed of mixotrophic dinoflagellates and silicoflagellates. In terms of read abundance, mixotrophic dinoflagellates were mainly present from the Younger Dryas to the beginning of the Early Holocene, whereas mixotrophic silicoflagellates, represented by Pseudopedinella sp., was present throughout the entire core (Fig. S4) (Sects. S1 and S2 in the Supplement).

The planktonic heterotrophic protists group included radiolarians, pelagic ciliates, dinoflagellates, and silicoflagellates. These groups were identified as being present at specific time periods, e.g. ∼ 12.4 to ∼ 10.2 kyr, and ∼ 2.3 to 1.3 kyr BP (Fig. S5 in the Supplement, Table S3). The mesozooplankton group comprises small metazoans and was dominated by arthropods (Copepoda and Malacostraca) and larvaceans (Appendicularia). The copepod Calanus spp. represented the majority of the mesozooplankton around the study area (Fig. S5) (Sects. S3 and S4 in the Supplement).

The zoobenthos was recorded as the most diverse group, primarily representing macrobenthic species (Fig. 4). This group included annelids, ascidiacean, molluscs, cnidarians, and echinoderms (Fig. S6 in the Supplement, Table S3). Zoobenthos taxa were most abundant around ∼ 9.3 kyr BP, as well as between ∼ 12.3 kyr BP (Fig. 4) (Sect. S5 in the Supplement).

The parasites were represented by 6 classes: Dinophyceae (Syndiniales), Gregarinomorphea, Paragregarea, Peronosporea, Hyphochytrea, and Enoplea (Fig. 4, Table S3). Among them, Syndiniales have the highest abundance and diversity, with 18 ASVs (mainly uncultured) detected throughout the studied time interval (Fig. S7 in the Supplement) (Sect. S6 in the Supplement). The nekton group included two ASVs assigned to Arctic cod (order Gnathostomata), which were detected only during the Younger Dryas and Early Holocene.

Among other heterotrophs, Cercozoa were dominant, accounting for a significant proportion of reads throughout the study period (Fig. 4, Table S3). Five classes of Cercozoa were identified: Ascetosporea, Phytomyxea, Granofilosea, Thecofilosea, and Imbricatea. Thecofilosea exhibited the highest richness with 51 ASVs (Fig. S8 in the Supplement, Table S3). Ecologically, cercozoans can be classified as parasitic, predatory, and bacterivorous. Cryothecomonas spp. were identified as predatory, while those in the classes Ascetosporea and Phytomyxea were classified as parasites. Other cercozoan ASVs were identified only to the family level, limiting precise ecological interpretations.

The other heterotrophs also include the MAST (Marine Stramenopiles), Labyrinthulea, Choanoflagellata and Apusomonada. The MAST group included 16 ASVs, representing four main sub-clades: MAST-1, MAST-3, MAST-9, and MAST-12. MAST-9 and MAST-12 dominated throughout the studied time period and exhibited high richness (Fig. S7, Table S3).

The Labyrinthulea included saprotrophic Thraustochytriaceae, and Aplanochytriidae, revealing a dominant presence and high richness around ∼ 11.8 to ∼ 2.2 kyr BP, and ∼ 1.7 kyr BP (Fig. 4). Most choanoflagellate ASVs belonged to environmental clades, except for Calliacantha sp., which was dominated around ∼ 12.7 kyr BP to ∼ 12.2 kyr BP (Fig. 4). Apusomonada, represented by the class Apusomonadea, appeared during certain time intervals, especially the Early Holocene (∼ 10.4 kyr BP), and Late Holocene (∼ 1.7 kyr BP) (Fig. 4, Table S3) (Sect. S7 in the Supplement).

4.5 Indicator taxa for Arctic and Atlantic water conditions

A total of 44 ASVs were identified as potential indicator taxa using three analytical approaches (sPLS, Spearman correlation, and DESeq) (Figs. S9 and S10 in the Supplement, Table S5 in the Supplement). Of these, 16 ASVs were identified as potential AW indicators, and belonged to the following groups: phytoplankton (1), planktonic heterotrophic protists (1), cercozoans (6), MAST (3), zoobenthos (1), labyrinthulean (1), and unassigned ASVs (1). The AW plankton indicators included a green algae Pyramimonas sp. (ASV55), heterotrophs-Silicoflagellate Pteridomonas sp. (ASV105), and a pelagic ciliate Cyclotrichium sp. (ASV265), while the AW benthic indicators comprised a polychaete Tharyx sp. (ASV278) (Fig. S11 in the Supplement, Table S5).

In contrast, 28 ASVs were associated with ArW (Table S5), primarily parasites (6), zoobenthos (3), phytoplankton (4), choanoflagellates (1), planktonic heterotrophic protists (3), mesozooplankton (1), mixoplankton (1), benthic ciliate (1), nekton (1), cercozoans (1), and unassigned ASVs (6) (Fig. S12 in the Supplement, Table S5). Potential ArW ASV-based indicators were identified among both planktonic and benthic taxa. Planktonic ASVs comprised autotrophic dinoflagellates (Prorocentrum sp., ASV64 and ASV153; Gymnodinium sp., ASV29), diatoms (Chaetoceros gelidus, ASV18), mixotrophic dinoflagellate (Heterocapsa rotundata, ASV204), pelagic ciliate (Strombidium sp., ASV123), radiolarian (Heteracon sp., ASV54; Acanthoplegma sp., ASV677), and mesozooplankton (Oikopleura sp., ASV500). Parasites included three Syndiniales dinoflagellates (ASV213, ASV341, and ASV391), an apicomplexan (Paralecudina sp., ASV958), and a parasitic nematode (Mermithidae sp., ASV101). Benthic indicators included benthic ciliate (Holosticha sp., ASV87), echinoderm (Ctenodiscus sp., ASV104), and bivalve (Tridonta sp., ASV137). Additionally, Calliacantha natans (ASV92) was identified as an ArW indicator in our study (Table S5).

Figure 5 Potential ASV-based indicators of AW and ArW conditions from the study (correlation > 0.4, p < 0.05, positive results from at least two method (displayed data as relative abundance %, and included the ASVs recorded in more than four samples).

Eight ASVs were identified as the most robust ASV-based indicators of AW and ArW conditions, based on their consistent occurrence through time (detected in more than four samples) and strong statistical support (correlation coefficient > 0.4; p < 0.05) (Fig. 5).

Figure 6 Heatmap of co-occurrence based on spearman rank coefficient analysis between eukaryote families represent in the study (illustrated only the positive correlation ≥ 0.4, and p value adjusted (BH) < 0.05).

4.6 Ecological interactions among eukaryote families

Spearman correlation analysis (r > 0.5, adjusted p < 0.001) was used to explore potential ecological interactions among eukaryotic families inhabiting similar environmental niches (Fig. 6, Table S6 in the Supplement). Parasitic cercozoans (Ascetosporea and Phagomyxidae) strongly correlated with algae families (Phaeocystaceae, Thalassiosiraceae, Pyramimonadaceae, and Prasinodermataceae), and dinoflagellates (Suessiaceae, and Gymnodiniaceae). Other cercozoans (CCW10-lineage, Novel-Clade-2, Cryothecomonas-lineage, and Ventricleftida) revealed significant correlations with algal and dinoflagellate groups (Fig. 6, Table S6). Among MAST groups, MAST-12 showed positive associations with algae families (Prasinodermataceae, Thalassiosiraceae, and Chaetocerotaceae) and parasites (Pirsoniaceae), while MAST-9 correlated with multiple phytoplankton families (Fig. 6, Table S6). Parasitic dinoflagellates (dino-group-II) showed strong correlation with haptophyte algae (Phaeocystaceae), diatom (Thalassiosiraceae), and dinoflagellates (Suessiaceae, and Gymnodiniaceae). Parasitic alveolates of the family Lecudinidae showed a strong positive correlation with the mesozooplankton (Malacostraca), radiolarians and ophiuroids (Fig. 6, Table S6). Another parasitic superfamily of Stramenopiles, the Pirsoniaceae displayed strong associations with various taxa, including MAST-12, cercozoans (Thaumatomonadidae), green algae (Prasinodermataceae, and Chlamydomonadales), haptophytes (Phaeocystaceae), diatoms (Thalassiosiraceae), silicoflagellates (Actinomonadaceae), dinoflagellates (Amphidiniopsidaceae), and polychaetae (Chaetopteridae) (Fig. 6, Table S6).

5 Discussion

This study expands our knowledge of eukaryotic communities' patterns in the Storfjordrenna over the past 13 300 years by providing high-resolution sedaDNA records of both fossilized and non-fossilized groups. We demonstrate how these communities responded to major climatic shifts since the Bølling–Allerød and highlight key ecological interactions among major taxonomic groups. These findings enhance our understanding of how environmental changes have shaped eukaryotic biodiversity in southern Svalbard.

5.1 Impacts of oceanographic changes on the eukaryotic community in Storfjordrenna

5.1.1 Bølling–Allerød (13.30 to 12.80 kyr BP)

The eukaryotic community in our record during the Bølling–Allerød reflected oceanographic conditions resembling those observed today in glacier-proximal areas of Arctic fjords, characterized by high turbidity due to meltwater discharge and the presence of colder, fresher waters (Łącka et al., 2019; Zajączkowski, 2008). Previous studies support this interpretation: the grounding line of the Svalbard Barents Ice Sheet (SBIS) retreated from Storfjordrenna before 13.95 kyr BP (Łącka et al., 2015), coinciding with SST reaching modern-like values (Łącka et al., 2019). However, despite elevated SST, primary productivity remained low likely due to the suppressive effect of turbid meltwater input from the retreating ice sheet (Łącka et al., 2015). Furthermore, biomarker data indicated a dominance of fresher ArW over AW, which has been linked to reduced primary productivity (Łącka et al., 2019). These conditions favored the development of a eukaryotic community dominated by heterotrophs, capable of thriving in such extreme environments. The most abundant taxa, in terms of both sequence reads and ASV richness, were bacterivorous cercozoans (Fig. 4). The dominant cercozoan was Limnofila sp., a genus primarily found in fresh and brackish waters (Mylnikov et al., 2015; Nikolaev et al., 2003), the presence of Limnofila sp. may reflect either local ecological conditions or allochthonous input via riverine transport or ice-rafted debris (Andruszkiewicz et al., 2019; Jo et al., 2025; Nguyen et al., 2026). Other important bacterivorous heterotrophs were MAST, particularly MAST-9D and MAST-12A, which are known to be adapted to extreme environmental conditions (Labarre et al., 2021; Lin et al., 2022; Obiol et al., 2024).

Despite lower read abundance, phytoplankton, mixoplankton, and planktonic heterotrophic protists showed high ASV richness during this period (Fig. 4). The phytoplankton community was dominated by autotrophic, sea-ice associated taxa, such as dinoflagellates Polarella glacialis (Harðardóttir et al., 2024) and Gymnodinium spp. (Kubiszyn and Wiktor, 2016), and diatom Thalassiosira spp. (Luddington et al., 2016). The mixoplankton was represented by the silicoflagellate, Pseudopedinella elastica (Fig. S4), which has been described as bacterivorous under conditions of limited light and nutrients (Gerea et al., 2016). The mesozooplankton community was primarily composed of the herbivorous Calanus spp. and the omnivorous Metridia longa. These species have been previously observed in the Svalbard region, with Calanus spp. dominating in terms of biomass (Daase et al., 2008). Despite unfavorable conditions caused by meltwater influx and low nutrient availability, both primary and secondary productivity persisted, likely concentrated in ice-proximal and frontal zones where the stratification enhanced nutrient retention and water column stability (Łącka et al., 2019, 2015).This suggests that environments near retreating ice sheets can act as biological hotspots, supporting productivity through ice-associated blooms and complex interactions within the microbial food web.

5.1.2 Younger Dryas (12.80 to 11.70 kyr BP)

The eukaryotic community during the Younger Dryas reflected the dramatic environmental changes that occurred at that time. The most notable change was the rapid decrease in biodiversity during the Bølling–Allerød and Younger Dryas transition, when alpha diversity indices reached near-zero values (Fig. 4). The reorganization of oceanographic conditions most likely caused a temporary slowdown of Atlantic meridional overturning circulation (AMOC) and a reduction in AW inflow (Łącka et al., 2020). This led to strong stratification, formation of perennial ice cover, and anoxic conditions at the bottom (Łącka et al., 2020). The presence of perennial ice cover led to a significant reduction in primary productivity in Storfjordrenna (Łącka et al., 2019). However, the sedaDNA record also revealed the presence of phyto- and mixoplankton during this period, especially presence of phytoplankton Thalassiosira spp. and Gymnodinium sp., and silicoflagellate P. elastica. This may suggest that, although limited, primary productivity still occurred under the ice. The detection of herbivorous mesozooplankton Calanus spp., and predatory Cryothecomonas spp., also coincided with the presence of phytoplankton. Notably, the early Younger Dryas also revealed a short-term increase in the relative abundance and diversity of zoobenthos, primarily polychaetae (Barantolla sp.) and molluscs (Tridonta sp. and Talochlamys sp.) (Fig. S6).

During the latter part of the Younger Dryas (after ∼ 12.4 kyr BP), increasing advection of AW and SST warming led to the replacement of perennial sea ice by seasonal ice cover (Łącka et al., 2019, 2020). This shift was followed by the development of a more diverse benthic foraminifera community (Łącka et al., 2015). Similarly, the sedaDNA record displayed the increase in the richness and abundance of zoobenthic taxa, mainly annelids, molluscs, and echinoderms (Fig. S6). However, the alkenone record suggested that the warming was may associated with low primary productivity, probably due to the continuous input of turbid meltwater from the decaying SBIS (Łącka et al., 2015).

In contrast, the sedaDNA record indicated a sudden phytoplankton bloom in the late Younger Dryas, dominated by M. polaris, Thalassiosira spp., Chaetoceros gelidus, and Gymnodinium spp. Micromonas polaris is typically associated with Arctic sea-ice environments (Bachy et al., 2022), and was recorded during the periods of low SST. Chaetoceros gelidus is known for its high tolerance under variable light and ocean acidification conditions (Biswas, 2022; Ribeiro et al., 2024) and may play a key role in plankton blooms and primary productivity, particularly during the Younger Dryas (Fig. S4). Phytoplankton blooms stimulated the development of secondary producers, mainly pelagic ciliates, and radiolarians, as well as mesozooplankton copepods (Calanus spp.) (Fig. S5). Altogether, these findings indicate that the latter part of the Younger Dryas (after ∼ 12.4 kyr BP) was characterized by periods of accelerated AW inflow and higher SST (Risebrobakken et al., 2010; Wollenburg et al., 2004), which promoted phytoplankton and zooplankton growth, and enhanced development of the benthic community.

5.1.3 Early Holocene (11.70 to 9.20 kyr BP)

The transition from the Younger Dryas to the Early Holocene was characterized by a significant decrease in biodiversity and the dominance of mixoplankton, primarily silicoflagellate P. elastica, and dinoflagellates such as Biecheleria sp. and Gotoius sp. According to Łącka et al. (2020), the onset of the Early Holocene was associated with a short-term decrease in SST and a decrease in foraminiferal fauna abundance. The low biodiversity and dominance of mixotrophic plankton observed in the sedaDNA record might be a consequence of this short-term deterioration in environmental conditions.

However, the further development of the Early Holocene was driven by an increasing influence of AW in the area, which was followed by an increase in SST and productivity (Devendra et al., 2023; Telesiński et al., 2018). Moreover, Arctic Front was located close to the Spitsbergen coast, leading to the formation of a highly productive frontal zone (Łącka et al., 2019). The amelioration of environmental conditions during the Early Holocene (Łącka et al., 2015) was reflected in a sudden peak of alpha diversity of overall eukaryotic community, accompanied by a notable increase in the richness and abundance of key ecological groups, including phytoplankton, zoobenthos, parasites, and other heterotrophs such as cercozoans and MAST (Fig. 4). However, taxa associated with sea-ice were an important component of the assemblage, suggesting that sea-ice formation still occurred in Storfjordrenna. Despite AW dominance, the presence of a brackish water cercozoan Limnofila sp., green algae M. polaris and Pyramimonas sp., as well as sea-ice-indicator P. glacialis may suggest episodic presence of sea-ice, and transition of Arctic Front (Figs. S4 and S8). The overall high eukaryotic biodiversity in the Early Holocene, particularly the diversity of phytoplankton, mixoplankton, mesozooplankton, and the gradual increase in MAST-9 species related to warm water further support the establishment of warm-water conditions with high nutrient availability (Łącka et al., 2019) (Figs. S4, S5, S7, and Table S5).

5.1.4 Mid Holocene (9.20 to 3.40 kyr BP)

Samples from the Mid Holocene, spanning the period between 4.0 and 7.5 kyr BP, show low eukaryote DNA recovery, likely due to extensive degradation, and are dominated by fungi and amoebozoan DNA, particularly Acanthamoeba, which constitutes a large proportion of the recovered sequences. The dominance of amoebozoan DNA may reflect elevated microbial activity, potentially accelerating post-depositional DNA degradation (Anderson, 2017; Boere et al., 2011; Butler and Rogerson, 1995). As a result, lack of eukaryotic sedaDNA data from the period between ∼ 7.5 and 4.0 kyr BP (see Results), restricting interpretation to the early Mid Holocene and the transition toward the Late Holocene. Therefore, this interpretation should be treated with caution due to the limited number of samples analyzed in this interval. The beginning of the Mid Holocene (9.2 kyr BP) in Storfjordrenna was marked by a significant drop in biodiversity, followed by an increase after 8.0 kyr BP. The species composition was predominantly composed of cercozoans, mainly Limnofila sp. and Cryothecomonas spp. (Fig. S8). Another important component of the eukaryotic assemblage was MAST species, including MAST-9D and MAST-12B, which had previously been recorded in the north Atlantic region (Lopez-Garcia et al., 2007; Newbold et al., 2012) (Fig. S7). The community composition resembled the one from Bølling–Allerød, dominated by heterotrophic taxa adapted to unfavorable environmental conditions. This aligns with evidence of a minor cooling event between 9.0 kyr BP and 8.0 kyr BP, as proposed in previous studies (Łącka et al., 2015).

In contrast, despite their relatively low abundance, the communities of phytoplankton, mesozooplankton, and planktonic heterotrophic protists displayed relatively high diversity in the early Mid-Holocene. This period in Storfjordrenna was characterized by limited ice rafting, variable SST, and interplay between the AW and ArW water masses rather than a continuous impact of AW (Łącka et al., 2019, 2015). Furthermore, the low alkenone flux indicated low primary productivity throughout the mid-Holocene (Łącka et al., 2019), which is consistent with the reduced representation of both phyto- and zooplankton taxa in the sedaDNA record and supports an interpretation of sustained low productivity–related community signals during this interval. Low productivity was also observed at that time in the Norwegian and Svalbard shelves, potentially due to the limited nutrient availability. According to Łącka et al. (2019), the reduction in primary productivity resulted from enhanced vertical stratification, which reduced vertical mixing in the water column, and thus, limited the nutrient resuspension to the surface. An alternative explanation is the early spring bloom, that occurs in the ice-free waters, and the subsequent development of mesozooplankton that graze on phytoplankton, thereby reducing the flux of organic matter to the bottom. However, the presence of low abundance sequences assigned to both phytoplankton and planktonic heterotrophic protists, as well as the increase in bacterivorous taxa, likely supports the first scenario. Overall, the lack of sea ice, and the variability in water masses and SST observed at the beginning of the Mid Holocene, created an unstable environment, which favored tolerant heterotrophic eukaryotes such as cercozoans or MAST.

5.1.5 Late Holocene (3.40 to 1.30 kyr BP)

The onset of the Late Holocene was marked by an increase in eukaryotic biodiversity, followed by a sharp decrease around 2.0 kyr BP. During this period, eukaryotic communities were predominantly composed of cercozoan and MAST (Fig. 4). Cercozoan abundance and richness exhibited an increasing, yet variable trend throughout the Late Holocene, whereas MAST decreased progressively over time.

Both phytoplankton and planktonic heterotrophic protists exhibited high richness, but variable abundance throughout the Late Holocene (Fig. 4). Furthermore, the presence of parasitic species, including the Syndiniales dinoflagellate (dino-group-I and dino-group-II) as well as the diatom-associated parasitic Pirsonia sp. (Schweikert and Schnepf, 1997), co-occurred with the phytoplankton suggesting that parasitic interactions may have influenced phytoplankton dynamics during the Late Holocene. (Fig. S7, Table S6). The Late Holocene coincided with the so-called Neoglacial cooling, which spanned the last 4.0 kyr BP. This period was characterized by a decline in SST (Risebrobakken et al., 2010), limited AW inflow and strengthening of ArW flow, which led to the formation of extensive ice cover (Berben et al., 2014; Devendra et al., 2023; Martrat et al., 2003). Records from Storfjordrenna also showed a cooling in the area, associated with enhanced ice rafting (Łącka et al., 2019, 2015). Thus, the increased abundance of phytoplankton in general, and ice-associated species P. glacialis in particular, is probably an effect of the cooling of surface waters and the formation of sea-ice, which launched convective water mixing and nutrient resuspension to the surface. In consequence, primary productivity increased, stimulating the development of the planktonic heterotrophic protists community (Fig. S5).

5.2 sedaDNA environmental indicators

This study identified 46 potential eukaryotic indicator taxa associated with AW and ArW conditions. Several of these taxa exhibited consistent temporal patterns that aligned with paleoenvironmental proxies, supporting their potential for long-term reconstructions. In contrast, others appeared only sporadically, which reduces their reliability as potential indicators. AW-associated taxa were primarily represented by cercozoans and MAST, while ArW-associated taxa included diatoms, dinoflagellates, choanoflagellates, Arctic zoobenthos, and zooplankton.

Bacterivorous cercozoans, including the Ventricleftida (ASV46), the Protaspa-lineage (ASV83, and ASV257), Ventricleftida spp. (ASV46) and the Cryothecomonas spp. (ASV377), were identified as potential AW indicators (Fig. S11, Table S5). However, their identification is currently based exclusively on molecular data, limiting ecological and biogeographical context and weakening their use in environmental reconstructions (Labarre et al., 2021; Obiol et al., 2024). Similarly, members of the bacterivorous MAST-9 group, notably MAST-9A and MAST-9D, were exclusively detected in AW conditions (Fig. 5, Table S5), which is consistent with their known tropical-to-temperate distribution. Within the phytoplankton communities, Pyramimonas parkeae (a green alga that prefers higher temperature regions; Bock et al., 2021), and the planktonic heterotrophic pelagic ciliate Cyclotrichium sp. (commonly found in warmer waters; Dirmenci et al., 2010; Xu et al., 2005), were also identified as potential AW indicators. However, they only occurred in a brief temporal window near the end of the Early Holocene, so their reliability as indicators needs to be verified by further studies (Fig. S11).

The relatively high number of cold-water species recorded was probably due to favorable overall conditions in the study area. The autotrophic Prorocentrum spp. are known to be toxin-producing, bloom-forming species with broad global distributions, including polar regions (Cen et al., 2020; Goncharenko et al., 2021; Stoecker and Lavrentyev, 2018; Tillmann et al., 2022). In the present study, both taxa were primarily detected at the onset of the Younger Dryas (Table S6). Their limited distribution suggests that they may be unreliable as long-term environmental indicators. Similarly, Heterocapsa rotundata a mixotrophic dinoflagellate commonly associated with harmful algal blooms in Arctic and North Atlantic waters (Rintala et al., 2010; Wu et al., 2022), was identified as a potential ArW indicator (Fig. S12). The genus Holosticha is a widespread benthic ciliate, associated with sea ice (Berger, 2003; Petz et al., 1995; Wilbert and Song, 2008), and known to feed on diatoms and flagellates (Lei et al., 2005), which was also identified as a potential ArW indicator (Fig. S12). However, their presence was confined to the Bølling–Allerød and Younger Dryas or the Younger Dryas and Early Holocene intervals, respectively, limiting their reliability as a long-term proxy for ArW conditions (Fig. S12).

In contrast, the ArW-associated diatom species Chaetoceros gelidus was consistently abundant (Fig. 5), contributing to bloom formation and demonstrating adaptability to low light conditions (Biswas, 2022; Hoppe et al., 2018). Among the zooplankton, two taxa were identified as potential indicators: the planktonic heterotrophic radiolarian species Heteracon sp. and the mesozooplankton filter feeder appendicularian Oikopleura vanhoeffeni (Deibel, 1986, 1988) (Fig. 5). Likewise, the cold-water bivalve species Tridonta sp. which is commonly found in the North Atlantic and Arctic region (Marincovich Jr et al., 2002; Petersen, 2001), demonstrated strong potential as an indicator species. Within parasitic dinoflagellates, three potential indicators belonging to dino-group-I, mostly associated with sea-ice conditions (Clarke et al., 2019), were identified (Fig. S12). Choanoflagellate recorded in the study can be identified as a sea-ice associated group due to the presence of potential ArW indicator taxa (Buck and Garrison, 1988; Thomsen and Østergaard, 2017) (Fig. 5). These taxa were all consistently present throughout the study period, suggesting a stable association with cold marine conditions.

C. gelidus, O. vanhoeffeni, C. natans, Dino-group-I-clade-I spp., and Tridonta sp. exhibited the strongest ArW indicator potential, reflected by their consistent detection across samples and their concordant associations with multiple independent paleoenvironmental proxies. In contrast, MAST-9D sp., MAST-9A sp., and Cryothecomonas spp. showed the strongest indicators of AW influence (Fig. 5). However, their Spearman correlation coefficients between environmental variables ranged from 0.3 to 0.6 (p < 0.05), indicating a weak to moderate association. This may be due to a combination of interspecific competition and the influence of multiple external environmental variables in the study area. These interacting factors contribute to the complexity of the ecosystem and limit the effectiveness of using single-proxy approaches when interpreting the responses of indicator species in paleoenvironmental reconstructions. Further studies of eukaryotic communities in other Arctic regions are therefore needed to validate these taxa as robust indicator species.

5.3 Interactions within eukaryotic community structure in Storfjordrenna

The biodiversity of eukaryotic communities in Storfjordrenna was previously influenced by the interplay between ArW and AW masses, as well as sea-ice coverage over the past 13.30 kyr BP. Throughout the study period, eukaryotic alpha diversity remained relatively stable, with a notable exception during the transitions between major climatic intervals suggesting that these changes were driven primarily by species replacement, rather than by loss of richness or evenness (Figs. 4 and S2). This may imply that key trophic interactions and ecosystem functions persisted, reflecting functional resilience of Storfjordrenna eukaryotic communities during periods of environmental change. Notably, biodiversity peaks coincided with the presence of sea-ice margins and frontal zones, environments known to promote phytoplankton growth and primary productivity (Fig. 4).

Analysis of phytoplankton diversity revealed a consistent presence of green algae throughout the study period, except during the Bølling–Allerød interstadial, when diatoms dominated. Taxonomic abundance showed dynamic fluctuations, with gradual declines observed during the transitions between major climatic intervals. This suggests that environmental instability may have influenced the structure of the phytoplankton community.

Key contributors to primary productivity in the Storfjordrenna were diatoms (Thalassiosira spp., and Chaetoceros spp.), green algae (M. polaris), and autotrophic dinoflagellates (P. glacialis, and Gymnodinium spp.). Additionally, Spearman rank correlation analysis showed that the family Actinomonadaceae, mainly represented by Pseudopedinella sp., was positively associated with diatoms. Spearman rank correlation analysis also revealed a positive association between parasitic cercozoans and phytoplankton communities (Fig. 6, Table S6), indicating that the presence of parasitic cercozoans may play a significant role in shaping ecological interactions within phytoplankton assemblages (Bass et al., 2009; Cavalier-Smith and Chao, 2003; Hartikainen et al., 2014). Conversely, the parasitic apicomplexan family Lecudinidae was associated with zoobenthos (e.g., Heteroconchia and Ophiurida) and mesozooplankton (e.g., Malacostraca), may highlight their parasitic relationships with marine invertebrates (Rueckert et al., 2015) (Fig. 6, Table S6). Parasitic dinoflagellates (dino-group-II) were positively associated with haptophytes and diatoms. The parasitic nanoflagellate Pirsoniaceae: Pirsonia sp., demonstrated a positive correlation with autotrophic microbes, including, dinoflagellates, diatoms, and silicoflagellates (Kühn et al., 2004; Schweikert and Schnepf, 1997). This raises questions about the nature of their ecological interactions, and whether they are strictly parasitic or co-occur under similar environmental conditions. These findings emphasize the need for further investigation to understand the mechanisms driving these interactions and their broader implications for microbial community dynamics.

Cercozoans have emerged as the most dominant group within the eukaryotic community, in terms of both abundance and species richness. Cercozoan community included taxa previously recorded in various habitats including fresh and marine environments (Chantangsi and Leander, 2010; Irwin et al., 2019). Their occurrence across a wide range of environmental conditions highlights their ecological flexibility and broad tolerance. Although cercozoans as a group exhibit high richness, only a few lineages, such as the Imbricata-novel clade 2, Protaspa spp., Cryothecomonas spp., and Ascetosporea, persisted consistently throughout the study period, while most of the others were restricted to specific time intervals (Fig. S8). This suggests that, although cercozoans as a whole group may not be sensitive to environmental shifts, individual lineages are likely to be more responsive.

MAST species constituted a major microbial group within the eukaryotic community, with MAST-9 dominant overall and MAST-12 particularly prevalent during the late Holocene. These two MAST subgroups are commonly associated with temperate regions or extreme environments such as cold seeps (Lin et al., 2022; Obiol et al., 2024). The statistical analysis identified potential warm-water indicator species within these groups (Fig. 5, Table S5). Based on the co-occurrence relation, the MAST-12, and MAST-9 subgroup revealed a positive correlation (> 0.6) with the parasitic family of Pirsoniaceae and the phytoplankton families (Fig. 6, Table S6), providing insight into their ecological activities within the eukaryotic community. Further analysis of the ecological traits and distribution of micro eukaryotic taxa, mainly cercozoans and MAST, could provide deeper insights into their ecological responses within the Storfjordrenna ecosystem.

6 Conclusions

Using sedaDNA metabarcoding, we reconstructed the paleoecology of eukaryotic communities in Storfjordrenna over the last 13.30 kyr BP, elucidating their sensitivity and adaptability to environmental variables. Most eukaryotic ASVs were assigned ecological functional roles via higher taxonomic classification, considered reliable due to extensive ribosomal reference databases. However, taxa with complex life cycles (e.g., dinoflagellates, ciliates) required genus/species-level identification for more accurate functional assignment, necessitating a cautious approach. While some classification errors are inevitable, their proportion in the dataset is expected to be minimal. Overall, the eukaryotic biodiversity in Storfjordrenna remained relatively stable, except during transitions between major climatic intervals, indicating that community changes were primarily driven by species replacement rather than loss of richness or evenness, which maintained ecosystem function. Peaks of biodiversity coincided with the presence of sea-ice margins and frontal zones; environments known to foster favorable conditions for phytoplankton development. Cercozoans and MAST emerged as dominant groups, demonstrating their ecological flexibility and broad tolerance. This study revealed that primary productivity in the Storfjordrenna region was mainly driven by phytoplankton, including diatoms (Thalassiosira spp., Chaetoceros spp.), green algae (Micromonas spp.), and autotrophs dinoflagellates (P. glacialis) as well as mixoplankton species such as Pseudopedinella elastica. Our approach revealed that, despite significant species turnover, functional diversity and ecosystem functions remained largely stable, highlighting the resilience of Arctic planktonic communities. Several potential ASV-indicators were identified through multi-method analyses, including taxa associated with specific water masses. Our findings also underscore the complex interplay of environmental drivers shaping community composition, revealing both positive and negative associations among key microbial taxa. Our findings highlight the potential of sedaDNA for reconstructing past eukaryotic communities and detecting environmental change. However, to fully unlock the potential of sedaDNA approach in palaeoecological studies, it is essential to develop reference databases for accurate and precise taxonomic identification of sequences, and to provide modern reference data on the ecology and distribution of Arctic microbial eukaryotes. Therefore, improving taxonomic resolution and validating indicator taxa remain essential for establishing robust palaeoecological indicators in Storfjordrenna and the broader Svalbard region.