A 2,000‐Year Record of Eelgrass (Zostera marina L.) Colonization Shows Substantial Gains in Blue Carbon Storage and Nutrient Retention

Assessing historical environmental conditions linked to habitat colonization is important for understanding long‐term resilience and improving conservation and restoration efforts. Such information is lacking for the seagrass Zostera marina, an important foundation species across cold‐temperate coastal areas of the Northern Hemisphere. Here, we reconstructed environmental conditions during the last 14,000 years from sediment cores in two eelgrass (Z. marina) meadows along the Swedish west coast, with the main aims to identify the time frame of seagrass colonization and describe subsequent biogeochemical changes following establishment. Based on vegetation proxies (lipid biomarkers), eelgrass colonization occurred about 2,000 years ago after geomorphological changes that resulted in a shallow, sheltered environment favoring seagrass growth. Seagrass establishment led to up to 20‐ and 24‐fold increases in sedimentary carbon and nitrogen accumulation rates, respectively. This demonstrates the capacity of seagrasses as efficient ecosystem engineers and their role in global change mitigation and adaptation through CO2 removal, and nutrient and sediment retention. By combining regional climate projections and landscape models, we assessed potential climate change effects on seagrass growth, productivity and distribution until 2100. These predictions showed that seagrass meadows are mostly at risk from increased sedimentation and hydrodynamic changes, while the impact from sea level rise alone might be of less importance in the studied area. This study showcases the positive feedback between seagrass colonization and environmental conditions, which holds promise for successful conservation and restoration efforts aimed at supporting climate change mitigation and adaptation, and the provision of several other crucial ecosystem services.

Plain Language Summary This study investigated the historical colonization of eelgrass (Zostera marina), an important marine vascular plant in cold-temperate coastal regions.Sediment cores from eelgrass meadows on the Swedish west coast dating back as far as 14,000 years were examined to understand the timecourse of eelgrass colonization and the subsequent modification of the environment.We found that eelgrass colonization began approximately 2,000 years ago, coinciding with the development of shallow, sheltered conditions that favored eelgrass growth.As eelgrass became established, this led to substantial habitat and sediment changes, including a 20-and 24-fold increase in carbon and nitrogen accumulation, respectively.This highlights the crucial role of eelgrass as a provider of important ecosystem services, such as regulation of climate, nutrient retention, and sediment protection.We also examined the potential effects of climate change on eelgrass growth and health, predicting that decreased water clarity and altered water flow pose the greatest risks.

Introduction
Coastal vegetated habitats, including seagrass meadows, are highly productive and diverse marine ecosystems (Barbier et al., 2011).Seagrass meadows provide a myriad of ecosystem services, such as biodiversity maintenance, provision of nursery habitats, contribution to climate change mitigation and coastal protection, and water purification (Cullen- Unsworth & Unsworth, 2013;de los Santos et al., 2020;Macreadie et al., 2021).Thus, seagrasses contribute important ecosystem services for the wellbeing of people and the planet.Through the reduction of hydrodynamic forces within the seagrass canopy and stabilization of sediment within the rootrhizome system, seagrasses trap and embed organic and inorganic particles (Lei et al., 2023;Samper-Villarreal et al., 2016) of both allochthonous and autochthonous origin (Asplund et al., 2021;Kennedy et al., 2010;Oreska et al., 2018) leading to accumulation of thick organic-rich sediment deposits, which can remain stable over millennia (Mateo et al., 1997).Through continuous sediment accretion, seagrass meadows support long-term blue carbon storage, nutrient retention and reduced turbidity (through sediment stabilization) in coastal areas (Lima et al., 2020;Mazarrasa et al., 2018).
Eelgrass (Zostera marina L.) is one of the most prevalent seagrass species globally with a geographical distribution spanning the entire cold-temperate zone of the Northern hemisphere (Short et al., 2007).The colonization of seagrasses in cold-temperate regions is highly dependent on the attenuation of irradiance with increasing water depth, and the lower depth limit of Z. marina distribution is usually regulated by water clarity and turbidity (Krause-Jensen et al., 2011;Nielsen et al., 2002), which is often negatively impacted by eutrophication, coastal exploitation (e.g., sedimentation) and climate change (through increased precipitation and land runoff) (Rabalais et al., 2009).The colonization, growth and distribution of Z. marina are governed by environmental conditions, such as hydrodynamic-driven turbidity and water depth, and human-induced disturbances, affecting water quality.To understand the historical establishment and subsequent colonization of seagrass plants, paleoreconstruction of sedimentary records has been successfully performed for seagrass meadows (Serrano et al., 2020), providing insights into the dynamics of ecosystem health status in response to long-term environmental change (Leiva-Dueñas et al., 2021;López-Sáez et al., 2009;Mateo et al., 2010).As seagrass loss and recolonization is related to environmental dynamics, information on past environmental conditions linked to seagrass ecosystem integrity could be of vital help for coastal managers guiding seagrass restoration efforts as well as to predict impacts on seagrass distribution and provisioning of ecosystem services from future climate change.
Climate change is considered a major threat to seagrass habitats (Short & Neckles, 1999), potentially leading to dramatic changes in the coastal environment, including loss of seagrass meadows (Orth et al., 2006).Global temperatures have drastically increased during the twenty-first century (IPCC, 2021), which entails the thermal expansion of the oceans and melting of land ice.Therefore, sea levels are predicted to increase globally (Masson-Delmotte et al., 2021), although regionally the rising sea levels could offset local sea level regressions from ongoing isostatic uplift (Meier et al., 2004).In addition, precipitation in the cold-temperate zone has increased and even higher precipitation is expected, with stronger rainfall events in the near future (Pörtner et al., 2022), which can lead to decreased salinity in coastal waters due to higher freshwater runoff from land (Cheng et al., 2020) that may in turn negatively influence seagrass distribution (Stevens & Lacy, 2012).
For cold-temperate seagrass species, such as Z. marina, using paleoreconstructions for describing past environmental conditions and to predict future changes is lacking (but see Kindeberg et al., 2019;Leiva-Dueñas et al., 2023).There are a number of seagrass-directed paleorecords of the genus Posidonia (López-Merino et al., 2017;Macreadie et al., 2012;Mateo et al., 1997Mateo et al., , 2010;;Serrano et al., 2011Serrano et al., , 2016)), but these are geographically constrained to Australia and the Mediterranean.Therefore, understanding the establishment and colonization of the widespread Z. marina is of high relevance and can likely apply to other seagrass species growing under similar environmental conditions.Here we aim to reconstruct millenary changes in habitat condition in response to environmental changes based on biogeochemical analyses of two sediment cores from Z. marina meadows to (a) establish a time frame for seagrass colonization, (b) assess the climatic and geomorphological conditions that favored seagrass colonization, (c) describe the biogeochemical changes in the sediment following the establishment of the seagrass meadows, and (d) predict potential future climate change impacts (from sea level rise [SLR], decrease in salinity and increase in runoff and sea surface temperature [SST]) on the seagrass meadow health status.

Study Context
The Swedish coastal zone has been severely impacted by eutrophication and sediment supply from land over the last century, resulting in a recession in the depth limit distribution of Z. marina (Boström et al., 2014;Rosenberg et al., 1996).Historical reports for Z. marina on the Swedish west coast showed depths of 15-20 m (Loo & Isaksson, 2015;Petersen, 1893), while today eelgrass plants are rarely found below 9 m (Skåne County Administrative Board, 2016; pers.obs.Asplund and Svensson) and generally growing in water depths from 0.5 to 6 m (Svensson et al., 2021).In addition to light conditions and turbidity levels, seawater temperature and salinity can negatively affect Z. marina performance (Marsh et al., 1986;Rasmusson et al., 2020;Xu et al., 2016).
Along the Swedish coast, fluctuations in sea level, temperature and salinity have occurred over the Holocene (the last ∼11.7 cal ka BP) due to isostatic and eustatic changes following the retreat of the Scandinavian ice sheet (Berner et al., 2011;Mörner, 1969;Törnqvist & Hijma, 2012) leading to geomorphological transformations along the coastline.Warming in the early Holocene caused the ice sheet to retreat and the westward outlet opened the drainage of the Baltic ice lake at 11.7 cal ka BP (Andrén et al., 2011;Björck, 2008;Jakobsson et al., 2007).This was followed by a series of dynamic changes in hydrographic conditions affecting sedimentation (Gyllencreutz et al., 2006) and at about 8.5 cal ka BP, there was a shift in the ocean circulation resulting in a higher inflow of marine waters (Conradsen & Heier-Nielsen, 1995).During the last ∼2 cal ka BP, the Swedish Skagerrak coast has been dominated by Atlantic water conditions (Erbs-Hansen et al., 2012).Periods of increased temperatures occurred during the Roman Warm Period (RWP) (∼2.3 to 1.6 cal ka BP) (Hass, 1996;Neukom et al., 2019) and the Medieval Climate Anomaly (MCA) (∼1.2 to 0.9 cal ka BP) (Mann et al., 2009), and these warmer periods were linked to less frequent storms and lower sedimentation runoff in the Skagerrak region (Hass, 1996).

Study Area and Sediment Core Sampling
Field sampling was conducted in two Z. marina meadows at Gåsö (58°14′10″N, 11°23′40″E), situated in the vicinity of the Gullmar Fjord mouth on the Swedish Skagerrak coast (Figure 1a), in August 2020.Gåsö is a small island (∼3 km 2 ) mainly composed of granitic bedrock (Figure 1b).In the middle of the island, connecting the western and the eastern parts, there is a low-elevation (less than 1 m above sea level [m a. s. l.]) postglacial sand deposit.The island has a relatively low exposure to human activities with only a few country residences in the eastern part.On the southern coast of the island, there is an up to 3 m deep elongated cove stretching toward the center of the island, whereas on the northern side there is a shallow (up to 4 m depth) embayment surrounded by smaller rocks and islets creating a hydrodynamically sheltered environment.The main wind directions are from the W and SW.Within both the northern and southern embayments, seagrass meadows are found growing from about 0.5 to 4 m water depth.The distribution of eelgrass was determined by encircling the meadows by boat repeatedly crossing the edge of the meadow with support of sonar and a sub-surface camera for validation of Z. marina vegetation.GPS waypoints were collected at the intersections and corrected for geometry errors caused by drifting prior to vectorization.The mapped seagrass area was approximately 7 and 5 ha at Gåsö S and N, respectively.In the middle of each of the seagrass meadows at Gåsö South (S) and Gåsö North (N), a sediment core (2 m long and 7.5 cm in diameter) was collected using SCUBA by hammering a PVC-corer.The sediment cores in Gåsö S and N were collected at 1.8 and 3 m water depths, respectively (Figure 1c), and the top ∼50 cm of the sediment cores have previously been described by Dahl et al. (2023).Sediment core compression was assessed once during coring by measuring the inner and outer portions of the corer before core extrusion from the seagrass meadow.From these values, a compression factor was calculated, which was similar for the two cores (18% and 17% in Gåsö S and N, respectively).The compression factors were applied to sediment depth and dry bulk density (DBD) to calculate decompressed values, which were used in subsequent analysis.The cores were stored in a cold room and were cut lengthwise into two hemi-cores.One half was sliced into 1 cm intervals and the other half was kept intact for X-ray fluorescence (XRF) measurements.

Sediment Chronology
The concentrations of 210 Pb were analyzed by measuring its decay product ( 210 Po) in equilibrium using alpha spectrometry (Sanchez-Cabeza et al., 1998), which aimed to obtain age models for the two cores.Details of the 210 Pb analysis are presented in Dahl et al. (2023).Ten radiocarbon ( 14 C) dates were obtained along the cores (6 in Gåsö N and 4 in Gåsö S).Prior to 14 C measurements in a single-stage Accelerator Mass Spectrometer (AMS), samples were digested using the standard acid-base-acid method (Molnár et al., 2013).The background values of measurements were estimated to be 0.25 pMC using phthalic anhydride with the NIST-OXII (134.06 pMC) standard used as reference material.The 14 C/ 12 C ratio was measured with an accuracy of >0.3%.For the isotopic fractionation correction, the ratio of 13 C/ 12 C was used.The 14 C ages were calibrated using the Marine20 calibration curve (Heaton et al., 2020) and corrected for a local marine reservoir effect (ΔR = 208 ± 57 years) based on map nos.77, 681-684, and 674 (Håkansson, 1970(Håkansson, , 1987;;Olsson, 1980) in the Marine20 reservoir database.Core age models based on both 210 Pb-and 14 C-dates were estimated using the R-package Bacon (Blaauw & Christeny, 2011).The sediment accumulation rates (SAR) and mass accumulation rates (MAR) were calculated using the Constant Rate of Supply (CRS) models for the top 210 Pb-dated sediment cores, and from the Bacon models based on cumulative mass versus age (Belshe et al., 2019) for the older sediment layers.SAR was later calculated by dividing the Bacon-derived MAR by the sediment DBD.

Sediment Biogeochemical Analyses
All sediment slices (n = 166 and 148 for Gåsö N and S, respectively) were dried at 60°C until constant weight was obtained to estimate DBD (g cm 3 ) (see data set [DOI: 10.5281/zenodo.10559287]for details on biogeochemical analysis performed along the cores).Samples (n = 21 and 22 for Gåsö N and S, respectively) for organic carbon (OC) and nitrogen (N) contents and C and N stable isotopic composition (δ 13 C and δ 15 N) were analyzed using a Carlo Erba NC2500 elemental analyzer connected to a Thermo Scientific Delta V Advantage Isotope Ratio Mass Spectrometer (EA-IRMS).Acetanilid (δ 13 C = 26.14 ± 0.15‰, δ 15 N = 0.38 ± 0.12‰) was used as the reference material and the standard deviation was approximately 1% for C and N concentrations and 0.1‰ for both C and N stable isotope ratios.The isotopic compositions were expressed in delta notation (per mil) relative to the VPDB (Vienna PeeDee Belemnite) for δ 13 C and to the atmospheric nitrogen standard for δ 15 N. Prior to the OC and N analysis, the sediment was ground using a mixer mill (Retsch MM 400), weighted in silver capsules and treated with 1 M HCl to remove inorganic carbon (direct addition using a pipet) (Dahl et al., 2016).Carbon accumulation rates (CAR) and nitrogen accumulation rates (NAR) were calculated using the weighted mean of OC and N contents and MAR (Arias-Ortiz et al., 2020).Calcium carbonate (CaCO₃) contents were estimated from an aliquot of dry sediment (n = 28 and 14 for Gåsö N and S, respectively) through loss on ignition combustion at 450°C for 6 hr followed by 2 hr at 950°C following Bengtsson and Enell (1986) and Heiri et al. (2001).Grain size distribution following the classification of Wentworth (1922) was performed using a laser diffraction particle size analyzer (Mastersizer 2000 MALVERN).Prior to analysis, the samples (n = 28 and 24 for Gåsö N and S, respectively) were sieved by 2 mm and the fraction <2 mm was treated with 30% H₂O₂ to remove organic matter.Mean grain size (μm) and sorting coefficient were calculated using the Gradistats program (Blott & Pye, 2001) based on Folk and Ward (1957).
Lipid analysis was carried out to trace seagrass plant-derived material in the sediment record (Jaffé et al., 2001;Sikes et al., 2009).For the analysis, 14 samples from each core as well as above-(n = 3) and belowground biomass (n = 3) of seagrass were freeze-dried and milled.A mixture of dichloromethane (DCM) and methanol (MeOH) (9:1 v/v) was added to the dried sediment, and samples were placed in a sonic bath for 15 min (Poynter & Eglinton, 1990).The samples were then centrifuged at 800 rpm for 10 min and the lipid extract was placed in a glass tube.This was repeated three times and the lipid extracts were combined.Activated copper powder was added to the lipid extracts for removal of elemental sulfur, and the lipid extracts were dried under a N 2 blowdown system and later re-dissolved in <1 mL DCM.Deactivated 95% silica gel was added to adsorb the lipid extract and the adsorbed samples were placed in glass pipettes that had been packed with pre-combusted (400°C, overnight) and deactivated (addition of 5% H 2 O by weight) silica gel, and the fractions of non-polar hydrocarbons were eluted using hexane.These fractions were dissolved in 400-1,000 μL hexane depending on the expected concentration on a Shimadzu QP2010 Ultra GC-MS equipped with an AOC-20i auto sampler and a split-splitless injector operated in splitless mode.A Zebron ZB-5HT Inferno GC column (30 m × 0.25 mm × 0.25 μm) was used for separation.An external standard containing a mixture of C 20 -C 40 n-alkanes with known concentration was analyzed in conjunction with the samples and used for quantification based on peak areas.The average chain length (ACL) (for homologs between C 17 and C 35 ) was calculated following ACL = ∑ (n × C n )/∑ C n , where n is the number of carbon atoms and C n is the concentration (mg g DW 1 ) of the n-alkane.
Magnetic susceptibility (MS) was used as a proxy to assess the abundance of ferromagnetic minerals derived from terrestrial soil erosion to assess the relative relationship between marine and terrestrial influences on sedimentation (López-Merino et al., 2017).MS was analyzed in 37 and 40 samples from Gåsö N and S, respectively, using a Multi Sensor Core Logger (LABCORE, University of Barcelona).Samples were packed in a 7 cm 3 cubicle and analyzed with a low frequency (∼0.1 mT) volumetric MS.Each sample was analyzed five times and calibration was manually performed before starting each measurement.Mass specific MS (χ; cm 3 g 1 ) was calculated as χ = κ/ρ, where κ is the average of the replicates and ρ is sample density (cm 3 g 1 ) (Hatfield et al., 2013).
XRF measurements were performed on the intact hemi-cores using an ITRAX XRF core scanner from Cox Analytical Systems, which produces digital imagery of the cores and μ-XRF elemental profiles (Si, S, Cl, K, Ca, Ti, Mn, Fe, Ni, Cu, Zn, Br, Rb, Sr, Zr and Pb).The core scanner used a Mo tube with the setting of 30 kV and 50 mA, and a step size of 500 μm with a dwell time of 25 s.Before data handling, the μ-XRF elemental profiles were normalized using centered log-ratio transformation, which also accounts for non-linear effects of the elemental matrix (Weltje & Tjallingii, 2008).

Climate Modeling
Future changes in sea level were extracted from the IPCC AR6 SLR Projections (Fox-Kemper et al., 2021;Garner et al., 2021).The medium-confidence projections of the integrated SLR over all components, including the Antarctic and Greenland ice sheets, glaciers, land water storage, ocean dynamics (including thermal expansion), and vertical land motion (non-climatic processes), were utilized.Although the projections were available for the period CE 2020-2150 at a 10-year resolution, with the baseline period of 1995-2014, we limited our projection to 2100.Historical simulations and future projections of SST, sea surface salinity, and surface runoff were derived from global climate models (GCMs) in the Coupled Model Intercomparison Project Phase 6 (CMIP6) (Eyring et al., 2016) (Table S1 in Supporting Information S1).The historical simulations cover the period 1850-2014 and all forcings were included.The future projections for the period 2015-2100 are from the Scenario Model Intercomparison Project (ScenarioMIP) for CMIP6 (O'Neill et al., 2016), in which future greenhouse gas emission scenarios are derived from the Shared Socioeconomic Pathways (SSPs) with different climate policies.Projections from the four SSPs in Tier 1 of the ScenarioMIP for CMIP6 are adopted in the current work, namely SSP5-8.5,SSP3-7.0,SSP2-4.5, and SSP1-2.6.The four SSPs cover a major range of emission scenarios, from sustainability to fossil-fueled development.The time series of sea-level changes and the three selected variables are from the GCM grid closest to the study area.

GIS Analyses
Coastal geomorphology reconstruction based on eustatic sea level change and isostatic change for the latest 11 cal ka BP was analyzed using a shore displacement model from the Geological Survey of Sweden (SGU) (freely available at https://www.sgu.se/produkter-och-tjanster/geologiska-data/oppna-data/jordarter-oppna-data/strandforskjutningsmodell/) based on calculations from Påsse and Andersson (2005) and Påsse and Daniels (2015).Data from the model were extracted for Gåsö and analyzed in ArcGIS pro (v.2.9).From this data, bay exposure in terms of effective fetch (km) was calculated for the seagrass sites following the calculations outlined in Rogala (1997) by creating sight lines each 12°for a full circle (360°) with a maximum distance of the lines set to 3 km.The effective fetch was calculated for every 1,000 years over the past 11 cal ka BP.Based on the predicted SLR for the four different climate change scenarios (i.e., SSP5-8.5,SSP3-7.0,SSP2-4.5, and SSP1-2.6), the future sea levels and the resulting changes in shoreline were modeled in ArcGIS pro and calculated from digital elevation models (based on 0.1 altitude resolution LiDAR data from the Swedish National Land Survey) using the raster calculator function.Based on the predicted future sea levels, the effective fetch was also calculated for the different climate change scenarios.The water depth curves of Gåsö were extracted from publicly available nautical charts obtained from the Swedish Maritime Administration.

Statistical Analyses
All statistical analyses were performed in R (v. 4.2.2).Principal component analysis (PCA) was used to explore common patterns in the different biogeochemical properties (i.e., OC, C/N, δ 13 C, δ 15 N, Br, Ca, Cl, Fe, K, Mn, Rb, Sr, Ti, Zr, MS, ACL, mean grain size and sediment sorting) analyzed along the cores.Prior to the PCA, the missMDA-package was used to impute missing values in the sediment profiles as not all biogeochemical analyses were performed in each sediment layer.We assessed the uncertainty associated with the data included in the PCA through multiple imputations with the function MIPCA.The results showed that the over-imputed values were acceptable and the PCA results were worth interpreting (Figure S1 in Supporting Information S1).The first two principal components (PC1 and PC2) explained a high proportion of the variance but showed a "horseshoe effect" that is an artifact typical of PCA when the PC-curves correlate to each other (Goodall, 1954), which was addressed by calculating the arc-length distance of PC1 and PC2.Finally, change point modeling (CPM) was used to identify change points along the arc-length distance profiles using the Beast-package (Zhao et al., 2019).

Sediment Age Models
The mean (±SD) calibrated radiocarbon age at the bottom of the cores were 14.6 ± 0.2 (188 cm) and 12.3 ± 0.1 (165 cm) cal ka BP for Gåsö S and N, respectively.The mean (min-max) MAR was 0.0101 g cm 2 year 1 (0.0094-0.0107) and SAR was 0.0151 cm year 1 (0.0140-0.0159) for Gåsö S and MAR was 0.0060 g cm 2 year 1 (0.0057-0.0062) and SAR was 0.0114 cm year 1 (0.0108-0.0118) for Gåsö N (Figure S2 and Table S2 in Supporting Information S1).For Gåsö S, a large stone was encountered in the sediment core at a cumulative mass of around 24 g cm 2 (corresponding to 70 cm sediment depth; Figure 3) and this section was therefore excluded in the Bacon model.In Gåsö N, one 14 C date was omitted in the age model due to inversed age (Table S2 in Supporting Information S1).

Stratigraphy and Lithological Units
The bottom lithological units in both Gåsö S and N cores were composed of homogeneous brown clay (from 195 to 140 cm in Gåsö S and from 173 to 141 cm in Gåsö N) (Figure 2).In Gåsö S, brownish clay with shell fragments and pebbles was found from 140 to 104 cm, whereas homogeneous gray clay (without shell fragments) was found in Gåsö N from 141 to 106 cm.A layer of gray clay with shell fragments was found from 106 to 69 cm, which was also present in Gåsö S (104-85 cm).A 7-12 cm layer of gravel mixed with shells and clay was found between 72 and 60 cm in Gåsö S and between 69 and 60 cm in Gåsö N. The uppermost section in both cores (from ∼60 cm to the surface), corresponding to the latest ∼2.0 cal ka BP, contained organic-rich and dark brown sediment with seagrass rhizomes (Figure 2).S1 in Supporting Information S1).

Physical and Biogeochemical Sediment Properties
The oldest sediment layers (14-8 cal ka BP) had relatively high MS (80-100 χ cm 3 g 1 ) associated with lithogenic material and δ 13 C values below 20‰ indicating terrestrial origin, which transitioned toward more marine signatures in shallower sediment layers (Figure 3a).At the bottom of the cores, the sediment was finer grained and well sorted between 14 and 8 cal ka BP and shifted toward coarser and more poorly sorted after 8 cal ka BP (Figure 3a), likely reflecting changes in the source material and deposition of finer grain-sized material (Figure 3).In the mid-section of the sediment profiles, the sediment was to a large degree composed of biogenic carbonate shells, and the highest carbonate content was found between 8 and 4.5 cal ka BP (Figures 3 and 4a).For the last 2 cal ka BP (from 60 to 62 cm to the surface), there was a change in the biogeochemical profiles with higher OC, N and silt-clay content, and lower DBD (Figure 3b), which was noticeable during the last ∼1 cal ka BP (Figure 4b).From ∼2 cal ka BP until present, the n-alkane depth profiles showed a higher contribution of C 17 , C 19 and C 21 homologs, which are characteristic of Z. marina (Figure S3 in Supporting Information S1).The introduction of marine macrophytes was also reflected in an initial increase in average chain length, ACL, of the n-alkanes (Figure 3b).

Principal Component Analysis and Change Point Modeling
Two principal components explained 58% of the variability for which PC1 accounted for 36% and PC2 for 22% of the variance (Figure 4a).PC1 showed positive loadings (>0.6) for the XRF-elements Rb, Fe, K, Mn and Ti, and negative loadings (< 0.6) for δ 13 C, OC, Cl and Br, whereas PC2 showed positive loadings (>0.6) for OC, Br and Zr, and negative loadings (< 0.6) for CN, Sr and Ca (Table 1).The relationship of different sediment properties varied along the cores with strong negative correlations between sediment depth and OC, δ 13 C and Br, respectively, while elements such as Fe, Ti and MS were positively correlated to sediment depth (Figure 4a).CPM applied to the arc-length distance profiles for the two sites identified three change points in the core from Gåsö S at ∼0.7 (48 cm), ∼2.7 (68 cm) and ∼10.2 cal ka BP (164 cm).At Gåsö N, five change points were suggested at ∼1.1 (52 cm), ∼2.6 shows the different shapes of the PC curves for the two cores, which is due to differences in sediment composition as seen in the PCA biplot (lower plot in a).The PCA biplot presents the relationship of sediment characteristics for PC1 and PC2, with colors indicating core depth.The arc-length distance in panel (b) represents the distance of PC1 and PC2 to the modelled curve, and was applied due to correlation of the two PC-axes.The black line in panel (b) indicates the first appearance of seagrass-derived sediment in the paleorecord.Section 1 shows the time frame prior to seagrass establishment (as a baseline), whereas Section 2 presents the period around seagrass colonization and Section 3 presents the seagrass stabilization phase (see Table 2).
(68 cm), ∼8.6 (128 cm), ∼9.6 (134 cm) and ∼11.5 cal ka BP (146 cm) (Figure 4b).Based on identified change points of the arc-length distance of PC1 and PC2, three main sections were distinguished that reflect the time frame of main interest for seagrass colonization (Figure 4 and Table 2).Section 1 is the time period prior to seagrass colonization, Section 2 marks the period for seagrass colonization, and Section 3 is defined as the seagrass stabilization period.Section 3 had up to 6-fold higher SAR than Section 1, whereas MAR was about 30%-50% faster in Section 3 than in Section 1 (Table 2).The CAR and NAR were also up to 24-fold higher in Section 1 (CAR: 8.2 ± 1.2 and 5.5 ± 1.7 g OC m 2 year 1 , NAR: 0.75 ± 0.11 and 0.50 ± 0.31 g N m 2 year 1 ) compared to Section 3 (CAR: 0.5 ± 0.1 and 0.3 ± 0.1 g OC m 2 year 1 , NAR: 0.05 ± 0.01 and 0.02 ± 0.04 g N m 2 year 1 ).

Coastal Geomorphology and Isostatic Changes During Holocene
Based on the GIS (Geographic Information Systems) analysis, the effective fetch (km) for the Gåsö sites decreased during the latest 6 cal ka BP due to the isostatic uplift (Figure 5a), and in particular over the latest 2 cal ka BP when Gåsö was above the sea level creating a relatively sheltered environment (with 58%-72% lower effective fetch compared to the conditions prevalent about 7-11 cal ka BP).Over the last millennium, the effective fetch and shore-level have stabilized, and Gåsö has a similar coastal geomorphology as today (Figure 5b).

Modeling the Impacts of Climate Change to 2100
The climate models predict a sea-level increase ranging from 0.20 to 0.53 m for the different SSP scenarios by 2100, which could result in the immersion of the center of the Gåsö Island where postglacial sands are deposited (Figures 6a and 6b).In the high emission scenarios (SSP3-7.0 and SSP5-8.5), the SLR will create an opening in the middle of the island by 2100 and potentially change the hydrodynamics of the southern embayment.The predictions revealed a maximum of 3% increase in effective fetch by 2100 due to SLR (based on SSP5-8.5).The SST is predicted to increase in all SSP scenarios by 1.3-4.5°C(mean over the year) in 2100 compared to the present SST Note.Bold numbers show factor loadings of >0.6 for positive loadings and < 0.6 for negative loadings for PC1 and PC2.The data included in the PCA were organic carbon (OC), C:N-ratio (CN) and stable isotopes of C (δ 13 C) and N (δ 15 N), XRF-elements (Cl, Br, Zr, Sr, Ca, Rb, Fe, K, Mn and Ti), n-alkanes (average chain length, ACL), mean sediment grain-size distribution and sorting, and magnetic susceptibility (MS).(calculated as an average of the SST for 2004-2014) (Figure 7a).The salinity will decrease in all SSP scenarios from around 27-28 at present down to 25-26 (Figure 7b).However, large variation and no clear trends were seen in the regional models for surface runoff (Figure 7c).

Discussion
The reconstruction of the past ∼14-12 cal ka BP revealed that the colonization of eelgrass around 2,000 years ago transformed the rate and shape of the biogeochemical sink through ecosystem stabilization and sediment accretion, leading to higher concentrations of finer (mud) grain size particles, and an increase in accumulation rates of up to 24-fold for sedimentary OC and N. The onset of seagrass establishment was likely related to the change in coastal geomorphology during this time period, resulting in decreased water depth and hydrodynamic exposure of Gåsö, which created a favorable environment for seagrass growth.After the initial colonization period, the most drastic biogeochemical changes occurred between 1.1 and 0.7 cal ka BP.This time period, corresponding to the MCA, was characterized by warmer temperatures and lower storm activity (Hass, 1996), and these changes in climate could also have further stimulated seagrass productivity and meadow stability, which in turn enhanced the accumulation of organic-rich sediment even further.This underscores the significance of both abiotic and biotic factors in facilitating seagrass establishment and sustaining the essential ecosystem services provided by healthy seagrass meadows.(Eyring et al., 2016;O'Neill et al., 2016).
For information on GCMs included in the historical and scenario ensembles, see Table S1 in Supporting Information S1.

Environmental Conditions Conductive to Seagrass Colonization
The deglaciation of Northern Europe during the late Pleistocene and early Holocene resulted in an isostatic uplift (Rosentau et al., 2021;Stroeven et al., 2016) that led to coastal geomorphology changes at Gåsö.Sediment derived from the melting ice sheet and outflow of freshwater and associated clay-rich sediments from the Baltic basin through south-central Sweden deposited in coastal settings (Bergsten & Dennegård, 1988;Gyllencreutz, 2005;Gyllencreutz & Kissel, 2006), which was observed in the deepest layers of the Gåsö cores composed of brown clay of terrestrial lithogenic material rich in mud, as indicated by elemental composition, δ 13 C isotopic signatures below 20‰ and increased MS values.The freshwater outflow decreased around 10.7 cal ka BP (Gyllencreutz & Kissel, 2006) and the opening of the English Channel at about 8.5 cal ka BP (Conradsen & Heier-Nielsen, 1995) resulted in higher hydrodynamic activity and a larger inflow of marine waters from the Atlantic that carried coarser sediments into Skagerrak (Gyllencreutz, 2005).This coarsening of sediment was also clearly visible in the cores from Gåsö at that time, with more poorly sorted sediment reflecting a higher hydrodynamic activity.There were also two change points identified in Gåsö N at ∼9.7 and ∼8.6 cal ka BP, which might reflect this change in oceanic circulation with the inflow of Atlantic water.Gåsö N also had a more pronounced shift in MS values compared to Gåsö S, indicating a change from terrestrial to marine influence (López-Merino et al., 2017).
A general gradual increase in the marine influence from 14 to 12 cal ka BP until present was marked by shifts in both cores around 2.7-2.6 cal ka BP (Section 2) that resulted in enhanced sediment, OC and N accretion as well as higher carbonate content likely due to the isostatic uplift with the shore-level becoming similar to present conditions that created a more sheltered and stable environment.Seagrass colonization was detected around 2 cal ka BP based on shifts in δ 13 C, n-alkanes (with the appearance of the seagrass-related C 17 , C 19 , and C 21 nalkanes; Chevalier et al., 2015;Rosenbauer et al., 2006; which were absent prior to this time period), silt-clay content, elemental composition (including higher Br and Cl associated to the increased organic matter content) and sediment density, which shows alteration of the sediment characteristics similar to pedogenesis in terrestrial soils (Piñeiro-Juncal et al., 2020).The transformation of the sediment was likely triggered by the higher depositional rate linked to the reduction in hydrodynamic exposure.These levels of hydrodynamic exposure are prevalent within the range of contemporary Zostera spp.distribution (Dahl et al., 2020;Prentice et al., 2019;Short et al., 2007).The colonization period of seagrass also coincided with the RWP, which was characterized by regionally warm temperatures (Neukom et al., 2019).Zostera marina can tolerate a large range of temperatures with an upper limit of approximately 25-30°C (Nejrup & Pedersen, 2008), and while the optimum water temperature for plant growth is likely site-specific, it has been shown that productivity of Z. marina increases with temperature (within the upper thermal limit) (Rasmusson et al., 2019), which could have stimulated the spread of seagrasses.During the time of seagrass establishment, the oceanographic conditions were more dominated by Atlantic waters and less turbidity (Erbs-Hansen et al., 2012).After the colonization phase (Section 2), a shift in the CPM was identified for both cores due to the higher presence of marine-derived organic matter and more wellsorted sediments of finer grain size composition (Section 3).This change in the biogeochemical trends with enhanced OC and N accumulation occurred at 1.1-0.7 cal ka BP, which coincides with the MCA, a period of higher temperature (Mann et al., 2009) that may have favored seagrass productivity and growth even further.The last millennium constitutes a stable period suitable for seagrass growth, which is suggested by the more consistent hydrodynamic exposure and shore-level at Gåsö.

Ecosystem Services Related to Seagrass Colonization
The shifts seen in the sediment record following eelgrass establishment and stabilization phases consisted of enhanced sediment and mud accretion, and organic matter accumulation that could be related to both changes in the coastal geomorphology (a shallower and more sheltered environment) and through enhanced seagrass biomass and productivity across the region.Indeed, the effect of the seagrass canopy likely reduced water turbidity, resulting in an overall increase in ecosystem stability and productivity, and associated increases in the accumulation of OC and N. The establishment of the seagrass was a crucial requirement for the modification of sediment biogeochemical properties, with approximately 47%-48% of the organic matter originating from the seagrass biomass itself, while most of the remainder (43%-45%) was of macroalgae origin (Dahl et al., 2023).However, the carbon-(CAR: 6-8 g OC m 2 year 1 ) and nitrogen accumulation rates (NAR: 0.5-0.8g N m 2 year 1 ) at Gåsö for the latest ∼0.7 to 1 cal ka BP (Section 3) are still lower than accumulation rates in Zostera spp.meadows in general (Martins et al., 2022;Prentice et al., 2020) and in the lower range compared to levels found in other studies in the Skagerrak-Kattegat region, which estimated CAR between 6 and 134 g OC m 2 year 1 and NAR from 0.7 to 14 g N m 2 year 1 (Dahl et al., 2023;Leiva-Dueñas et al., 2023).However, these studies only assessed accumulation over shorter time periods (the last decades to century), while long-term accumulation rates tend to be lower for seagrass in general due to diagenesis and remineralization of the organic matter over time (Belshe et al., 2019).This indicates that while the organic matter has been gradually accumulating since the onset of seagrass colonization (Section 2), it has done so at a comparatively faster rate (20to 24-fold) since seagrass stabilization (Section 3 in comparison to Section 1).Although sediments tend to have a lower density at the surface due to condensation with depth, the clear decrease in DBD and increased sorting of finer sediment components following the onset of the organic-rich seagrass sediment build-up after the seagrass establishment shows the effect of the sediment stabilization by the seagrass structure.

Climate Change Impacts on Seagrass Plant Health Status
The predicted SLR by 2100 alone will likely not impact the present seagrass distribution at Gåsö.In fact, with an increase of about 0.5 m (for SSP5-8.5 in 2100), the seagrass meadow in the southern bay (Gåsö S) has the potential to colonize an extended area, primarily the inner part of the bay that is currently too shallow (less than 0.4 m) for Z. marina.The lateral expansion of the rhizomes from an established Z. marina meadow has been estimated at 16 cm year 1 (Olesen & Sand-Jensen, 1994), which would allow the seagrass meadow to adjust to the predicted SLR and likely thrive along the deeper part of the bay based on the current depth range of Z. marina on the Swedish Skagerrak coast.In contrast, the GIS-modeled SLR leads to a submerging of the middle part of the island, creating a channel connecting the southern and northern bays that will likely remobilize postglacial sand that could result in smothering of the seagrass meadow.This center part of the island surfaced around 800 to 700 years ago due to the isostatic uplift, which occurred during the time period of seagrass stabilization for Gåsö S (Section 3), and an opening of the channel in the future scenarios may change the hydrodynamics of Gåsö and increase the hydrodynamic activity.Though, SLR alone would not affect the hydrodynamic exposure with a calculated maximum increase of 3% in effective fetch by 2100 for the highest climate scenario (SSP5-8.5),increased sediment transport in a north-south direction due to the opening of the channel could lead to increased turbidity locally, which might negatively affect the contemporary distribution of seagrass at Gåsö.However, the potential impact from turbidity and sedimentation might be counteracted by increasing the available area for seagrass expansion through the deepening of the southern bay.
The projected regional SST increase until 2100 (with a predicted annual mean ranging from 11 to 14°C and mean summer maximum temperatures of 18-21°C in 2100 for the different climate scenarios) is within the growth range for Z. marina, which generally has a high tolerance for increased temperatures (Rasmusson et al., 2020).Historical periods of increased sea water temperature (i.e., the RWP and the MCA) were still lower compared to today (Lapointe et al., 2020) and therefore present temperatures are likely representing the warmest period during the colonization history of the seagrass meadows at Gåsö.Nejrup and Pedersen (2008) identified the optimum water temperature for Z. marina to range from 10 to 20°C, with an increase in mortality and a decrease in photosynthetic activity and growth occurring at 25-30°C.However, the SST data do not account for extreme events such as heat waves, which are expected to increase in the future.Such heat wave events could substantially increase the water temperature for short periods of time and result in major seagrass loss and negatively affect ecosystem services (Aoki et al., 2020(Aoki et al., , 2021)).Although an increase in water temperature, either as a short-term single event or as a gradual increase over time, might not strongly affect seagrass productivity and growth in the Skagerrak region, a higher temperature can lead to reduced water quality through stimulation of filamentousor microalgae growth (Moore et al., 2012), which likely has a stronger impact on the seagrass health status than the increased SST alone.
This decrease in modeled salinity (reaching around 25 for SSP3-7.0 and SSP5-8.5 in 2100) is also within the tolerance limit for Z. marina, which shows large plasticity in response to salinity as the species can grow in brackish environments down to around 5-6 (Boström et al., 2014).However, the salinity tolerance is also population-specific with Z. marina growing in the marine environment (as the seagrass at Gåsö) being less adapted to low salinity (Salo et al., 2014).The projected salinity decrease is still within the tolerance limit and will likely not affect the seagrass plant health to a large extent (Hellblom & Björk, 1999;Zhang et al., 2022).Reconstruction of salinity and sea water temperature in the Skagerrak region (Jiang et al., 1998;Ni et al., 2021) shows that these environmental parameters have likely remained within the tolerance limit for Z. marina growth throughout the Holocene.

Opportunities for Seagrass Conservation and Restoration
Although the environmental conditions thousands of years ago were in many ways different, reconstructing the interactions between environmental conditions and the colonization and establishment by seagrasses together with the prediction of potential climate change impacts provides important information for seagrass conservation and restoration, which could help guide management.This provides further evidence that the sediment characteristics and somewhat sheltered conditions (leading to reduced turbidity) are of importance for successful Z. marina replantation and re-colonization, as shown by several restoration projects (Orth et al., 2010).Positive feedback likely occurred following the establishment of seagrasses through enhanced water quality and irradiance by reducing resuspension and increasing the deposition of suspended particles, which resulted in enhanced clay and silt, and OC and N accumulation.Moksnes et al. (2021) proposed that suitable restoration sites in the Skagerrak-Kattegat area should have mud contents of less than 34% and OC levels below 5%.If the mud content and resuspension of sediment are high and the seagrass meadow is in an early stage of establishment, likely the seagrass plants cannot stabilize the sediment enough for sufficient irradiance to support plant growth and reproduction (Moksnes et al., 2016) and the Z. marina meadows at Gåsö established in coarse sediment with low mud (∼5-6%) and OC contents (0.5%-1.5%).This finding aligns with the ongoing replantation projects on the Swedish Skagerrak coast, where several tons of sand have been added to the bottom within a km-wide project area in order to stabilize the sediment prior to replantation of Z. marina shoots.
This study shows critical coastal geomorphology requirements for the success of seagrass establishment and colonization, including a low hydrodynamic exposure (with an effective fetch of <1.3 Lf [km]) and overall stability of the environment.This confirms the findings of previous studies showing that site selection for restoration projects needs to be carefully evaluated for a successful outcome of the restoration effort (Short et al., 2002;van Katwijk et al., 2009).As long-term monitoring of coastal habitats is commonly lacking, paleoreconstruction could help identify areas that are environmentally stable, which could be more suitable for habitat conservation, and to assess carbon and nitrogen abatement following seagrass establishment (i.e., estimated at between 63 and 68 Mg OC ha 1 ; Dahl et al., 2023;and 5 and 7 Mg N ha 1 over the past 2 cal ka BP) related to for example, Nationally Determined Contributions.The historical paleorecord in combination with regional climate models is useful in the assessment of risk related to climate change for cold-temperate Z. marina and provides insights into potential impacts linked to climate change scenarios (e.g., SLR, precipitation and storminess), such as the enhanced build-up of sedimentary organic matter during relatively warmer periods of increased temperature in the late Holocene.Zostera marina in Northern Europe inhabits the lower end of its thermal tolerance range and therefore, it has the potential to thrive under predicted warming scenarios (Rasmusson et al., 2019).However, the Z. marina meadows in Skagerrak are subjected to multiple threats (both from single events and through cumulative impacts; Baden et al., 2003;Moksnes et al., 2008Moksnes et al., , 2018) that might weaken their resilience against climate change threats (Björk et al., 2008;Unsworth et al., 2015).This multiplicity of potential disturbances needs to be considered to ensure successful habitat conservation.
Successful restoration and re-vegetation of seagrass meadows has been shown to increase ecosystem services derived from the build-up of sediment deposits (Greiner et al., 2013;Marbà et al., 2015) and enhance biodiversity, among other benefits for the people and the planet (Orth et al., 2020).However, many restoration projects have failed owing to the environmental conditions not being suitable for seagrass growth (van Katwijk et al., 2016).
This study provides further evidence of the importance of environmental conditions for natural seagrass colonization and this knowledge could thus help to identify sites suitable for successful seagrass conservation and restoration efforts.

Figure 1 .
Figure 1.Maps showing (a) the study locations (Gåsö South [S] and Gåsö North [N]) in the vicinity of the Gullmar Fjord, (b) substrate types at Gåsö, and (c) water depth ranging from 0 to 6 m and current distribution of eelgrass.

Figure 2 .
Figure 2. Sediment stratigraphy of the cores.The shifts between the distinct layers were identified by visual inspection and from the digital imagery obtained from the XRF-core scanning.The dashed lines indicate the calibrated 14 C-dates (see TableS1in Supporting Information S1).

Figure 3 .
Figure 3. Summary of the physical and biogeochemical properties of the sediment, including organic carbon (OC), nitrogen (N), stable carbon isotopes (δ 13 C), proportion of carbonate content (CaCO 3 ), dry bulk density (DBD), proportion of silt-clay (mud) content, degree of sorting (sorting coefficient), magnetic susceptibility (MS), average chain length (ACL) of n-alkanes (for homologs between C 17 and C 35 ) and seagrass biomass associated n-alkanes (i.e., C 17 , C 19 , and C 21 ) (a) along core depths and ages (shown as calibrated thousand years before present), and (b) zoomed in on the transition period of seagrass colonization (upper ∼80 cm).Note that panel (b) shows sediment depth on the primary axis for better representation of the last ∼2 ka cal BP.The black line in panels (a) and (b) indicates the first appearance of seagrass-derived sediment in the paleorecord, as evident by the introduction of Zostera marina biomass associated n-alkanes (Figure S3 in Supporting Information S1).

Figure 4 .
Figure 4. Changes in sediment characteristics shown as biplot of PC1 and PC2 scores in relation to sediment depth indicated by colors (a), and arc-length distance with age (calibrated thousand years before present) (b).The horizontal orange (Gåsö S) and green (Gåsö N) lines in panel (b) show the change point probability for the arc-length distance of PC1 and PC2 for each core.The top panel in (a) shows the different shapes of the PC curves for the two cores, which is due to differences in sediment composition as seen in the PCA biplot (lower plot in a).The PCA biplot presents the relationship of sediment characteristics for PC1 and PC2, with colors indicating core depth.The arc-length distance in panel (b) represents the distance of PC1 and PC2 to the modelled curve, and was applied due to correlation of the two PC-axes.The black line in panel (b) indicates the first appearance of seagrass-derived sediment in the paleorecord.Section 1 shows the time frame prior to seagrass establishment (as a baseline), whereas Section 2 presents the period around seagrass colonization and Section 3 presents the seagrass stabilization phase (see Table2).

Figure 5 .
Figure 5. Changes in effective fetch (km) (a), and land uplift for 1, 2, 4, and 11 cal ka BP (calibrated thousand years before present) (b).Modeled data on shore displacement were downloaded from the Swedish Geological Survey and show the isostatic uplift and shore-level changes at Gåsö.The black line in panel (a) shows the first appearance of seagrass-derived sediment in the paleorecord.

Figure 7 .
Figure 7. Predicted annual mean sea surface temperature (SST) (a), salinity (b), and surface runoff (c) changes until 2100 for four Shared Socioeconomic Pathway (SSP) scenarios.Data are ensemble means from CMIP6 global climate models (GCMs), retrieved for the closest grid point to Gåsö(Eyring et al., 2016;O'Neill et al., 2016).For information on GCMs included in the historical and scenario ensembles, see TableS1in Supporting Information S1.

Table 1
Factor Loadings of the Principal Component Analyses (PCA) Based on the Biogeochemical Data From Both Sediment Cores (Gåsö S and N) Combined

Table 2
Sediment Layer Thickness, Age, Sediment Accretion Rates and Organic Carbon and N Accumulation Rates for Three Sediment Sections Identified Based on Change Point Modeling Analyses of Arc-Length Distance for PC1 and PC2 in Gåsö S and N Cores Note.Section 1 shows the period before establishment of the seagrass (as a baseline), whereas Section 2 shows the time around seagrass colonization and Section 3 presents the period for seagrass stabilization.Cal ka BP, calibrated thousand years before present; SAR, sediment accumulation rate; MAR, mass accumulation rate; CAR, organic carbon accumulation rate; NAR, nitrogen accumulation rate.All accumulation rates are presented as mean ± SD.DAHL ET AL.