Introduction

Seagrasses occupy only 0.1% of the ocean surface but are considered one of the largest carbon sinks worldwide1,2,3,4. Unlike terrestrial ecosystems, which store organic carbon (Corg) mainly in living biomass, Corg stocks in seagrass meadows are mainly found in their soils, where it can accumulate over millennia5,6. Given their ability to capture and retain Corg and elevate the seafloor, seagrass meadows play a significant role in climate change mitigation and adaptation by sequestering greenhouse gases and buffering the impacts of rising sea levels4.

Seagrasses comprise more than 70 species of seagrasses7 that have a wide variation in traits including differences in primary production rates, below ground biomass, the recalcitrance of the Corg in their organs, and the ability of their canopies to trap allochthonous carbon8. Seagrasses occur across different depositional environments9, which results in highly variable soil Corg stocks among seagrass habitats (up to 18-fold10). Such differences are the result of the interactions among multiple biotic and abiotic factors (e.g. species composition, geomorphology, hydrodynamic energy and water depth) acting in the water column, canopy and the soils, that can affect rates of Corg deposition and accumulation over millenary time scales11,12,13,14,15. However, global estimates of Corg stocks in seagrass ecosystems have been calculated based on a very narrow data set, based on few species and habitats mainly from sites within the Mediterranean, Northern Atlantic, and eastern Indian Oceans6,16. Indeed, the paucity of Corg burial estimates in seagrass ecosystems is limiting our understanding of their carbon sequestration capacity and the assessment of differences in Corg storage across seagrass habitats. In order to improve global estimates of Corg stocks in seagrass meadows, research is needed that expands the current databases to include assessments in regions representing geomorphological and biological characteristics thus far underrepresented in the literature, such as arid regions.

The Red Sea, a tropical, arid region lacking riverine inputs, provides a case for a region currently underrepresented in the assessment of seagrass Corg stocks. Although Corg stocks within terrestrial ecosystems of the Red Sea arid bioregion are relatively low, coastal environments may serve as critical and important carbon sinks within this region17. While seagrass communities in the Red Sea are mainly dominated by opportunistic and small species18, even these small seagrasses have been documented to significantly enhance sediment stabilization and accumulation19 and potentially contributing to Corg storage10,20. Given that there are over 370 km2 of seagrasses in Saudi Arabia18,21, there is potential for these meadows to contain substantial stores of Corg. This study aims to provide estimates of soil Corg stocks and accumulation rates in seagrass meadows from the Red Sea, thus contributing to achieve a more balanced representation of variation among seagrass habitats by considering underrepresented regions10, and to place our results within a larger context by comparing these data to estimates from global datasets.

We combine estimates of soil Corg density down to 1 m depth with soil chronologies derived from 14C age dating to estimate (a) Corg stocks within the top meter of soil, and (b) the accumulation rate of Corg over the last millennia. We also estimate the contribution of autochthonous and allochthonous sources to the seagrass soil Corg pool and determine soil grain-size composition to examine their relationships with soil Corg pools.

Results

Soils in Red Sea seagrass meadows were mainly constituted of clay and silt particles (37 ± 0.7% on average), with a relatively high abundance of very fine sands (21 ± 0.4%) compared to fine sands (16 ± 0.4%), medium sands (12 ± 0.3%) and coarse sands (14 ± 0.7%) (Table 1). The average (±SE) dry bulk density in 1 m-thick Red Sea seagrass soils was 1.05 ± 0.01 g cm−3 and the average Corg concentration was 0.33 ± 0.01% and 3.05 ± 0.07 mg cm−3. The soil accretion rates in Red Sea seagrass meadows over centennial time scales ranged broadly, between 0.02 and 1.57 cm yr−1 (Table 2). In 1 m-thick soil deposits and over 500 to 2,000 years of accumulation, seagrass meadows accumulated, on average, 3.4 ± 0.3 kg Corg m−2 at a rate of 6.8 ± 1.7 g Corg m−2 yr−1 (Table 2).

Table 1 Average ± standard error (SE) of (a) dry bulk density (in g cm−3), organic carbon (Corg) content (in % and mg cm−3), δ13C signatures of soil organic matter and (b) sediment grain-size content at Red Sea seagrass soil cores (for the total length of core sampled; see Supporting Information Table B for further details).
Table 2 Soil accretion rates (cm yr−1) and organic carbon (Corg) stocks (kg m−2) and accumulation rates (g m−2 yr−1) in Central Red Sea seagrass meadows (in 1 m-thick soils).

The δ13C values of sedimentary organic matter in seagrass soils averaged −16 ± 0.2‰ (Table 1). The δ13C values of potential organic sources into seagrass soils collected at the four study sites are presented in Table 3. The mixing models applied indicated that seagrass detritus was the most important source of soil Corg in Red Sea seagrass meadows (41%), while mangrove plus halophytes and seaweed plus seston inputs were less important (32 and 27%, respectively; Fig. 1 and Supporting Information Table A).

Table 3 Mean (±SE) of isotopic carbon values (‰) of potential organic sources into seagrass soils collected at the four study sites.
Figure 1
figure 1

Boxplot showing the results of the isotopic mixing model (IsoSource software). Proportion (25%, 50% and 75% quantiles) of autochthonous (i.e. seagrass matter) and allochthonous (i.e. mangroves plus halophytes, seaweed and seston) Corg to the seagrass soil Corg pool (top 60 cm of the cores) based on study site (Economic City, Khor Alkharar, Petro Rabigh and Thuwal Island) and seagrass species (Enhalus acoroides, Halodule uninervis, Halophila stipulacea, Thalassia hemprichii and Thalassodendrum ciliatum).

Soil biogeochemical characteristics (dry bulk density, Corg content (in % and mg cm−3), δ13C of soil organic matter, and sediment grain-size composition) differed significantly among study sites and meadows with distinct species composition (P < 0.05; Table 4). Soil accretion rates (cm yr−1) and Corg accumulation rates (g m−2 yr−1) also differed significantly among locations (P < 0.01; Fig. 2A) but did not differ among the seagrass species tested (P > 0.05; Fig. 2B), while Corg stocks (kg m−2 in 1 m-thick soils) did not differ either among locations or meadows with distinctive species composition (P > 0.05; Table 4; Fig. 2A,B).

Table 4 Results of General Linear Models. Soil dry bulk density, soil accretion rates, organic carbon (Corg) concentration, stable carbon signatures (δ13C) of sedimentary organic matter, and sediment grain size fractions in response to Site (Thuwal Island, Economic City, Petro Rabigh and Khor Alkharar) and Seagrass species (Halophila stipulacea, Thalassia hemprichii, Enhalus acoroides, Thalassodendrum ciliatum and Halodule uninervis).
Figure 2
figure 2

Soil accretion rates (cm yr−1) and organic carbon (Corg) content (%), stocks (kg m−2 in 1 m-thick soils; extrapolated when necessary, see methods) and accumulation rates (g m−2 yr−1), stable carbon isotope signatures of soil organic matter (δ13C; in ‰) and sediment grain-size fractions in Central Red Sea seagrass meadows based on study site (Economic City, Khor Alkharar, Petro Rabigh and Thuwal Island) and seagrass species (Enhalus acoroides, Halodule uninervis, Halophila stipulacea, Thalassia hemprichii and Thalassodendrum ciliatum). Mean ± SE values are reported. The results of Tukey HSD posthoc tests to assess pairwise differences are indicated: different letters (a,b,c,d) indicate significant differences (P < 0.05).

Soil dry bulk density in meadows at Economic City (0.9 ± 0.01 g cm−3) was lower compared to those at Petro Rabigh and Khor Alkharar (1.1 ± 0.01 g cm−3 in both cases), while seagrass soils at Thuwal Island had the highest soil dry bulk density values (1.3 ± 0.01 g cm−3). The soil Corg content was significantly higher at Economic City (0.4 ± 0.01% Corg) and Khor Alkharar (0.3 ± 0.01% Corg) compared to meadows at Thuwal Island (0.2 ± 0.01% Corg) and Petro Rabigh (0.3 ± 0.03% Corg; P < 0.001; Fig. 2A). At Petro Rabigh, seagrass soil accretion rates (1.3 ± 0.3 cm yr−1) and Corg accumulation rates (27.3 ± 1.5 g Corg m−2 yr−1) were significantly higher (P < 0.001) than those from the other sites (averaging 0.14 ± 0.1 cm yr−1 and 4.7 ± 1.0 g Corg m−2 yr−1; Fig. 2A). Soil Corg stocks in 1 m-thick deposits were not significantly different among locations (ranging from 2.7 ± 0.15 to 4.0 ± 1.9 kg Corg m−2).

Seagrass soils at Thuwal Island had higher amounts of fine soil particles (78% of clay and silt and very fine sands) compared to the other locations (ranging from 51 to 63%; Fig. 2A). Meadows at Economic City had a relatively larger proportion of coarse sands (25%) compared to the other locations (ranging from 1 to 14%). Meadows at Thuwal Island and Petro Rabigh were 13C-depleted (averaging −17.8 and −22.7‰, respectively) compared to the other locations (δ13C ranging from −12.7 to −15.4‰; Fig. 2A). The mixing models applied indicated that seagrass detritus was a relevant source of soil Corg in meadows at Economic City (51%) and Khor Alkharar (45%), but a relatively minor contributor at the rest of locations (ranging from 12 to 31%; Fig. 1A). Mangrove plus halophyte constituted the main soil Corg sources in Petro Rabigh meadows (73%) and Thuwal Island (38%), while the contribution of seston ranged between 15 and 32% among the study sites.

Dry bulk density was significantly lower in T. hemprichii meadows (0.7 ± 0.02 g cm−3) compared to the other seagrass species (ranging from 0.9 to 1.3 g cm−3). The soil Corg content was significantly higher in T. hemprichii (0.7 ± 0.02% Corg) and E. acoroides (0.4 ± 0.01% Corg) compared to H. stipulacea, T. ciliatum and H. uninervis meadows (ranging from 0.2 ± 0.01 to 0.3 ± 0.01% Corg; Fig. 2B). The soil Corg concentration (% and mg Corg cm−3) decreased down core. Soil Corg stocks and accumulation rates were not significantly different among meadows with distinct species compositions (ranging from 2.8 ± 0.4 to 4.8 ± 0.3 kg Corg m−2 and 5.5 ± 1.3 to 12.47 ± 6.1 g Corg m−2 yr−1, respectively; P > 0.05). Seagrass meadows formed by H. stipulacea, T. hemprichii and E. acoroides contained higher amounts of clay and silt particles (ranging from 42 to 48%) compared to those of T. ciliatum and H. uninervis (ranging from 21 to 22%; Fig. 2B).

Meadows formed by T. hemprichii, E. acoroides and H. uninervis were 13C-enriched (ranging from −11.2 to −14‰) compared to H. stipulacea (−16.0‰) and T. ciliatum (−20.0‰) meadows (Fig. 2B). The trends in δ13C values with soil depth remained stable in T. hemprichii, E. acoroides and H. uninervis meadows (Supporting Information Fig. A), indicating that either the organic matter decomposition in the soils or the inputs of organic matter remained stable. However, the δ13C values become more negative with depth/ageing in H, stipulacea and T. ciliatum meadows, in particular below ~60 cm soil depth. The 13C-depletion in soil organic matter could either indicate the lack of seagrass matter inputs below cm ~60 in the cores or the decomposition of carbohydrates and proteins with ageing, which have more positive δ13C values than the remaining organic matter (e.g. lignin22). In order to constrain the uncertainties mentioned above, we run the mixing models using average δ13C values within the top 60 cm of the cores only. Despite the δ13C values only remained stable within the top 60 cm in E. acoroides (R2 = 0.05) and H. stipulacea (R2 = 0.14) meadows (i.e. the δ13C values significantly increased with soil depth in T. hemprichii meadows (R2 = 0.60) and significantly decreased with soil depth in H. univervis (R2 = 0.30) and T. ciliatum (R2 = 0.34)), their variability is relatively small within the top 60 cm of soil compared to 1 m-thick soils, thereby constraining the uncertainties associated with diagenetic effects. The mixing models applied indicated that seagrass detritus was a significant source of soil Corg in all meadows studied: T. hemprichii (58%), E. acoroides (59%), H. uninervis meadows (45%) and H. stipulacea (40%), except T. ciliatum (21%; Fig. 1B). Mangrove plus halophyte constituted the main soil Corg sources in T. ciliatum meadows (54%), while the contribution of seston ranged from 24 to 29% in the meadows among meadows with different species composition.

Discussion

Seagrass meadows represent one of the most important vegetated communities in the otherwise arid and oligotrophic Red Sea region. Seagrasses are widely distributed along the Kingdom of Saudi Arabia Red Sea coast18, and the soils underneath those seagrass meadows contain considerable Corg stocks. The soil Corg content of Red Sea seagrass in 1 m-thick soils (3.4 kg Corg m−2 on average) is well below the values from global estimates (ranging from 12 to 83 kg Corg m−26). The relatively low Corg sink capacity of Red Sea seagrasses could be due to the extreme environmental conditions such as nutrient limitation and high temperature, reducing the growth rates of the seagrasses and increasing the rate of respiration in the soil21,23.

This disconnection between Red Sea seagrass Corg stocks and the global estimates are likely linked to the very limited data set used to produce global estimates, which was biased by the extremely high Corg content of soils from Mediterranean Posidonia oceanica meadows6. Recent estimates of soil Corg stocks in low biomass seagrass species from Abu Dhabi (H. uninervis, Halophila ovalis and H. stipulacea; ranging from 0.2 to 10.9 kg Corg m−224), Asia (Zostera spp. and T. hemprichii; ranging from 3.8 to 8.9 kg Corg m−220), and Australia (H. uninervis and T. hemprichii; ranging from 2.3 to 5.0 kg Corg m−210) are within the range of Corg stocks estimated for Red Sea seagrasses. Hence, it is likely that Corg stocks in Red Sea seagrass soils tend to be in the lower range but are not necessarily below those in seagrass soils in all other meadows, suggesting a need to update the global estimate of Corg content in seagrass soils using a more balanced geographical distribution of seagrass meadows, but also accounting for habitat variability (i.e. diversity of morphological traits across species and geomorphology).

Moreover, long-term Corg accumulation rates in Red Sea seagrass (7 g Corg m−2 yr−1 on average) are lower than previous estimates for large and long-living Posidonia spp. meadows in the Mediterranean (84 ± 20 g Corg m−2 yr−1) and Australia (12 ± 7 g Corg m−2 yr−1)25, but similar to estimates from low biomass and fast-growing seagrass species (i.e. Zostera spp. and T. hemprichii) from Japan (ranging from 1.8 to 10.1 g Corg m−2 yr−120).

Nevertheless, owing to the low Corg storage in desert land areas and the relatively large seagrass habitat in the Red Sea coast of Saudi Arabia (370 km218), seagrasses constitute hotspots of Corg storage in this extremely arid region. Multiplying this area by the average Corg stocks in 1 m-thick soils and the average Corg accumulation rates, yields a total estimate of 1.2 ± 0.1 Tg Corg at 1,700 ± 337 Mg Corg yr−1 in Saudi Arabia meadows. The total Corg stored in Saudi Arabia’s seagrass meadows is roughly equivalent to 7 years of total CO2 emissions from fossil-fuel burning, cement production, and gas flaring by Saudi Arabia, while seagrass meadows sequester annually around 0.8% of these emissions (Saudi Arabia emissions estimated at 0.16 Tg C at 2014 rates26).

Long-term (i.e. based on 14C) soil accumulation rates in Red Sea seagrass (ranging from 0.2 to 16 mm yr−1; 2.4 mm yr−1 on average) are within the range of 14C-derived values reported in previous studies from Australia (ranging from 0.15 to 2.5 mm yr−125,27), Japan (ranging from 0.37 to 1.3 mm yr−120), Spain and Italy (ranging from 0.6 to 4.9 mm yr−15,25,28,29). The capacity of seagrass to elevate the seabed through sediment accretion has been previously recognized as a major component of their role in climate change adaptation4, as it helps mitigate against sea level rise. The results obtained in this study confirm that seagrass meadows in the Red Sea play a significant role in climate change adaptation through the protection against sea level rise, despite this being an arid region with very limited supply of terrestrial sediment via run-off. The sea level rise in the coast of Saudi Arabia has been estimated at 2.2 ± 0.5 mm yr−1 30, hence seagrass ecosystems along the Saudi Arabian coast have been playing a key role in offsetting relative sea level rise.

The large variability in Corg concentrations, stocks and accumulation rates among seagrass habitats (i.e. species composition and location) support the hypothesis that Corg storage in seagrass soils is influenced by interactions of biological, chemical and physical factors within the meadow14,15,31. Despite that no significant differences in Corg stocks among locations existed (at 95% confidence), the biogeochemical characteristics of the cores allowed the reconstruction of the processes and drivers involved in Corg storage (Fig. 3A). Soil Corg was negatively correlated to soil dry bulk density (R2 = 0.34), as previously shown in a range of sediments including seagrass soils32,33, which could explain the significant differences found in Corg concentration (%) and the lack of differences in Corg stocks (kg m−2).

Figure 3
figure 3

Biplots showing the relationships among the variables studies in the seagrass cores from the Red Sea based on study site and seagrass species. <0.125 mm (%) indicates the percentage of clay and silt and very fine sands within the bulk soil.

The relatively high soil Corg concentrations (%) at Economic City could be related to the relatively high accumulation of seagrass detritus and abundance of fine sediments. These results support the hypothesis that the seagrass plants themselves play a key role in determining the amount of Corg available for burial14, while the presence of fine sediments tends to reduce remineralization rates due to lower oxygen exchange and redox potentials11,34,35. The mechanisms behind Corg accumulation and preservation in seagrass meadows at Petro Rabigh appear to be mainly related to the relative high soil accumulation rates together with large fluxes of Corg from adjacent mangrove and tidal marsh ecosystems. Previous studies have shown that high soil accumulation rates in seagrass meadows, linked to the capacity of their canopy to tap and retain sediment particles12,36, the hydrodynamic energy and the production of biogenic carbonates within the meadow37,38, contribute to higher accumulation and preservation of Corg after burial14. Petro Rabigh is an enclosed environment surrounded by mangrove forests, which has been shown to largely contribute to soil Corg storage in adjacent seagrass meadows15,39. The relatively low soil Corg storage of seagrass meadows at Thuwal Island could be explained by the relatively low contribution of seagrass detritus to the soil Corg pool and the low soil accumulation rates (Fig. 3A).

Clear differences were observed among meadows with distinct species composition, with the highest soil Corg concentrations (%) found in meadows composed of the largest seagrass species T. hemprichii and E. acoroides (Fig. 3B). However, the relatively low soil dry bulk density found in these meadows led to similar Corg stocks among all meadows studied. The results obtained in this study show that soil Corg concentration was influenced by the relative contribution of seagrass detritus to the soil Corg pool and the amount of fine sediments, which support the results obtained in previous studies14,15,25. The relatively high soil Corg concentration and seagrass contribution to the soil Corg pools in T. hemprichii and E. acoroides could be explained by the highest above- and below ground biomass of stands formed by these species (ranging from 72 to 87 g DW m−2 and 210 to 392 g DW m−2, respectively) compared to the other seagrass species studied (ranging from 2.3 to 27 g DW m−2 and 2.6 to 61 g DW m−240. This study supports previous research reporting that the intrinsic properties of the seagrass themselves (e.g. canopy structure, below- and above-ground biomass, and productivity) can influence soil Corg storage14,15,31. Moreover, the relative constant C stable isotope signatures along the cores confirm the stability of organic sources to soil Corg pools, except for H. stipulacea and T. ciliatum meadows (i.e. δ13C values decreased below cm 60), which may indicate that seagrass meadows have only been present for the last centuries at these locations. The presence of coarse soil fibers throughout ~14 cores indicated that seagrasses were present at the coring sites throughout the period reconstructed or the soil depth studied. However, in half of the cores seagrass fibers disappeared at 25–60 cm depth, which could be either due to seagrass absence or the decomposition of coarse organic matter with ageing. Indeed, with the proxies analyzed here it is not possible to assure that the seagrass species occurring at present have remained the same through time.

Our results contribute to gaps in the existing global database on seagrass meadow Corg stocks and accumulation rates, which were thus far lacking information from seagrass species in arid environments and suggest that even meadows comprised of ephemeral seagrass species can play an important role in Corg sequestration.

Material and Methods

Study site and sampling

This study was conducted in four locations (Thuwal Island, Petro Rabigh, Economic City and Khor Alkharar) along 80 km of the Kingdom of Saudi Arabia coastline in the Central Red Sea (Fig. 4). Seagrass meadows are found along the Saudi coast, mainly composed by H. stipulacea, T. hemprichii, E. acoroides, T. ciliatum and H. uninervis18.

Figure 4
figure 4

Location of seagrass meadows sampled in Saudi Arabia, Central Red Sea. The map was produced with ArcMap Version 10.2. Background map credits: the World Administrative Divisions layer provided by Esri Data and Maps, and DeLorme Publishing Company. Redistribution rights are granted http://www.esri.com/~/media/Files/Pdfs/legal/pdfs/redist_rights_103.pdf?la=en. The seagrass species present within the meadows surveyed at each study site are indicated.

Seagrass meadows at Thuwal Island grow on shallow soil of weathered coral and are located near the fisherman city of Thuwal41. Petro Rabigh is a major industrial and petrochemical complex, whereas the Economic City is a newly developed city and harbor complex subject to intense coastal development41,42. The Khor Alkharar lagoon encompasses a relatively undeveloped coastal plain and is permanently connected to the Red Sea.

Twenty-seven soil cores were sampled in 1 to 7 m-deep mono-specific seagrass meadows using manual percussion and rotation (PVC pipe with an inner diameter of 60 mm; Supporting Information Table B). Three to four replicate cores were sampled within 100 m2 of each mono-specific seagrass meadow at each site (three cores at Thuwal Island, 10 cores at Economic City, four cores at Petro Rabigh and 10 cores at Khor Alkharar). One third of the cores collected at each site were kept inside the PVC corers and transported to the laboratory (hereafter referred to as ‘whole cores’). The other cores from each study site were sampled in the field using a corer consisting of a PVC pipe with pre-drilled holes in the sidewall (3 cm wide and 3 cm apart; hereafter referred to as ‘port cores’), allowing sub-sampling of soil samples along the core in the field by inserting 60 ml syringes into the pre-drilled holes along the PVC pipes (Supporting Information Table B).

The total length of the core barrel used, the empty space inside the barrel before retrieval, the length of barrel outside the soil before retrieval, and the length of retrieved seagrass soil were recorded in order to correct the core lengths for compression effects and all variables studied here are referenced to the corrected, uncompressed depths (Supporting Information Table B). The volume of each subsample retrieved from the port cores was recorded in the field. The whole cores were sealed at both ends, transported vertically and stored at 4 °C before processing in the laboratory.

Plants, seaweeds and seston (>0.7 μm) were sampled across the study sites for isotopic characterization of the potential sources of organic matter in seagrass soils. Seagrasses included E. acoroides, T. hemprichii, T. ciliatum, H. stipulacea and H. uninervis. Mangroves and halophytes included Avicennia marina, Salicornia spp., Zygophyllum cocenium, Anabasis setifera and Suaeda monoica. Seaweeds included Padina spp., Colpomenia sinuosa, Sargassum spp. and Turbinaria ornata; seston >0.7 μm. The full dataset and details on methods can be found in Almahasheer et al.17.

Laboratory procedures

The whole cores were opened lengthwise and cut into 1 cm-thick slices, and each slice together with the sub-samples from the port cores were oven-dried at 60 °C until constant weight to determine the dry bulk density (g cm−3). All samples from the port cores and every second slice of the whole cores were then grounded in an agate mortar and subdivided for analysis.

For the analyses of soil organic carbon (Corg) and stable isotope composition (δ13C), 1 g of ground sample was acidified with 4% HCl to remove inorganic carbon, centrifuged (3400 revolutions per minute, for 5 min), and the supernatant with acid residues was carefully removed by pipette, avoiding resuspension. The sample was then washed with Milli-Q water, centrifuged and the supernatant removed. The residual samples were re-dried and then encapsulated for C analyses using a Thermo Delta V Conflo III coupled to a Costech 4010 at the UH Hilo Analytical Laboratory, USA. The content of Corg was calculated for the bulk (pre-acidified) samples. Carbon isotope ratios are expressed as δ values in parts per thousand (‰) relative to the Vienna PeeDee Belemnite standard. Replicate assays and standards indicated measurement errors of 0.01% for Corg content and 0.06‰ for δ13C.

For sediment grain-size analyses, a Mastersizer 2000-Malvern laser-diffraction particle analyzer was used following sieving (1 mm) and digestion of <1 mm soil samples with 30% hydrogen peroxide. Grain size fractions were categorized following Wentworth scale: clay and silt particles (<0.063 mm), very fine sand (>0.063 <0.125 mm), fine sand (>0.125 <0.25 mm), medium sand (>0.25 <0.5 mm), and coarse sand (>0.5 <0.75 mm)43.

A total of 58 radiocarbon analyses were conducted in the 27 cores sampled (1–5 analyses per core) at the AMS Direct Laboratory (USA) following standard procedures44. Samples consisted of pooled shells and bulk soil (Supporting Information Table C). Shells were partially digested with 10% HCl, rinsed in ultrapure MQ water in order to remove fine sediment particles, inspected under a stereomicroscope for absence of attached reworked materials, and dried at 60 °C to a constant weight before radiocarbon dating. The 14C age-depth models were produced using the R routine “Bacon” for Bayesian chronology building45, after 14C calibration using the marine13 radiocarbon age calibration curve46 taking into account a local Delta R of 110 ± 38 years47. From the Bacon routine output, the mean age was used to produce an age-depth weighted regression model forced through 0 (0 cm is cal. BP: 1950), using as weight the sum of the Euclidean distance of the minimum and maximum ages. In four cores, the 14C results indicated either that the samples dated were modern (younger than ~400 years) or that the core was mixed (Supporting Information Table C), and therefore we did not produce age-depth models for these four cores. The relatively unknown marine reservoir effects at our study sites (and changes through time) remains a big assumption when calibrating 14C ages48. All 14C results used to model core age-depth chronologies in this study are older than ~400 years, and therefore, the burning of fossil fuels did not affect our 14C-derived soil accumulation rates. All dates reported in this paper are expressed as radiocarbon calibrated years.

Numerical procedures

Corg density (g Corg cm−3) was calculated for each soil depth in each core by multiplying the sediment dry bulk density (g cm−3) by the Corg concentration (%). For soil depths where Corg content (%) was not analyzed, we extrapolated the %Corg (i.e. by averaging the %Corg between above and below depths) and multiplied the %Corg by the dry bulk density (g cm−3) to obtain Corg density (g Corg cm−3). To allow direct comparisons among locations, the soil Corg standing stocks per unit area (cumulative stocks; kg Corg m−2) were standardized to 1 m-thick deposits. The total soil depth sampled was higher than 100 cm in 13 cores out of 27 cores sampled and therefore, no extrapolation was required for these cores. However, the soil Corg stocks in 1 m-thick soil deposits were inferred in 14 cores (soil depths sampled ranged from 44 to 64 cm) to 1 m, by extrapolating linearly integrated values of Corg content (cumulative Corg stock; kg Corg cm−2) with depth. Correlation between extrapolated Corg stocks from 44 cm to 1 m and measured Corg stocks in 1 m soil cores was r = 0.80 (P < 0.001; Supporting Information Fig. B). Note that scaling Corg stocks to 1 m using this method could either lead to over- or underestimates of Corg stocks.

Soil accretion rates (expressed in cm yr−1), soil accumulation rates (expressed in g DW m−2 y−1) and soil Corg accumulation rates (expressed in g Corg m−2 y−1) for the last millennia were estimated using 14C age-depth models (Table 2). Accumulation rates of Corg were calculated in 24 out of the 27 cores sampled by multiplying the Corg inventories in 1 m-thick soil by the average 14C soil accretion rate.

Analyses to test for differences in the variables studied among sites were performed using General Linear Model procedures in SPSS v. 14.0. General Linear Models were used to test for differences in dry bulk density (g cm−3), soil accretion rates (cm yr−1), soil Corg concentration (in %), soil Corg density (mg cm−3), soil Corg stocks (kg m−2 in 1 m-thick soils) and soil Corg accumulation rates (g m−2 yr−1), δ13C signatures of organic matter, and sediment grain size fractions among study sites and species composition (Table 4), followed by Tukey HSD posthoc tests to assess pairwise differences (Fig. 3). All response variables were square-root transformed prior to analyses and had homogenous variances. Study site (Thuwal Island, Petro Rabigh, Economic City and Khor Alkharar) and seagrass species (H. stipulacea, T. hemprichii, E. acoroides, T. ciliatum and H. uninervis) were treated as fixed factors in all statistical models (probability distribution: normal; link function: identity).

Stable Isotope Mixing Models were used to estimate the proportion of the autochthonous and allochthonous Corg to the seagrass soil Corg pool using δ13C and a one-isotope three-source mixing model49,50. The δ13C signatures within the top 60 cm of each core were pooled and analysed for the probability of relative organic matter contribution to soil stocks using Stable Isotope Mixing Models in R (‘simmr’ and ‘rjags’ packages)51. The δ13C signatures of potential Corg sources (seagrass was considered as autochthonous Corg, while mangroves plus halophytes, and seaweed plus seston were considered allochthonous Corg) in the four study sites were obtained from Almahasheer et al.17. The ‘simmr_mcmc()’ function works by repeatedly producing potential values of the proportional contribution of source material through a Markov chain Monte Carlo, with initial burn-in iterations (1,000) discarded and subsequent iterations (10,000) stored for use in the posterior distribution and analyses of the data51,52. Model convergence was confirmed using diagnostic plots and upper confidence intervals, while no overlap between source δ13C signature means ± standard deviations were observed. The simmr package allows for the incorporation of δ13C uncertainty into mixing models, while producing a Bayesian quantification of the most likely source contributors where there is a greater than n +1 sources when matching against n isotopes52. The dataset generated for this manuscript is provided as Supporting Information.