Multi-millennial genetic resilience of Baltic diatom populations disturbed in the past centuries

20 Little is known about the genetic diversity and stability of natural populations over millennial time 21 scales, although the current biodiversity crisis calls for heightened understanding. Marine 22 phytoplankton, the primary producers forming the basis of food webs in the oceans, play a pivotal role 23 in maintaining marine ecosystems health and serve as indicators of environmental change. This study 24 examines the genetic diversity and shifts in allelic composition in the diatom species Skeletonema 25 marinoi over ~ 8000 years in the Baltic Sea by analyzing chloroplast and mitochondrial genomes. 26 Ancient environmental DNA (aeDNA) from sediment cores demonstrates stability and resilience of 27 genetic composition and diversity of this species across millennia in the context of major climate 28 events. Accelerated change in allelic composition is observed from historical periods onwards, 29 coinciding with times of intensifying human activity, like the Roman Empire, the Viking Age, and the 30 Hanseatic Age, suggesting that anthropogenic stressors have profoundly impacted this species for the 31 last two millennia. The data indicate a very high natural stability and resilience of the genomic 32 composition of the species and underscore the importance of uncovering genomic disruptions caused 33 by human impact on organisms, even those not directly exploited, to better predict and manage future 34 biodiversity. 35 36

37 Human activities have had a profound impact on natural environments, resulting in the endangerment 38 of a multitude of ecosystems, including marine environments 1. Phytoplankton forms the basis of 39 marine food webs, and, as crucial components of marine ecosystems, changes in this group are central 40 to understanding ecosystem shifts2. Among the phytoplankton, diatoms, which play a significant role 41 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 2 in global biogeochemical cycles, are sensitive bioindicators 3, as the different species have specific 42 requirements. Skeletonema marinoi , a prominent marine diatom species 4, is influenced by 43 temperature, migration, and human activities5. For this species, changes in bloom timing6 and optimal 44 growth temperature7 have already been observed as a result of recent climate change. Investigation 45 of marine phytoplankton response to current global changes are mostly based on taxonomic diversity, 46 but recent studies show that diatoms also quickly adapt to environmental changes through acclimation 47 and genetic adaptation on a population genomic level 8. In fact, intraspecific genomic variation is an 48 important factor in the response and resilience of diatoms to environmental perturbations 9. Thus, 49 investigating the genomic composition of diatom populations and its changes is crucial for 50 understanding their adaptability to environmental changes. 51 The Baltic Sea, despite its relatively brief history, has undergone significant transformations, making it 52 suitable for studies addressing organism responses to environmental changes. For the last ~ 10,000 53 years, these transformations include a transition between fresh and brackish water around 8,500 -54 8,000 years ago 10 and major climate changes since then 11. Due to its shallow enclosed nature and a 55 steep salinity gradient, modern biodiversity is relatively low 12. Current pressures such as pollution, 56 eutrophication and global warming present significant challenges to the biota13. 57 Sediments can serve as ecological archives when undisturbed, offering insights into past 58 environmental changes and biodiversity shifts 14, and provide long time series of phytoplankton 59 dynamics. Phytoplankton archives in sediments include biomarkers, microfossils, (living) resting stages, 60 and sedimentary ancient DNA (sedaDNA). These remains contain information on biodiversity and can 61 reveal changes in adaptive traits 7,15. Here, we develop a time series of population -level responses of 62 the diatom S. marinoi in the Baltic Sea, by employing a targeted approach to enrich chloroplast and 63 mitochondrial genomes using hybridization enrichment. This allows us to analyze genetic diversity and 64 shifts in allelic composition across complete organellar genomes. We investigate 1) the degree of 65 stability, response and resilience of S. marinoi to natural environmental changes in the Baltic Sea over 66 the last ca. 8000 years, and 2) identify changes in the genetic composition of populations in recent 67 centuries of heightened anthropogenic activities. 68 69

70 We analyzed sediment samples of two cores located in the Baltic Sea, the Eastern Gotland Basin (EGB) 71 and the Gulf of Finland (GOF) (Fig. 1A), covering a time span up to approximately the Ancylus Lake 72 (extrapolated age: 8148 cal yr BP) (BP; with present = 1950 Common Era) (Fig. 1B). We investigated 73 the population structure of S. marinoi over this period at the two locations by focusing on Single 74 Nucleotide Polymorphisms (SNPs) on the organelle genomes (Fig. 1C). This was achieved by target 75 enrichment of the chloroplast and the mitochondrion. 76 This study provides an analysis of the genetic dynamics of S. marinoi populations in the EGB and GOF 77 over millennia, revealing insights into the effects of both natural environmental changes and 78 anthropogenic activities. The use of RNA baits for hybridization enrichment resulted in a higher degree 79 of resolution, enabling the identification of a greater number of SNPs in both mitochondrial and 80 chloroplast genomes. The analysis revealed that climate events exert an influence on genetic 81 variations, with specific climate events aligning with temporal variations in the population’s genetic 82 makeup, but that the genetic composition reverted back to a state that was mostly stable across 83 millennia in between these events. Longer lasting patterns of genetic change over time were observed 84 at both sites in the past two millennia, coinciding with periods of intensified anthropogenic activities. 85 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 3 These findings imply high resilience of S. marinoi populations across time and climate events across 86 millennia, along with more recent changes. 87 88 Fig. 1. A) Location of the coring sites and corresponding cores in the central Baltic Sea. Gulf of Finland (GOF; core 89 EMB262/12-3GC) and Eastern Gotland Basin (EGB; core EMB262/6 -30GC). B) Sample distribution (in cm) in the sediment 90 cores and corresponding ages. C) Number of SNPs called per 100 base pairs across the genome for each sample. The 91 mitochondrion and chloroplast genomes of S. marinoi are shown. The gaps indicate the repeat masked regions. 92 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 4 General data assessment 93 DNA vs. RNA baits 94 In order to study past S. marinoi population dynamics in the EGB and GOF, we employed DNA and RNA 95 baits for hybridization enrichment of full chloroplast and mitochondrial genomes. The use of RNA baits 96 resulted in a higher degree of resolution. Following trimming, DNA and RNA baits yielded 34.9 million 97 and 51.8 million sequences, respectively. Of these, 34.8 million (DNA) and 35.3 million (RNA) reads 98 were successfully mapped to reference organelles (Fig. 2A and Supplementary Table S1). In the case 99 of RNA baits, 6.41 million reads were mapped to the mitochondrion, while 28.9 million were mapped 100 to the chloroplast. For DNA baits, the numbers were 6.27 million and 28.6 million, respectively. The 101 RNA baits identified a greater number of SNPs in both the mitochondrial (301 vs. 92) and chloroplast 102 (1716 vs. 403) genomes. Consequently, the analysis focuses on the results obtained from RNA baits. 103 The results of the DNA bait analysis can be found in the supplementary section 2.1. 104 The average coverage achieved for the experimental samples was generally high, with an average 105 mitochondrial coverage of 78.77 and an average chloroplast coverage of 575.5. This coverage 106 demonstrates the efficiency of the hybridisation capture approach and ensures reliability of the 107 genomic data collected. Furthermore, the Extraction Blanks and Library Blanks produced low to non -108 meaningful coverage (see Supplementary Table S2). All information about the samples, corresponding 109 age, location and other metadata is available in Supplementary Table S3. 110 Ancient DNA – damage pattern 111 We analyzed damage patterns in ancient DNA (aDNA) samples, focusing on C-to-T substitutions. These 112 substitutions are key aDNA authentication markers due to their pervasiveness in aDNA and their 113 function in indicating cytosine deamination, a prevalent form of DNA damage. No significant 114 correlation (R = 0.18, p = 0.315) was found between mapping coverage and substitution frequency, 115 suggesting that damage pattern does not affect the mapping coverage. The degree of damage was 116 found to significantly increase with the age of the samples (R = 0.65, p = 0.0003, Fig. 2B). A greater 117 number of sequences were mapped to the chloroplast due to its larger genome size. This allows a 118 greater number of reads to be used to detect damage patterns, increasing the overall value of the 119 data. 120 Genetic diversity and spatial differentiation 121 A comparison of the genetic differences between S. marinoi retrieved from the EGB and the GOF 122 respectively, revealed that while the two locations exhibited some differences, they also shared a 123 considerable degree of similarity. This was evidenced by the FST value, a measure used in population 124 genetics to quantify the genetic differentiation between populations, which ranged between 0.05 and 125 0.1 (see Supplementary Material, Fig. S6). Notably, the genomic data from EGB shows a higher level of 126 genetic diversity than the one from GOF. This is presented by different measures of distinct variant 127 diversity, with significant differences observed (p-values: 0.002, 0.005, 0.006, see Fig. 2 C-E). 128 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 5 129 Fig. 2: Data assessment: A) Sequencing Yield and Single Nucleotide Polymorphisms (SNPs): Post-trimming sequencing yield 130 for DNA and RNA baits, and the number of SNPs detected in mitochondrial and chloroplast genomes for each bait set. B) 131 C-to-T Substitution: Analysis of C -to-T substitutions over time for each organelle from all data (GOF & EGB). C -E) 132 Comparative analysis of normalized nucleotide diversity and diversity indices. C) Relative Nucleotide Diversity (PI), D) 133 Shannon, and E) Simpson between Eastern Gotland Basin (EGB) and Gulf of Finland (GOF). Each bar graph shows data for 134 both sites with significant p-values. 135 136 Effects of environmental change 137 Our study on S. marinoi organelles at two locations in the Baltic Sea reveals that across several 138 millennia the genetic composition of the population remained stable, punctuated by differences in 139 specific periods, putatively coinciding with environmental changes. A Principal Component Analysis 140 (PCA), conducted to investigate differences in allelic composition between sediment layers, showed 141 that the majority of samples are situated within a primary cluster, but a number of smaller clusters are 142 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 6 found that correspond with periods marked by specific environmental events (Fig. 3A). Together, PC1 143 and PC2 describe 42.6% of the data variance. 144 The retrieved specific clusters (Fig. 3A) correspond to specific climate periods, including the Holocene 145 Thermal Maximum (HTM, 10,000-6,000 cal yr BP), Little Ice Age/Late Antique Little Ice Age (LIA, 550 -146 100 cal yr BP/LALIA, 1,414 -1,290 cal yr BP), Medieval Climate Anomaly (MCA), and Modern Warm 147 Period (MWP). We verified that existing variation in allelic composition is not a result of the C -to-T-148 substitution rate by performing a PERMANOVA, which showed that the C -to-T-substitution rate 149 accounted for approximately 6.7% of the variation in the PCA space (R² = 0.067, p = 0.19), though this 150 was not statistically significant (see Supplement, Fig. S5). Additionally, we performed a PERMANOVA 151 to assess the impact of Total Organic Carbon (TOC) on the variation in the PCA space, which showed 152 that TOC accounted for approximately 4.5% of the variation (R² = 0.045, p = 0.328), also not statistically 153 significant. 154 In the EGB, the period under consideration extends to the Ancylus Lake (extrapolated age: 8,148 cal yr 155 BP), and responses are observed during the HTM, LALIA, MCA, and LIA periods. The GOF record, which 156 covers the last ca. 4,300 years, displays changes during the LIA and MWP periods. During periods 157 without profound climate events, we observe a consistency in allelic composition, which points 158 toresilience of the algal populations over time. However, within warm periods (e.g. HTM and MCA), 159 there is a notable degree of variability, particularly for EGB. This pattern of long -term stability and 160 resilience, interrupted by distinct variation, is visualized in Fig. 3B, depicting changes along PC1 in both 161 cores. This shows temporal variations that align with specific climate events, suggesting that the 162 population’s genetic makeup is undergoing changes in response to these events. Additionally, 163 nucleotide diversity increased significantly during the transition from the Littorina Sea to the Modern 164 Baltic Sea, around 4,000 cal yr BP, and then successively decreased again (Supplementary Material, 165 Fig S7). 166 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 7 167 Fig. 3: Principal component analyses of the allelic composition of Skeletonema marinoi organelles, demonstrating the 168 genetic diversity and population structure. A) PC1 and PC2 are categorized by climate events, total organic carbon (TOC) 169 and the age of each sample is given for both sites. B) PC1 versus time for both sites. The red dashed line marks the transition 170 from the Littorina Sea to the Modern Baltic Sea. Climate events include the Holocene Thermal Maximum (HTM), Late 171 Antique Little Ice Age (LALIA), Medieval Climate Anomaly (MCA), Little Ice Age (LIA), Modern Warm Period (MWP) and NA 172 = No event 173 174 Allele turnover, a measure of the rate at which new alleles replace old ones, provides insights into 175 genetic variation over time. Figure 4 displays patterns of genetic change over time. A Generalized 176 Additive Model (GAM) was used to analyze allele turnover as a function of time. 177 At EGB, an increase in allele turnover is observed around 1,500 -1,000 cal yr BP, indicating a potential 178 for heightened genetic change from the Medieval period onwards. At GOF, turnover remains 179 consistent until around 94 to -34 cal yr BP. While industrialization and anthropogenic impacts may be 180 associated with recent genetic shifts, establishing a causal link between these changes and earlier 181 historical periods requires caution. 182 183 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 8 184 Fig. 4: Temporal Analysis of Allele Turnover in Skeletonema marinoi. This figure presents the Generalized Additive Model 185 (GAM) fitted to the allele turnover over time (BP). The data includes both Eastern Gotland Basin (EGB) and Gulf of Finland 186 (GOF) sites. Key environmental events, such as the Holocene Thermal Maximum (HTM), Late Antique Little Ice Age (LALIA), 187 Medieval Climate Anomaly, Little Ice Age (LIA), and Modern Warm Period (MWP), are noted. Additionally, anthropogenic 188 periods and associated shipping activity estimated from literature sources are shown. 189 190

191 The millennial-long stability, spanning from the Ancylus Lake (approx. 8,148 calibrated years BP), to 192 Littorina Sea and Modern Baltic Sea to 1,500 –1,000 cal yr BP in the Eastern Gotland Basin (EGB), and 193 from around 4300 to 94 –34 cal yr BP in the Gulf of Finland (GOF), aligns with findings from other 194 phytoplankton studies, albeit with a different temporal perspective. The findings of these studies 195 indicate that long -term genetic stability is frequently maintained in marine populations due to their 196 substantial population size and high dispersal capacity, even in the face of significant environmental 197 changes 16. These factors contribute to a buffering effect against genetic drift and local extinctions, 198 allowing populations to maintain genetic diversity over long periods of time. The ability of S. marinoi 199 to cope with changing environmental conditions through reversible shifts in genetic composition as 200 shown here for the Baltic Sea further supports this stability. This resilience is a common trait among 201 phytoplankton, allowing them to respond to environmental changes without losing overall genetic 202 diversity17. The observed stability is expected given the ecological and evolutionary characteristics of 203 phytoplankton populations. Large population sizes, high reproductive rates, and the potential for gene 204 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 9 flow between populations contribute to the maintenance of genetic diversity and stability over long 205 periods of time 18. Other findings, such as studies of diatoms and other phytoplankton species have 206 reported patterns of genetic stability and adaptability of phytoplankton to environmental change 19. 207 These findings imply that the observed stability in S. marinoi populations is not only expected, but also 208 indicative of the broader ecological and evolutionary dynamics of marine phytoplankton. Our data 209 demonstrates that this stability can be upheld across several millennia. 210 However, this stability is not without punctuated interruptions. Our results show that climatic events 211 have an influence on the allelic composition of S. marinoi populations. For example, shifts in the 212 genetic composition of the EGB population were observed during the Holocene Thermal Maximum, 213 the Late Antique Little Ice Age, and the Medieval Climate Anomaly. Similarly, the GOF population 214 showed changes during the Little Ice Age and the Modern Warm Period. These shifts highlight the 215 ability of S. marinoi to cope with different climatic events. The discrepancy between EGB and GOF may 216 be due to the slightly higher genetic diversity in the EGB (Fig. 2C -E), which may facilitate adaptation. 217 The differences between the two sites can also be attributed to their distinct geographical 218 characteristics. The EGB, located in the Baltic Proper, may experience an influx of migrants from other 219 populations, increasing genetic diversity and resilience. In contrast, the genetic composition of the GOF 220 population could be influenced by its unique environment, characterized by freshwater influx and 221 lower salinity levels20, as well as more extensive ice cover in the northern regions of the Baltic Sea21. In 222 particular, samples from the Little Ice Age and the Modern Warm Period show more pronounced 223 responses in the GOF. The allelic composition remains relatively constant in the absence of climatic 224 events, indicating a stable genetic structure under stable conditions and once a shift has manifested. 225 After this long period of punctuated genetic stability, recent centuries reveal a novel pattern of rapid 226 allelic turnover, and this coincides with increased human activity. This rapid turnover is particularly 227 evident in the GOF during periods such as the pre-Roman Iron Age, the Viking Age, the Hanseatic Age, 228 and the Industrial Revolution. These results highlight the possible impact of human activities on the 229 genetic dynamics of S. marinoi populations. This increased activity is manifested by increased shipping 230 activity in the Baltic Sea 22,23. Ballast water exchange 24 in the last two centuries could have been a 231 specific agent of increased changes in haplotype composition, causing rapid translocation of 232 populations. This rapid allelic turnover contrasts with the more gradual and reversible shifts observed 233 in response to natural climatic events, highlighting the influence of anthropogenic factors on the 234 genetic composition of these populations. Due to its coastal position, the GOF is more exposed to these 235 changes, and perhaps the phytoplankton is also more directly affected. 236 The rapid allelic turnover observed in recent centuries, particularly in response to human activities, 237 underscores the influence that anthropogenic factors may have on the genetic dynamics of 238 populations. Our study on the population dynamics of S. marinoi populations provides valuable insights 239 into the long-term resilience and adaptability of marine phytoplankton. The observed genetic stability 240 over several millennia, punctuated by reversible shifts in response to climatic events, underscores the 241 inherent resilience of these populations. However, the recent rapid allelic turnover highlights the 242 potential impact of human activities on population dynamics and structure, and thus presents an 243 opportunity for further investigation. 244 At the same time it is important to acknowledge the limitations of our study. The limited number of 245 samples included may affect the overall strength and generalizability of our findings. In addition, the 246 lack of precise tests limits our ability to draw definitive conclusions about long-term genetic trends and 247 the full extent of human impact. To date, we also lack a sufficient understanding of the distribution 248 and taphonomy of ancient eDNA, but in sediments, DNA of small organisms, such as phytoplankton, 249 has been shown to give a representative signal of a water body 25. The data of the two cores display a 250 high level of congruency, and despite the limitations, this indicates a link between environmental 251 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 10 change, human activity, and genetic dynamics. Future studies with larger sample sizes and e.g. more 252 locations are needed to further explore the complex interactions between environmental change and 253 genetic responses. 254 The results of our study provide a solid base for future research on the effect of current global change 255 on the genetic composition of natural populations across millennia. By understanding the genetic 256 stability and adaptability of S. marinoi populations, we can better understand ecological and 257 evolutionary dynamics of marine phytoplankton. This knowledge is crucial for predicting how these 258 populations may respond to current and future environmental changes, including climate change and 259 direct anthropogenic pressures, and thus offer a valuable foundation for future conservation efforts. 260 Our study provides novel insights into the genetic resilience and adaptability of S. marinoi populations. 261 It demonstrates that these populations maintain stability over millennia and respond dynamically to 262 both natural climate fluctuations and human-induced changes. In sum, it highlights the importance of 263 considering both natural and anthropogenic drivers and their relative impacts when assessing and 264 predicting genetic resilience of marine ecosystems. 265 266 267

& Methods 268 Study area 269 The Baltic Sea is a relatively young and shallow brackish water system currently facing a significant risk 270 of increasing hypoxia, a condition characterized by low oxygen levels and bottom water anoxia. This 271 phenomenon is primarily attributable to the reduction in dissolved oxygen and nutrient discharge 272 resulting from warmer water temperatures, which in turn promote the formation of algal blooms 26. 273 These blooms decompose and consume oxygen, leading to the development of hypoxia. In recent 274 history, the Baltic Sea has experienced a notable increase in eutrophication, a process driven by the 275 excessive input of nutrients such as nitrogen and phosphorus from agricultural runoff, wastewater 276 discharge, and industrial activities 13. The enrichment of nutrients has resulted in the proliferation of 277 phytoplankton and algal blooms, which, upon decomposition, have further depleted oxygen levels in 278 the water. The Baltic's unique position as a land-enclosed sea with high human and industrial activity, 279 coupled with its low biodiversity, makes it particularly vulnerable to these changes. 280 The Baltic Sea's history is well known. It has undergone substantial environmental changes in the past. 281 After the end of the last glaciation approximately 10,000 years ago, it went through several fresh- and 282 saltwater stages, starting as a large meltwater lake, then transitioning into various states of salinity 283 due to glacial retreat and land rebound27. These stages included the Baltic Ice Lake, the Yoldia Sea, the 284 Ancylus Lake, and the Littorina Sea 28. The Baltic Sea’s temperature and salinity have fluctuated over 285 time, with notable increases during the Holocene Thermal Maximum and the Medieval Climate 286 Anomaly, and a recent increase since 185027,29,30. 287 Within the Baltic Sea this study covers a time series of two distinct locations: The Eastern Gotland 288 Basin, located in the Baltic Proper with a profound water depth of max. 249 meters, and the 289 comparatively shallower Gulf of Finland, reaching a depth of 81 meters. At both locations, we collected 290 sediment cores to inform on the effects of environmental changes on the population genetic 291 composition of S. marinoi. 292 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 11 Sampling 293 Two sediment cores were taken (Fig. 1A) in April 2021, during expedition EMB262. 1) Eastern Gotland 294 Basin (EGB, 57°17.004'N, 020°07.244'E, 241 m water depth, core EMB262/6-30GC); 2) Gulf of Finland 295 (GOF, 59°34.443’N, 023°36.461’E, 81 m water depth, core EMB262/12-3GC). The cores (Fig. 1B, ca. 500 296 cm) were taken using a gravity corer (GC). The sampling for the EGB began at a depth of 32 cm. The 297 cores were sampled on board of the research vessel using sterile syringes according to Epp et al.31 and 298 immediately frozen for storage. 299 Core dating was performed as described in Schmidt et al. 32. This was done by correlating organic 300 carbon records from post -Littorina transgression sediments in the Baltic Sea. The chronology of the 301 cores is determined by the relative S content and the Br/K ratio. The Br/K ratio reflects changes in the 302 bulk organic carbon content of the sediments. The data are visually matched with XRF and organic 303 carbon data from dated sediment cores from the same or nearby locations for different historical 304 periods. For the sediment core from the Gulf of Finland, an independent Bayesian age model was 305 constructed based on radiocarbon dates of bulk organic matter in the sediments. This allows the 306 sediments to be accurately dated to around 2300 BC. The same core samples were used as in Schmidt 307 et al.32 and in this study. 308 309 Laboratory work 310 DNA extraction 311 The PowerSoil Pro Kit from Qiagen (Hilden, Germany) was used to extract DNA from both sediment 312 cores (n=26, see supplementary Table SX) covering various depths (ages). For the extraction 0.5 g of 313 sediment was used per sample, resulting in 26 samples and 2 extraction blanks. The extraction process 314 was carried out according to the manufacturer's protocol with some modifications including an 315 overnight incubation at 56 °C with the addition of 20 μL proteinase K (20 mg/ml) to enhance sample 316 lysis. The washing steps were performed using a Qiagen Vacuum Pump (Hilden, Germany). 317 Subsequently, samples were centrifuged for 3 minutes. After centrifugation, the elution process was 318 performed in two steps. Each part involved adding 75 μL of elution buffer to the sample, allowing it to 319 incubate for 5 minutes, and then collecting the eluate. Both eluates were collected in the same tube. 320 Library preparation 321 A single -stranded library approach according to Gansauge et al. 33 was performed to maximize the 322 retrieval of genetic information from these degraded sedaDNA samples. This method, optimized for 323 fragmented DNA, has been shown to significantly increase the sequence yield and preserves unique 324 molecules, thereby providing a more comprehensive genomic analysis. The library preparation was 325 conducted with 40 ng input DNA. All 26 samples and 2 extraction blanks were processed in the library 326 preparation. Additionally, 2 library blanks were created. 327 After the indexing PCR, each library was cleaned using the MinElute PCR purification Kit (50,250; 328 Qiagen) resulting in a final amount of 20 μL per sample (For more detailed info, see Supplement Section 329 1.1.). The samples were then analyzed on a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, 330 USA) using the Agilent Bioanalyzer High Sensitivity DNA Analysis Kit (Agilent Technologies, Santa Clara, 331 CA, USA). Due to adapters still being present in the samples, we performed an additional purification 332 using the HighPrep TM bead cleanup (MagBio Genomics Inc., Gaithersburg, MD, USA). 1.6x of 333 HighPrepTM PCR reagent was used. The DNA concentration was measured using ds -DNA HS Assay Kit 334 and the Qubit® 4 fluorometer (Invitrogen Thermo Fisher, Waltham, MA, USA). The libraries were 335 combined into pools of four, each with the same DNA amount. The extraction and library blanks were 336 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 12 incorporated in the same volume as the sample with the lowest concentration. The libraries were then 337 further used for hybridization enrichment of chloroplast and mitochondrial genomes. 338 Hybridization enrichment 339 In this study, a hybridization enrichment approach was employed. Hybridization enrichment allows for 340 targeted isolation of specific sequences, improving sequencing efficiency by reducing the presence of 341 non-target DNA. Thus, providing a better representation of the studied species, enabling a more 342 detailed and precise analysis of its genetic composition. 343 Bait design 344 The design of the bait set was carried out for two target reference organelle genomes: the 345 mitochondrion and the chloroplast (GenBank PRJNA493755). Together, these genomes consist of 346 170,788 nucleotides (nt), with 127,202 nt from the Chloroplast and 43,586 nt from the mitochondrion, 347 and an average GC content of 30.9%. To improve the accuracy of the baits, the contigs were 348 softmasked for simple and low-complexity repeats. 349 Subsequently, baits were designed with a length of 80 nt and a tiling coverage of 4x. This means that, 350 on average, each nucleotide in the target sequence is covered by four different baits, thereby ensuring 351 redundancy and enhancing the probability of successful hybridization. This process led to the 352 generation of 10,000 probes. The baits were then filtered based on softmasking for simple and low 353 complexity repeats. Only those with 35% or less softmasking retained. A total of 7,985 baits passed 354 this filtering step. 355 Enrichment 356 This bait design was then used for two different bait sets: 1) RNA baits (BioCat, Heidelberg, Germany); 357 2) DNA baits (Integrated DNA Technologies, IDT). The libraries were enriched with both RNA and DNA 358 baits, resulting in two sets of libraries: 1) RNA baits enriched; 2) DNA baits enriched. General steps are 359 explained in the following. More detailed information on the used protocols can be found in the 360 supplementary data (section 1.2.). 361 The enrichment with RNA baits was conducted according to the myBaits v.5.02 manual. A two -round 362 enrichment strategy was performed. The hybridization mix was prepared and transferred to a rotation 363 oven once the 24-hour incubation at a hybridization temperature of 63°C started. In the following bead 364 cleanup, the libraries were resuspended and amplified. Each 27.4 μL PCR reaction included the 365 following nine components: 1) 3.45 μL of DEPC treated H 2O; 2) 2,5 μL 10X HiFi PCR Buffer (Invitrogen 366 Thermo Fisher, Waltham, MA, USA); 3) 0.25 dNTPs (25mM); 4) 1 μL BSA (20 mg/ml); 5) 1 μL MgSO4 (50 367 mM); 6) 0.2 μL Platinum ™ Taq DNA -Polymerase High Fidelity (5 U/ μL) (Invitrogen Thermo Fisher, 368 Waltham, MA, USA); 7) 2 μL IS5_bridge_P5 (10 μM); 8) 2 μL IS5_bridge_P5 (10 µM), 9) 15 μL library 369 pool. The amplification was run using the following settings: 1 minute initiation at 94°C, 14 cycles of 370 15 seconds denaturation at 94°C, 20 seconds annealing at 60°C and 1 minute extension at 68°C, 371 followed by 2 minutes of final elongation at 68°C and afterward the sample was stored at a 372 temperature of 20°C. Afterwards, the samples were purified with beads and eluted in 10 μL. The 373 second round of enrichment was performed in a similar way. The PCR was conducted in 8 cycles and 374 the pools were eluted in 17 μL after the bead purification (see more detailed in supplement section 375 1.2.2.). 376 Enrichment with DNA baits was also performed in two rounds following the IDT (Integrated DNA 377 Technologies, Löwen, Belgium) xGen TM hybridization capture of DNA libraries protocol (option: 378 AMPure XP Bead DNA concentration protocol), with the following modification. The hybridization 379 temperature was set to 63°C. Once this temperature was reached on the thermal cycler, the tubes 380 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 13 were moved to a rotation oven and incubated for 16 hours (see more detailed in supplement section 381 1.2.1.). 382 All resulting pools were quantified using Qubit® ds -DNA HS Assay Kit and Agilent Bioanalyzer HS Kit. 383 Afterwards, the pools enriched with RNA baits were collectively combined into a single pool, while 384 those enriched with DNA baits were combined into another separate pool, each in equal DNA amounts. 385 This resulted in two distinct library pools. Both of these pools were sent for Illumina sequencing 386 (paired-end, 2 x 150 bp, MiSeq-v3) at Fasteris SA (Geneva, Switzerland). 387 388 Bioinformatics 389 The sequences were trimmed, mapped against references, and converted to finally calculate mapping 390 coverage and for variant calling. First, sequence file names were modified to fit the sample identifier 391 using python v3.10.13. Then raw sequence reads were processed using Autotrim v0.6.1, a tool that 392 automates the quality control and trimming of high -throughput sequence data 34. Autotrim employs 393 Trimmomatic v0.3935, to trim sequences based on quality and remove adapters, FastQC v0.12.1 36 to 394 assess the quality of the raw sequence data, and MultiQC v1.14 37 to aggregate the results into a single 395 report. Following the initial processing, the trimmed reads were mapped to the reference chloroplast 396 and mitochondrial genomes (GenBank: PRJNA493755) using BWA -MEM v0.7.1738. BWA is a software 397 package for mapping low-divergent sequences against a reference genome, which is particularly useful 398 for precisely aligning sequencing reads. 399 The mapped reads were then converted to bam format using Samtools v1.1539, a suite of programs for 400 interacting with high-throughput sequencing data. In addition, the per-sample mapping coverage was 401 calculated. After careful consideration, we decided against deduplicating our sequence data as 402 deduplication resulted in a significant loss of data. This is probably due to the nature of our short and 403 fragmented sequences. In addition, each variation, even those present in duplicates, may hold 404 biologically significant information. Therefore, removing duplicates may have eliminated critical 405 genetic diversity, compromising the integrity of our analysis. 406 VCFtools v0.1.1640 and BCFtools v1.1939 were used for calling and filtering variants from the alignment. 407 Variant calling is the process of identifying differences, such as single-nucleotide polymorphisms (SNPs) 408 and insertions/deletions (indels), between the sequenced sample and the reference genome. Variants 409 were filtered to include only those with a quality score (Q) of 30 or higher. The identified variants were 410 then phased into haplotypes using Beagle v5.441. Haplotyping, or phasing, is the process of determining 411 the distribution of alleles of multiple variants along each chromosome. In this context, it can be said it 412 accounts for unique combinations of alleles at variant sites (distinct variations) across the chloroplast 413 and mitochondrial genome. 414 The pattern of ancient DNA damage was assessed using mapDamage2 v2.2.2 42. This tool quantifies 415 patterns of DNA damage in ancient samples, which can provide insights into the preservation and 416 authenticity of ancient DNA. 417 418 Analyses 419 The bioinformatic output was converted into a table format and for further analyses with R v4.3.1: The 420 mapping coverage (i.e., the number of reads mapped against reference genome), allele frequencies, 421 distinct variants, nucleotide diversity and C -to-T-substitution rate were combined with the sample 422 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 14 metadata (Supplement, Tab. S3) using the R package “tidyverse” 43. The data was visualized using the 423 packages “tidyverse´”, “ggpubr” v0.6.044 and “tidypaleo” v0.1.345. 424 In order to facilitate a comparison of the genetic diversity between GOF and EGB, nucleotide diversity 425 (resulting from Samtools) and distinct variants (resulting from Beagle) were normalized based on the 426 mapping coverage for mitochondria and plastid. Second, Pearson's correlation was employed to verify 427 the consistency of nucleotide diversity between the mitochondrion and chloroplast. For comparison, 428 the dataset was subset to the overlapping time period of both locations. Subsequently, a t -test was 429 applied to compare and test the significance of the nucleotide diversity (mitochondrial and chloroplast 430 data combined) between the two sites. In addition, each distinct genetic variant was assigned a unique 431 identifier for easy reference. These identifiers were utilized to calculate two diversity indices: the 432 Shannon index and the Simpson index. These indices were compared between the EGB and GOF sites 433 using a t-test, 434 To further identify spatial differences between GOF and EGB, we analyzed allele frequencies (resulting 435 from VCFtools & BCFtools) and calculated the Fixation Index (FST) using corresponding samples from 436 both locations. The FST is a measure of population differentiation attributable to genetic structure. 437 The FST value ranges from 0 to 1, with 0 indicating complete interbreeding (i.e., the two populations 438 are freely intermixing) and 1 indicating that all genetic variation is explained by the population 439 structure (i.e., the two populations do not share any genetic diversity). We implemented a custom R 440 function to calculate the FST over time. This process involved sorting the data by age, subsetting 441 samples within a specified time window, and computing FST based on allele frequency variance. The 442 allelic composition, derived from parsing VCFtools output, was analyzed using Principal Component 443 Analysis (PCA) to examine the genetic structure.. We also performed Pearson’s correlation and 444 PERMANOVA to assess the impact of C-to-T substitution rates and Total Organic Carbon (TOC) on the 445 variation in the PCA space, specifically the first two principal components (PCs). These analyses helped 446 determine the extent to which these factors influenced genetic composition. For PCA and 447 PERMANOVA we also used the functions provided by the package “vegan” v2.6.446. 448 The C -to-T-substitutions at the first position from the mapDamage2 output were examined to 449 authenticate the ancient DNA (aDNA) and analyze the damage patterns. A Pearson's correlation test 450 was conducted to assess the influence of mapping coverage on these substitutions. This was done to 451 evaluate if coverage has an impact on the observed damage. The Baltic Sea Mn/Ti ratio data were also 452 integrated into the analysis. 453 To assess the genetic diversity and adaptability of the population, we calculated the rate of allele 454 turnover. This measure, derived from the ratio of significant allele changes (greater than 1%) to total 455 generations, reflects the estimated rate of genetic change over time based on the present data. Given 456 the low FST values, we inferred that the two sites, EGB and GOF, belong to a single population. The 457 calculation covered all of the data, from the oldest sample from EGB to the youngest sample from GOF, 458 providing insight into the evolving genetic makeup of the population. 459 To identify phases of change in the allele turnover across the entire dataset, a generalized additive 460 model (GAM) was fitted using the “gam” function from the “mgcv” package47. This model was applied 461 to the complete dataset, including periods where we did not have direct samples, as the allele turnover 462 was estimated based on the available data. We calculated the mean allele turnover for the present 463 data points, which were then included in the GAM plot. Allele turnover was used as a response variable 464 to fit a gam model as a function of “age”. 465 Corresponding metadata to this analysis can be found in supplementary material (Table S3). Data on 466 shipping activity in the Baltic Sea were obtained from three main sources: Kontny 48, Mägi 49, and 467 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 15 Rytkönen et al. 50. Kontny's work provided insights into maritime contacts during the Roman and 468 Migration Periods (1st -7th centuries AD), while Mägi's research focused on the role of the Eastern 469 Baltic in Viking Age communication across the Baltic Sea. Rytkönen et al.'s research on the statistical 470 analysis of Baltic Sea shipping revealed a significant increase in maritime traffic during the late 20th 471 and early 21st centuries. These sources were used to establish a timeline of shipping intensity from 472 1200 BP to the present. 473 474

475 This study was funded by the K314/2020 grant of the Collaborative Excellence Programme of the 476 Leibniz Association. The authors acknowledge support by the High Performance and Cloud Computing 477 Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden -478 Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 479 37/935-1 FUGG. AS acknowledges funding from the International Max Planck Research School for 480 Quantitative Behaviour, Ecology and Evolution. MB is supported by the Hessian Ministry of Higher 481 Education, Research and the Arts through the LOEWE Centre for Translational Biodiversity Genomics. 482 For help with preparing reagents for library preparation and enrichment experiments we thank Clara 483 Linde Schlimbach and Anna Weis. For setting up the ancient DNA clean laboratory in Frankfurt where 484 all DNA extractions were done by AS and JR, we thank Leonie Schardt. The authors declare no conflicts 485 of interest. 486 AI statement 487 In the preparation of this manuscript, we employed artificial intelligence (AI) tools to enhance the 488 quality of our work and to optimize our research process. 489 R Script Optimization: We utilized an AI -based code optimization to refine our R scripts. The tool 490 provided suggestions for code enhancement, identified potential bugs, and recommended more 491 efficient coding practices. This not only improved the performance of our scripts but also ensured the 492 reproducibility and reliability of our results. 493 We believe that the use of AI in our research process has improved the quality of our work. However, 494 we stress that the final decisions on manuscript content and the interpretation of the results were 495 made by the authors. 496 Data availability 497 The datasets generated and analyzed during this study are publicly available. Data and scripts can be 498 accessed through the GitHub repository at https://github.com/Alex-132/seda_enrichment. 499 Additionally, sequence data have been deposited in the European Nucleotide Archive (ENA) under 500 the accession number XXXXX (will be uploaded there). 501 502

503 1. Goulletquer, P., Gros, P., Boeuf, G. & Weber, J. The Impacts of Human Activities on Marine 504 Biodiversity. in Biodiversity in the Marine Environment 15–20 (Springer Netherlands, 505 Dordrecht, 2014). doi:10.1007/978-94-017-8566-2_2. 506 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 16 2. Panja, P., Kar, T. & Jana, D. K. Impacts of global warming on phytoplankton–zooplankton 507 dynamics: a modelling study. Environ Dev Sustain (2024) doi:10.1007/s10668-023-04430-3. 508 3. Benoiston, A.-S. et al. The evolution of diatoms and their biogeochemical functions. Phil. 509 Trans. R. Soc. B 372, 20160397 (2017). 510 4. Johansson, O. N. et al. Skeletonema marinoi as a new genetic model for marine chain-511 forming diatoms. Sci Rep 9, 5391 (2019). 512 5. Nakweya, G. Southern Ocean phytoplankton changes affect global climate regulation. Nat 513 Africa d44148-023-00242–9 (2023) doi:10.1038/d44148-023-00242-9. 514 6. Hjerne, O., Hajdu, S., Larsson, U., Downing, A. S. & Winder, M. Climate Driven Changes in 515 Timing, Composition and Magnitude of the Baltic Sea Phytoplankton Spring Bloom. 516 Frontiers in Marine Science 6, (2019). 517 7. Hattich, G. S. I. et al. Temperature optima of a natural diatom population increases as 518 global warming proceeds. Nature Climate Change 1–8 (2024). 519 8. Rynearson, T. A., Bishop, I. W. & Collins, S. The Population Genetics and Evolutionary 520 Potential of Diatoms. in The Molecular Life of Diatoms (eds. Falciatore, A. & Mock, T.) 29–57 521 (Springer International Publishing, Cham, 2022). doi:10.1007/978-3-030-92499-7_2. 522 9. Godhe, A. & Rynearson, T. The role of intraspecific variation in the ecological and 523 evolutionary success of diatoms in changing environments. Phil. Trans. R. Soc. B 372, 524 20160399 (2017). 525 10. Björck, S. The late Quaternary development of the Baltic Sea basin. in Assessment of 526 climate change for the Baltic Sea Basin 398–407 (Springer, 2008). 527 11. Zillén, L., Conley, D. J., Andrén, T., Andrén, E. & Björck, S. Past occurrences of hypoxia in the 528 Baltic Sea and the role of climate variability, environmental change and human impact. 529 Earth-Science Reviews 91, 77–92 (2008). 530 12. Ojaveer, H. et al. Status of Biodiversity in the Baltic Sea. PLoS ONE 5, e12467 (2010). 531 13. Reusch, T. B. H. et al. The Baltic Sea as a time machine for the future coastal ocean. Sci. 532 Adv. 4, eaar8195 (2018). 533 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 17 14. Smol, J. P. The power of the past: using sediments to track the effects of multiple stressors 534 on lake ecosystems. Freshwater Biology 55, 43–59 (2010). 535 15. Ellegaard, M. et al. Dead or alive: sediment DNA archives as tools for tracking aquatic 536 evolution and adaptation. Commun Biol 3, 1–11 (2020). 537 16. Rengefors, K., Kremp, A., Reusch, T. B. H. & Wood, A. M. Genetic diversity and evolution in 538 eukaryotic phytoplankton: revelations from population genetic studies. J. Plankton Res. 539 plankt;fbw098v1 (2017) doi:10.1093/plankt/fbw098. 540 17. Laso-Jadart, R., O’Malley, M., Sykulski, A. M., Ambroise, C. & Madoui, M.-A. Holistic view of 541 the seascape dynamics and environment impact on macro-scale genetic connectivity of 542 marine plankton populations. BMC Ecol Evo 23, 46 (2023). 543 18. Brand, L. E. Review of Genetic Variation in Marine Phytoplankton Species and the Ecological 544 Implications. Biological Oceanography (1980) doi:10.1080/01965581.1988.1074954. 545 19. Salmaso, N. & Tolotti, M. Phytoplankton and anthropogenic changes in pelagic 546 environments. Hydrobiologia 848, 251–284 (2021). 547 20. Myrberg, K. & Soomere, T. The Gulf of Finland, Its Hydrography and Circulation Dynamics. in 548 Preventive Methods for Coastal Protection: Towards the Use of Ocean Dynamics for 549 Pollution Control (eds. Soomere, T. & Quak, E.) 181–222 (Springer International Publishing, 550 Heidelberg, 2013). doi:10.1007/978-3-319-00440-2_6. 551 21. Leppäranta, M. History and Future of Snow and Sea Ice in the Baltic Sea. in Oxford Research 552 Encyclopedia of Climate Science (2023). doi:10.1093/acrefore/9780190228620.013.891. 553 22. Salimi, P. A. et al. A review of the diversity and impact of invasive non-native species in 554 tropical marine ecosystems. Mar Biodivers Rec 14, 11 (2021). 555 23. Wang, Q., Cheng, F., Xue, J., Xiao, N. & Wu, H. Bacterial community composition and 556 diversity in the ballast water of container ships arriving at Yangshan Port, Shanghai, China. 557 Mar Pollut Bull 160, 111640 (2020). 558 24. Andrés, J. et al. Environment and shipping drive environmental DNA beta‐diversity among 559 commercial ports. Molecular Ecology 32, 6696–6709 (2023). 560 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 18 25. Wang, Y., Wessels, M., Pedersen, M. W. & Epp, L. S. Spatial distribution of sedimentary DNA 561 is taxon-specific and linked to local occurrence at intra-lake scale. Communications Earth 562 & Environment 4, 172 (2023). 563 26. Carstensen, J. et al. Hypoxia in the Baltic Sea: Biogeochemical Cycles, Benthic Fauna, and 564 Management. AMBIO 43, 26–36 (2014). 565 27. Snoeijs-Leijonmalm, P. & Andrén, E. Why is the Baltic Sea so special to live in? in Biological 566 Oceanography of the Baltic Sea (eds. Snoeijs-Leijonmalm, P., Schubert, H. & Radziejewska, 567 T.) 23–84 (Springer Netherlands, Dordrecht, 2017). doi:10.1007/978-94-007-0668-2_2. 568 28. Andrén, T. et al. The Development of the Baltic Sea Basin During the Last 130 ka. in The 569 Baltic Sea Basin (eds. Harff, J., Björck, S. & Hoth, P.) 75–97 (Springer Berlin Heidelberg, 570 Berlin, Heidelberg, 2011). doi:10.1007/978-3-642-17220-5_4. 571 29. Kabel, K. et al. Impact of climate change on the Baltic Sea ecosystem over the past 1,000 572 years. Nature Climate Change 2, 871–874 (2012). 573 30. Seppä, H., Bjune, A. E., Telford, R. J., Birks, H. J. B. & Veski, S. Last nine-thousand years of 574 temperature variability in Northern Europe. Clim. Past 5, 523–535 (2009). 575 31. Epp, L. S., Zimmermann, H. H. & Stoof-Leichsenring, K. R. Sampling and Extraction of 576 Ancient DNA from Sediments. in Ancient DNA (eds. Shapiro, B. et al.) vol. 1963 31–44 577 (Springer New York, New York, NY, 2019). 578 32. Schmidt, A. et al. Decoding the Baltic Sea’s past and present: A simple molecular index for 579 ecosystem assessment. Ecological Indicators 166, 112494 (2024). 580 33. Gansauge, M.-T., Aximu-Petri, A., Nagel, S. & Meyer, M. Manual and automated preparation 581 of single-stranded DNA libraries for the sequencing of DNA from ancient biological remains 582 and other sources of highly degraded DNA. Nat Protoc 15, 2279–2300 (2020). 583 34. Waldvogel, A.-M. et al. The genomic footprint of climate adaptation in Chironomus riparius. 584 Molecular Ecology 27, 1439–1456 (2018). 585 35. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence 586 data. Bioinformatics 30, 2114–2120 (2014). 587 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 19 36. Andrews, S. FastQC: A quality control analysis tool for high throughput sequencing data. 588 (2021). 589 37. Ewels, P., Magnusson, M., Lundin, S. & Käller, M. MultiQC: summarize analysis results for 590 multiple tools and samples in a single report. Bioinformatics 32, 3047–3048 (2016). 591 38. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows–Wheeler transform. 592 Bioinformatics 25, 1754–1760 (2009). 593 39. Danecek, P. et al. Twelve years of SAMtools and BCFtools. GigaScience 10, giab008 (2021). 594 40. Danecek, P. et al. The variant call format and VCFtools. Bioinformatics 27, 2156–2158 595 (2011). 596 41. Browning, B. L., Tian, X., Zhou, Y. & Browning, S. R. Fast two-stage phasing of large-scale 597 sequence data. The American Journal of Human Genetics 108, 1880–1890 (2021). 598 42. Jónsson, H., Ginolhac, A., Schubert, M., Johnson, P. L. F. & Orlando, L. mapDamage2.0: fast 599 approximate Bayesian estimates of ancient DNA damage parameters. Bioinformatics 29, 600 1682–1684 (2013). 601 43. Wickham, H. et al. Welcome to the Tidyverse. JOSS 4, 1686 (2019). 602 44. Kassambara, A. ggpubr: ‘ggplot2’ Based Publication Ready Plots. (2023). 603 45. Dunnington, D. W., Libera, N., Kurek, J., Spooner, I. S. & Gagnon, G. A. tidypaleo : 604 Visualizing Paleoenvironmental Archives Using ggplot2. J. Stat. Soft. 101, (2022). 605 46. Oksanen, J. et al. vegan: Community Ecology Package. R package version 2.5-6. 2019. 606 Community Ecol 8, 732–740 (2020). 607 47. Wood, S. N. Fast Stable Restricted Maximum Likelihood and Marginal Likelihood Estimation 608 of Semiparametric Generalized Linear Models. Journal of the Royal Statistical Society Series 609 B: Statistical Methodology 73, 3–36 (2011). 610 48. Kontny, B. Maritime contacts across the Baltic Sea during the Roman and Migration Periods 611 (1st-7th centuries AD) in the light of archaeological sources from the central-European 612 perspective. (2023). 613 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint 20 49. Mägi, M. In Austrvegr: The Role of the Eastern Baltic in Viking Age Communication across the 614 Baltic Sea. Brill; Leiden. vol. 84 (2018). 615 50. Rytkönen, J., Siitonen, L., Riipi, T., Sassi, J. & Sukselainen, J. Statistical Analyses of the Baltic 616 Maritime Traffic. (VTT Technical Research Centre of Finland, 2002). 617 618 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted March 13, 2025. ; https://doi.org/10.1101/2025.03.10.642313doi: bioRxiv preprint