Marine ecosystem shifts with deglacial sea-ice loss inferred from ancient DNA shotgun sequencing

The study addresses how polar marine ecosystems respond over long timescales to sea-ice loss associated with climate warming. Despite extensive modern monitoring, long-term ecosystem responses, especially among zooplankton, fish, non-fossilizing algae, and benthic organisms, remain poorly understood. The Bering Sea has experienced declining sea-ice duration and extreme low-ice events, with documented short-term cascading effects on food webs. The authors aim to reconstruct organismal composition shifts from the Last Glacial Maximum through the Holocene using sedimentary ancient DNA (sedaDNA) to provide insights into long-term ecosystem dynamics, inform risk assessments, and identify potential trajectories under continued sea-ice decline.

Marine sediments preserve climate and sea-ice histories via proxies such as biomarkers (e.g., IP25 produced by sea-ice diatoms) and microfossils (e.g., diatom frustules). Alkenones (from haptophytes) inform late summer/early fall sea-surface temperatures. While these proxies have revealed sea-ice variability, many non-fossilizing groups (protists, zooplankton) and poorly represented taxa (fish) are not well captured. SedaDNA metabarcoding has shown past plankton responses to water mass characteristics and sea ice, but PCR-based approaches can suffer from biases related to fragment size and primer mismatches. Metagenomic shotgun sequencing expands taxa coverage, enables damage authentication, and reduces PCR bias. Co-occurrence networks help interpret multidimensional ecological data, but metagenomic sparsity can induce spurious correlations. Gaussian copula graphical models (GCGMs) can separate environmental effects from intrinsic taxa associations, handle compositionality, and reduce spurious links through mediator taxa. Such network approaches had not previously been applied to sedaDNA.

• Study site and paleoenvironmental context: Sediment core SO201-2-12KL was retrieved from the western continental slope of Kamchatka (53.993°N, 162.37°E; 2173 m water depth), ~70 km offshore. Prior multi-proxy reconstructions (diatom microfossils, IP25) indicate seasonal sea ice during the LGM, Heinrich Stadial 1, and Younger Dryas, and predominantly ice-free conditions during the Bølling-Allerød and Holocene.
• Sampling and chronology: In 2018, 25 one-centimeter samples were collected from a 9.05 m piston core archived at 4 °C since 2009. The established age model indicates sedimentation rates of ~80 cm kyr−1 during deglaciation and ~30 cm kyr−1 in the Holocene, giving mean temporal resolution of ~780 years; each 1-cm sample integrates ~15.5 years (late glacial) to ~28.7 years (Holocene) of DNA deposition. Five negative controls accompanied the samples.
• DNA extraction: Extractions used the DNeasy PowerMax Soil Kit from ~7.5 g sediment in a dedicated clean facility. Extracts were quantified (Qubit), concentrated as needed, diluted to 3 ng µL−1, and stored at −20 °C.
• Library preparation and sequencing: Single-stranded DNA libraries (optimized for degraded DNA) were prepared using 15 ng DNA. Double-indexing PCR (10–14 cycles) with indexed P5/P7 primers and AccuPrime Pfx was performed; library size distributions were checked (Agilent TapeStation). Libraries were pooled (negatives at 10:1 lower proportion) to 10 nM and sequenced as paired-end 2×125 bp: first pool on Illumina HiSeq 2500 (40.6 GB), second and third pools on Illumina NovaSeq SP flow cells (198.1 GB). Library blanks monitored contamination.
• Bioinformatics: Quality control by FastQC; duplicate removal by FastUniq and clumpify (BBMap); trimming and merging by fastp (quality, adapter, low-complexity filters; overlapping read merging with error correction). Taxonomic classification used Kraken2 against the NCBI nt database (downloaded June 2021), with confidence 0.2, run in single mode for merged and paired mode for unmerged reads. Family-level assignments were retained and merged across read types. To focus on ecologically relevant targets and avoid storage-related biases, retained groups included phototrophic bacteria, photo- and heterotrophic protists, marine macrophytes, and Metazoa of likely regional aquatic origin (fish filtered by FishBase occurrence in the subarctic North Pacific/Bering Sea). Potential contaminants and parasitic lineages (e.g., human pathogens, apicomplexans) were excluded due to database annotation issues and contamination signals in blanks. Families occurring in at least three samples with ≥10 counts were retained. Datasets were resampled 500 times to equalize sequencing depth (pelagic: 6593 counts/sample; benthic: 1839 counts/sample).
• Negative controls and authentication: After processing, blanks retained 0.37–2.14% of original reads; 1519 reads in total were assigned at family or lower levels, mostly Hominidae and common bacteria. Damage patterns indicating ancient DNA were assessed with HOPS v0.34 on key taxonomic groups, examining C→T substitutions in terminal positions.
• Environmental data and statistics: IP25 was available for the same core; an independent NE Pacific SST stack was used and both were interpolated to sample ages (with cautious handling of missing IP25 values for the last three samples). Pairwise Spearman rank correlations (Benjamini–Hochberg adjusted p-values) were computed, and positive correlations (ρ>0.4, adjusted p<0.1) formed undirected co-occurrence networks (igraph, Fruchterman–Reingold layout). To account for environmental covariation and mediator taxa, Gaussian copula graphical models were constructed using ecoCopula: stacked SDM with SST and IP25 as covariates, followed by graphical lasso (lambda=0.51). Network similarity was quantified via edge Jaccard index against randomized null-models. Temporal autocorrelation was tested (acf); stratigraphic plots and clustering (CONISS) were produced using rioja and vegan.

• Taxonomic recovery: Across 25 samples, 918,186,452 read pairs were generated; ~70.76% passed quality checks. 13,119,146 read pairs were taxonomically classified; ~98.5% remained unclassified, consistent with incomplete marine reference databases. Mean fragment lengths were 83–105 bp, consistent with sedaDNA.
• Pelagic composition: 167 families (456,058 read pairs after filtering/resampling) were identified among pelagic/sea-ice/demersal groups. Pooled across samples, read shares after resampling: fish 66.8%, phototrophic bacteria 12.8%, phototrophic protists 10.2%. Among phototrophs, Cyanobacteria 11.5%, Chlorophyta 3.6%, Chlorobi 0.8%, diatoms 4.7%, haptophytes 1.0%. Zooplankton assignments were lower (heterotrophic protists 2.1%; copepods 0.3%). Among metazoans, dominant fish families included Salmonidae 12.1%, Serranidae 8.0%, Gadidae 6.0%. Marine mammals detected included Delphinidae 0.3%, Balaenopteridae 0.3%, Otariidae 0.1%, Phocidae 0.4%, and Mustelidae (sea otter) 2.9%.
• Networks and environmental associations: The Spearman network (167 families with ≥10 counts in ≥3 samples) comprised two modules with 148 nodes and 446 positive edges (ρ>0.4, adjusted p<0.1): a sea-ice module (eight nodes positively correlated with IP25) and a sea-ice-free module (13 nodes showing positive trends with higher SST). The ecoCopula (GCGM) network had 167 nodes and 474 edges; modules reflected functional groups (e.g., diatom- and fish-dominated modules). Despite controlling for IP25 and SST, families positively correlated with IP25 in the Spearman network remained co-associated, implying underlying biogeographic or unmodeled environmental structuring. There was a 34% overlap of positive associations between networks, exceeding null expectations.
• Sea-ice module characteristics: Elevated during LGM/HS1/YD, including diatom families (Bacillariaceae 0.8%, Stephanodiscaceae 0.6%, Thalassiosiraceae 1.1%, Triparmaceae 0.3%, Phaeocystaceae 0.6%), cold-adapted chlorophytes (Mamiellaceae 0.5%, Bathycoccaceae 1.7%), copepods (Calanidae 0.03%, Metridinidae 0.03%), and Gadidae. The eight IP25-correlated families had 109 links to 28 families, with 42.2% of links involving diatoms (15.6% centric, 26.6% pennate), underscoring diatoms as central primary producers under seasonal sea ice.
• Ice-free module characteristics: Warmer Holocene phases were enriched in families typical of warmer waters (e.g., Hemiaulaceae, Rhizosoleniaceae, Chloropicaceae). The 28 SST-correlated families (adjusted p<0.1) connected via 114 links to 51 families across 13 functional groups; these nodes had fewer neighbors than IP25-correlated families (p<0.001), with neighbors mostly fish (43% of links), phototrophic bacteria (25.4%), Chlorophyta (11.4%), and centric diatoms (1.8%). Salmonidae and Clupeidae (Pacific herring) had higher read counts in warm phases (e.g., early Holocene), consistent with modern observations of increased distributions in warm years and known regional spawning grounds.
• Benthic composition and pelagic–benthic coupling: Benthic assignments totaled 263,193 reads (45,975 after resampling), lower than pelagic (476,058 originally; 210,800 after resampling). Detected benthic families typical of shelf/upper slope included Pectinidae 22.7%, Mytilidae 1.5%, Asteriidae 14.1%, Priapulidae 6.9%, Nephtheidae 1.0%, Pocilloporidae 2.8%, Strongylocentrotidae 4.5%; nearshore macrophytes were present (Zosteraceae 0.8%; brown algae 2.9%; red macroalgae 10.9%). The pelagic:benthic read ratio was highest in late glacial/YD and declined with warmer SSTs (Spearman ρ = −0.39, p=0.052). Zosteraceae increased in the Holocene and correlated with summer/early fall SSTs (Spearman ρ = 0.73, adjusted p=0.026). Kelp (Laminariaceae) were more abundant in the late glacial and co-occurred with Pacific cod, suggesting functional benthic–pelagic links.
• Ecosystem shift over deglaciation: A transition from a sea-ice-adapted system (diatoms–copepods–Gadidae) during late glacial to an ice-free Holocene system dominated by cyanobacteria and increased Salmonidae and Clupeidae. Results imply a shift from nano-/micro-phytoplankton to pico-sized phototrophs, with consequences for carbon export and benthic food supply.
• Methodological advance: Shotgun sedaDNA extended detection of fish and marine mammals back ~20,000 years, beyond previous records, demonstrating feasibility for millennial-scale biotic reconstructions.

The findings address the central question of how deglacial sea-ice loss restructured marine ecosystems in the western Bering Sea. Coherent associations between sea-ice proxies and taxa indicate that seasonal sea ice selected for diatom- and copepod-based food webs supporting codfish, whereas ice-free, warmer conditions favored cyanobacteria and green algae, along with higher abundances of salmon and herring. Network analyses, including GCGMs that control for environmental covariates and mediator taxa, substantiate robust co-association patterns across functional groups, suggesting shared environmental responses rather than short-term biotic interactions. The observed shift toward pico-sized phototrophs under warmer, stratified, lower-nutrient surface waters implies reduced efficiency of the biological carbon pump relative to diatom-dominated systems and potentially diminished benthic food quality/quantity. Benthic signals (e.g., increased seagrass with SST and reduced pelagic:benthic ratio under warmth) reinforce tighter pelagic–benthic coupling changes, with macrophyte expansions providing habitat/spawning benefits to salmon and herring. Conversely, cod and pollock, reliant on large lipid-rich copepods, may be disadvantaged under prolonged warm, low-ice regimes. Overall, the study demonstrates that sedaDNA shotgun sequencing can reconstruct long-term ecosystem dynamics across diverse taxa, informing ocean and cryosphere risk assessments by identifying likely future shifts under continued sea-ice decline.

This study reconstructs a deglacial ecosystem transition from a sea-ice-adapted, diatom–copepod–codfish system to an ice-free, pico-phytoplankton (cyanobacteria)–dominated system, with increased salmon and herring. By integrating sedaDNA shotgun metagenomics with proxy records and robust network analyses, the work extends the temporal detection of fish and marine mammals to ~20 ka and reveals consistent taxonomic shifts tied to sea-ice loss and warming. The results imply changes in carbon export pathways and benthic food supply, and forecast northward expansions of seagrass, salmon, and herring under ongoing warming, while cod/pollock may face disadvantages. Future research should improve reference databases for marine taxa, expand spatial and temporal sedaDNA records to refine biogeographic inferences, integrate quantitative biomarker and nutrient data to constrain productivity and export, and develop approaches to attribute organic matter sources and functional traits to past community changes.

• Taxonomic classification limitations: Approximately 98.5% of reads remained unclassified due to incomplete reference databases, especially for protists and deep-sea benthos, reducing taxonomic resolution (family-level focus) and potentially biasing diversity estimates.
• Temporal averaging: SedaDNA samples integrate decades, limiting inference of direct ecological interactions; networks reflect shared environmental responses rather than real-time co-occurrence.
• Proxy biases and environmental data: IP25 and SST stacks have inherent uncertainties and potential regional biases; some IP25 values required imputation for the youngest samples. Limited sample size advises caution in correlation interpretations.
• Potential storage effects: Long-term cold storage (since 2009) could bias prokaryotic and fungal composition; analyses mitigated this by focusing on phototrophs, protists, macrophytes, and metazoans.
• Low benthic assignment rates: Particularly for deep-sea taxa, many species are undescribed or missing from databases, potentially underrepresenting benthic diversity.
• Contamination control constraints: Although blanks showed low read counts and damage patterns support ancient origin, some background contamination (e.g., Hominidae, common bacteria) was detected and removed during filtering.