Introduction

The configuration of landscape mosaics strongly influences the strength and spatial patterning of ecological processes (Wiens 1995), which may be further affected in landscapes modified by human activities (Fischer and Lindenmayer 2007). Land-use changes are fundamentally restructuring, altering and redeveloping our biosphere while providing human benefits (Foley et al. 2005). In the highly productive coastal environment, habitat composition and connectivity across land–sea gradients control principal processes that in turn are subject to global environmental change (Grober-Dunsmore et al. 2009; Olesen et al. 2018). Today, human dependence on coastal resources is constantly growing as population density increases (Crossland et al. 2005; Pittman et al. 2019), inducing significant pressure from anthropogenic activities within the land–sea interface (Sloan et al. 2007). To accurately address the strength and spatial patterning of ecological processes in the complex coastal environment, it is hence particularly important to understand the joint effect of landscape configuration, land-use change and ecosystem degradation.

The coastal environment encompasses many vegetation-covered habitats that due to their high productivity and efficiency in trapping and storing organic material play an important role as natural carbon sinks, thereby safeguarding critical climate regulation services (Mcleod et al. 2011). Mangrove forests and seagrass meadows are key habitats widely distributed in tropical and subtropical coastal seascapes (Short et al. 2007; Hamilton and Casey 2016) and considered among the most important, so called, “blue carbon” habitats that mitigate climate change through a high sequestration of atmospheric CO2 (Mcleod et al. 2011; Serrano et al. 2019). The creation process for long-term storage of carbon (Smith 1981; Mateo et al. 1997; Gonneea et al. 2004) and the levels of organic matter accumulation in mangrove- and seagrass-sediments depend on factors operating on different scales, including primary production (Kristensen et al. 2008; Duarte et al. 2010), contemporary environmental conditions and benthic processes (Marbà et al. 2006; Kristensen et al. 2008), habitat configuration and composition (e.g. size, patchiness and edge perimeter; Ricart et al. 2015a, 2017; Oreska et al. 2017), allochthonous input from nearby environments, such as rivers and adjacent coastal habitats (Kennedy et al. 2010; Saintilan et al. 2013; Ricart et al. 2020), and removal of particulate and dissolved organic material due to export from coastal habitats to shelf and open-ocean waters (Alongi et al. 1989; Duarte and Krause-Jensen 2017; Najjar et al. 2018). Furthermore, the exchange of organic matter within and between coastal habitats in tropical and subtropical regions varies depending on seascape configuration (Bouillon and Connolly 2009; Ricart et al. 2015a; Gullström et al. 2018) and is driven by abiotic conditions such as tidal forces (Slim et al. 1996) and biologically transferred energy through e.g. animal movement or trophic relays (Kneib 1997). Thus, the pathways for carbon transfers within the coastal seascape are many and highly diverse (Bouillon and Connolly 2009; Regnier et al. 2013; Fan et al. 2020), making it a challenge to assess and understand the fate of organic carbon across the land–sea continuum. The exchange of organic carbon between mangroves and seagrass meadows has been verified in a few studies (e.g. Slim et al. 1996; Bouillon and Connolly 2009), but how the strength of this link varies spatially across a coastal landscape influenced by land-use change is less understood.

Continuing degradation and loss of vegetated coastal ecosystems is a significant threat to natural blue carbon sinks globally (Pendleton et al. 2012). The extent of mangroves and seagrass meadows has been reduced by approximately one third during the past few decades as a consequence of human-driven stress (Waycott et al. 2009; Spalding et al. 2010; Mcleod et al. 2011). Such landscape transformation may have profound consequences for carbon dynamics (e.g. Pastick et al. 2017; Trevathan-Tackett et al. 2018). For management purposes, the protection of hotspot carbon-sink environments could contribute to a nature-based solution against increasing greenhouse gas emissions. However, in dynamic seascapes influenced by high-energy water movements, spatial heterogeneity and structurally complex habitats, it is intricate to determine the dynamics and fate of carbon, and hence how to enforce efficient management strategies. In addition, degradation of carbon-rich habitats in the coastal interface may lead to erosion of their sediment (Dahl et al. 2016a) and increased export to the adjacent environments (Duarte and Krause-Jensen 2017). This could in turn alter the source-sink dynamics, exchange patterns and fate of organic matter, which may have ultimate implications for the habitats’ carbon sink capacity (Bouillon and Connolly 2009).

Using a multiscale landscape ecology approach (Turner 1989; Wu and Hobbs 2002), we estimated patterns of variability of sedimentary organic carbon (Corg) storage in seagrass meadows with the overall aim to assess the influence of seascape configuration and land-use change (based on mangrove conversion) on the dynamics and fate of carbon within a coastal land–sea interface. The surveyed area, Tsimipaika Bay in northwest Madagascar, is a diverse mangrove–seagrass coastal embayment, where extensive mangrove deforestation has occurred in localized areas of the bay. Based on a selection of fundamental landscape ecology-derived metrics, we predicted that the sedimentary Corg content (%), stocks (g m−2) and source composition (based on fractionation of δ13C) in seagrass meadows may depend on (a) the spatial arrangement of habitat patches and distances to surrounding seascape features (mangrove forest, water channels, deep water, deforested area and seagrass meadow types), (b) patch heterogeneity of habitats (continuous vs. patchy areas), (c) within-patch attributes (seagrass structural complexity, sediment properties and bioturbating infauna), and (d) functional connectivity (supply and movement) from surrounding seascape features (including also degraded mangrove forest). We hypothesized that (1) broad landscape-scale metrics will have a strong influence on the dynamics and fate of carbon in seagrass sediments, and (2) deforestation of mangrove will increase levels and shift the composition of sedimentary Corg in recent sediments of adjacent seagrass meadows towards being more mangrove-generated.

Material and methods

Study site

The study was conducted in a mangrove–seagrass coastal seascape in the southeastern part of Tsimipaika Bay (previously referred to as Ambanja Bay; Jones et al. 2016) in northwestern Madagascar, encompassing an area of approximate 100 km2 (13°31′26.52″S, 48°24′41.56″E–13°28′8.94″S, 48°29′48.91″E, Fig. 1), between September and December 2016. Seagrass meadows are distributed within a 7 km coastal shallow-water zone (Figs. 1c, 2). The dominating seagrass species in the bay is Enhalus acoroides (Ea), while mixed meadows composed of Cymodocea rotundata/serrulata and Thalassia hemprichii (Cym/Th) (primarily spread on the tidal flats) and meadows dominated by Thalassodendron ciliatum (Tc) (mainly growing further out towards the open ocean) are also common (Fig. 2). These three seagrass communities (i.e. Ea, Cym/Th and Tc) constituted the focal habitats of this study due to their wide distribution (Gullström et al. 2002) and high potential as carbon sinks in the western Indian Ocean region (e.g. Githaiga et al. 2017; Gullström et al. 2018; Juma et al. 2020). The bay also contained vast scattered sandy areas with low densities of intermixed small-sized seagrasses belonging to the genera Halodule spp. and Halophila spp. (Fig. 2), which were not further investigated in the current study on the assumption that they in this high-current environment (being particularly short, thin and sparsely distributed mainly in the intertidal areas) are expected to have limited potential to capture large amounts of allochthonous organic material from other sources. However, it is known from e.g. Australian sites that Halophila ovalis contains considerably high carbon stocks, but the most probable reason for these high levels is that they there thrive in accumulation bottoms and that the seagrass shoots themselves contribute little to the carbon storage capacity (Lavery et al. 2013). Syringodium isoetifolium was also present, but only to a very minor extent (Fig. 2). The area is influenced by a semi-diurnal tidal regime with a mean spring tide of 3.8 m (McKenna and Allen 2005), generating both intertidal- and shallow-subtidal meadows. Hydrodynamic forces in the area are influenced by a complex network of rivers, mangrove creeks and coastal channels between the vegetated zones (Fig. 1). On the landward side of the seagrass meadows, the coastline is fringed by extensive mangrove forests, where the dominant species are Rhizophora mucronata, Ceriops tagal, Avicennia marina, Sonneratia alba and Bruguiera gymnorhiza (Jones et al. 2016; Arias-Ortiz et al. 2020). The Tsimipaika Bay together with Ambaro Bay (north of Tsimipaika Bay) comprises Madagascar’s second most extensive mangrove area (Giri 2011), which since the 1990s has experienced extensive degradation for timber and charcoal extraction (Jones et al. 2016). The loss of mangroves in the entire Tsimipaika-Ambaro Bay complex has been estimated to about 1000 ha from 1990 to 2010 (Jones et al. 2014). Within the study area, substantial mangrove losses have occurred in the southeast, while the southwestern mangrove area is still relatively pristine (Jones et al. 2016, Fig. 2). The seascape area in closest vicinity to the deforested mangrove areas almost exclusively contained Ea-dominated meadows, unvegetated areas and channels, while the major area adjacent to the intact mangroves comprised habitats with a larger variety of seagrass species. The terrestrial primary producers in the areas landward from the mangroves in Tsimipaika Bay consist of a mosaic of grass, bushes, cultivated forests and agricultural areas dominated by rice production (Jones et al. 2016).

Fig. 1
figure 1

Maps showing the location of the study area in Tsimipaika Bay (a, b previously referred to Ambanja Bay) in northwest Madagascar. The encircled area in the satellite image (c) visualizes the study area

Fig. 2
figure 2

Illustrations of average proportion (%) of organic carbon (Corg) in the DW sediment (a), sediment density (g DW cm−3) (b), and accumulated organic carbon stocks (g Corg m−2) (c) in the topmost 25 cm of sediment. Abbreviations: Unveg unvegetated area, Tc is the Thalassodendron ciliatum dominated habitat, Cym/Th mixed habitat dominated by Cymodocea rotundata/serrulata and Thalassia hemprichii, Si is the Syringodium isoetifolium dominated area, Ea is the Enhalus acoroides dominated habitat, Hal spp. is the Halodule spp. and Halophila spp. dominated area, mixed meadow area dominated by undefined seagrass species, deep ocean area with an average depth > 10 m, mangrove intact forested mangrove area, mangrove disturbed deforested mangrove area, SAV undefined submerged aquatic vegetation. The bar plots to the right represent mean (± SE) proportion (%) of sedimentary (DW) organic carbon (Corg) (top graph), sediment density (g DW cm−3) (mid graph), and accumulated sedimentary (DW) organic carbon (g Corg m−2) (bottom graph) in the topmost 25 cm of sediment in the different seagrass habitats, mangrove and unvegetated areas, respectively

Habitat mapping and landscape metrics for correlation analysis

Mapping of habitats was conducted along transect lines (stretching from the mangrove margin to the offshore seagrass meadow edges) that were placed about 50–100 m apart and perpendicular to the shoreline across the bay, either by walking in low tide or by boat using an aquascope in deeper water. Along each transect, the start and end of each habitat were marked (using a handheld GPS) and type of habitat (i.e. unvegetated area, specific seagrass species, channel or mangrove) noted. This information was added to the Google Earth software, where polygons were created using satellite imagery based on data from the field ground-truthing. The polygons were transferred to ArcMap 10.7 (ESRI, Redland), where areas of each habitat were calculated, while habitat areas of intact and deforested mangroves were adopted from Jones et al. (2016). Combining the two sets of habitat data, we produced a single complete habitat map in the ArcMap software. From the final habitat map, spatial metrics, including both seascape- and distance-based measures, were calculated. Seascape-based units of 1 km in diameter buffer zones were set around each seagrass sample site, and within each unit, size, proportion and diversity of the different habitats and patches were calculated (Table 1). Distances to the edges of intact- and deforested-mangrove areas, channels, coastal baseline (parallel to the coastline in an east–west direction across the embayment) and open ocean (> 10 m depth) were measured from the center of each seascape unit (Table 1). The deforested mangrove distance metric was complemented by two additional measures of proximity to deforested mangroves (i.e. distances to midpoints of the eastern and central parts of the deforested area) to allow a spatially representative assessment of the associations to deforestation (Table 1).

Table 1 Metrics, and associated abbreviations, used as predictor variables in the partial least squares (PLS) regression analyses (Table 2) and PCA model (Fig. 4) for organic carbon (Corg) stocks

Sampling of sediment

To assess the spatial dynamics of sedimentary Corg, sediment cores (n = 4) were collected (at least 10 m apart) from a total of 36 sites in selected habitats, including Ea- (n = 12 sites), Cym/Th- (n = 6) and Tc- (n = 6) meadows, mangroves (n = 4) and unvegetated areas (n = 8) across the bay (for locations, see Fig. 2) using push corers (ø = 4.7 cm, h = 60 cm). The core sampling was performed by scientific SCUBA divers in submerged areas or by walking in low-tide water, depending on site location and contemporary tidal level (maximal depth: ~ 6 m). The sediment cores were extruded and sliced into six depth sections, including 0–2.5, 2.5–5, 5–12.5, 12.5–25, 25–37.5 and 37.5–45 cm. At least the sections down core to 25 cm depth could be retrieved in all cores and were later used in the carbon analysis. This is also a standard sediment depth used in several sediment carbon stock assessments (e.g. Lavery et al. 2013; Röhr et al. 2018). To assess any temporal dynamics and dating of the sediments, we collected nine long sediment cores using PVC tubes (ø = 6.2 cm, h = 1.5 m) in seven Ea meadows and two unvegetated sites. These long-core tubes were manually hammered into the sediments, extracted by hand and cut lengthwise. The sediments inside the long corers were sliced at 0.5-cm-thick intervals throughout the first 20 cm, and at 1-cm-thick intervals below this depth.

The bottom edge of the corers were sharpened to reduce sediment compression and to make it easier to shred the seagrass roots and rhizomes. Sediment compression was accounted for in all cores by measuring the distance from the top of the core to the sediment surface, inside and outside the corer after being inserted into the sediment (Glew et al. 2002). From this, the compression-correction factor (%) was calculated for each core following Howard et al. (2014), for which the sediment volume of the slices in the core was adjusted. Average compressions were also calculated (short cores: 21 ± 2.5%, long cores: 20 ± 5%, mean ± SE).

Seagrass biometrics and infauna collection

To assess the influence of seagrass structural complexity and bioturbating infauna on sedimentary Corg variables (i.e. content, stocks and source composition), seagrass meadow plant biometrics and abundance and biomass of infauna were recorded randomly in the vicinity to each core sampled in the seagrass meadows. Mean seagrass canopy height (based on maximum shoot lengths, n = 20) was measured with a measuring stick. Mean seagrass shoot density was estimated in 0.25 × 0.25 m quadrats (n = 5). Total cover (%) of seagrass as well as cover of the dominant seagrass species were estimated in 0.5 × 0.5 m squares (n = 10). Infauna was collected with a sediment push corer (ø = 10 cm) down to about 10 cm sediment depth and then sieved (1 mm mesh), counted (abundance) and weighted (wet weight biomass).

Sediment density and carbon analysis

Sediment slices from the cores were weighed wet, homogenized (by stirring) and a subsample of 60 mL (short cores) or the whole section (long cores) from each depth was dried at 60 °C for approximately 48 h until constant weight. For each sediment slice, dry bulk density (g DW cm−3; hereafter referred to as sediment density) was calculated and adjusted to the compression corrected volume for each core slice. Sediment samples were sieved (1 mm) to exclude large living plant material and shells before they were ground to a fine powder using a mixing mill (Retch 400) for subsequent carbon analysis.

The carbon content, i.e. % C of the DW sediment from both short- and long-cores, was analyzed using a carbon–nitrogen elemental analyzer (Flash 2000, Thermo Fisher Scientific, Waltham, MA, United States, precision error < 0.02% for carbon). Each sediment sample was divided into two subsamples for analysis of the total carbon (Ctot) and Corg contents (%). Prior to Corg analysis, the subsamples were treated with 1 M HCl (by direct pipetting) to remove inorganic carbon and dried for about 24 h. The inorganic carbon content was then determined from the difference between the Ctot and Corg contents. The weighted average of % Corg content down to each specific depth in the sediment was determined. Accumulated areal organic carbon stocks (g Corg m−2) were estimated by integrating the Corg density (i.e. Corg multiplied by the compression-corrected sediment density) with depth above the upper 25 and 100 cm in the short and long cores, respectively.

Stable isotope analysis

Bulk stable isotope signals of carbon (δ13C) and nitrogen (δ15N) were analyzed (at UC Davis Stable Isotope Facility, CA, USA) in sediment from 16 sites, including meadows dominated by either Ea (n = 8), Tc (n = 3) or Cym/Th (n = 2) as well as unvegetated sites (n = 3) at different locations within the bay. In total, 48 cores (three from each site) were analyzed for stable isotope signals. The measures were conducted in each sediment depth interval of each core and were thereafter used to calculate the weighted average for each core. Material for potential Corg sources (n = 4 from each source) was collected, including samples of above- and below-ground seagrass biomass (from seagrass meadows dominated by the different species), mangrove- and terrestrial-plant material and suspended particulate organic matter (SPOM) from offshore waters. SPOM was collected by filtering 1.5 L of water onto pre-combusted GF/F-filters and the plant material was washed, cleaned from epiphytes, dried and ground before analysis. To remove inorganic carbon prior to analysis, the sediment, seagrass biomass and SPOM were treated with low concentrated acid (1 M HCl) to minimize any δ15N fractionation according to recommendations in Kennedy et al. (2005). The isotope signals were analyzed with a micro cube elemental analyzer (Elementar Analysensysteme GmbH, Hanau, Germany) interfaced with an isotope ratio mass spectrometer (EA-IRMS). The samples were interspersed with four known laboratory references (bouvine liver, glutamatic acid, enriched alanine and nylon 6), which are calibrated against international reference materials (Sharp 2005), and the temporary isotope ratio generated is thus relative to the reference gas peak analyzed with each sample. These temporary ratios are consequently corrected for the entire batch based on the known laboratory references, generating the final ratios. The samples’ isotope ratios were expressed as δ values in parts per thousand (‰) relative to VPDB (Vienna Pee Dee Belemnite) for carbon and air for nitrogen (Sharp 2005). The long term standard deviation for the EA-IRMS is ± 0.2 ‰ for 13C and ± 0.3 ‰ for 15 N.

Specific 210Pb activities were measured in the upper 40 cm of six of the long cores from Ea-dominated meadows to assess potential changes in sediment accumulation rates and dynamics with time (i.e. last century to decades). For details of methods and results, see supplementary material (“210Pb analysis for sediment dating”). Total 210Pb and 226Ra specific activities (derived by alpha- and gamma-spectrometry, respectively) were within error of one another from the surface to the bottom of the analyzed cores, confirming agreement between alpha and gamma methods but resulting in negligible excess 210Pb concentrations, and thereby non-datable sediment profiles (Arias-Ortiz et al. 2018).

Data analysis

To elucidate the relative importance of the different predictors (i.e. distance-, seascape-, plant-, and sediment-metrics; Table 1) for the Corg content (%) in the topmost 25 cm of the sediment in the studied seascapes as well as for % source of origin of Corg (evaluating distance metrics only), the data were analyzed with projection to latent structures models by means of partial least squares (PLS) regression analysis (Wold et al. 2001) on log10(x + 1)-transformed data using the software SIMCA-P + 15.0.2 (UMETRICS). PLS analyses were conducted on data including all seagrass seascapes as well as on data from seascapes of each seagrass community (i.e. Ea, Cym/Th and Tc) separately. Principal component analysis (PCA) was performed in SIMCA-P + 15.0.2 (UMETRICS) to explore and visualize (in a biplot) the potential associations between the explanatory metrics, the sedimentary Corg content and the different seagrass sampling sites.

Variation in the sedimentary Corg content (%), sediment density and accumulated Corg stocks (g m−2) among the different seagrass habitats (i.e. Ea, Cym/Th and Tc) and unvegetated areas was tested with one-way ANOVAs. The mangrove samples were excluded in the statistical comparison due to the limited sample size (since they were not in focus in the current study). A posteriori multiple comparisons were conducted using Tukey’s HSD test. Prior to all analyses, data were tested to meet the assumption of homogeneity of variances (Levene 1960). To be able to compare the means in Corg among habitat types, we assessed if the sampled sites were randomly selected (not spatially autocorrelated) using a spatial autocorrelation model (Morans I) on the %Corg data (0–25 cm depth) for the 36 study sites in ArcGIS v. 10.7 (Supplementary Fig. S2). Variation in Corg content, sediment density and accumulated Corg stocks in the long sediment cores were compared between Ea-dominated meadows and unvegetated sites using the non-parametric Kruskal–Wallis test (since with the log10(x + 1) transformation, homogeneity of variance was not achieved). All univariate statistics were performed using Statistica v. 13. The sources contributing to the sedimentary Corg (based on the stable isotope signals) were analyzed with the Stable Isotope Analysis (SIAR) package in R (v. 3.5.3). The model was based on original concentrations of carbon and nitrogen and the discrimination factors used were 3‰ for δ15N and 1‰ for δ13C to cover the ranges of output signals in seagrass ecosystems as summarized by Lepoint et al. (2004). As the mangrove and terrestrial material had overlapping δ13C signatures (Lamb et al. 2006), these two sources were pooled prior to analysis, and the specific sources used in the final analysis were seagrass plants, mangrove/terrestrial material and SPOM (see Supplementary Fig. S1 showing a biplot of the δ15N and δ13C signals of the sources used in the mixing model).

Results

Sediment carbon content, bulk density and carbon stocks

The sedimentary Corg content (%), sediment density (g DW cm−3) and accumulated Corg stocks (g Corg m−2) did clearly vary in the seagrass surface sediments (0–25 cm depth) across sites in the embayment (ranges: 0.16–1.18% Corg, 0.49–1.32 g DW cm−3, 496–2383 g Corg m−2, respectively, Fig. 2). There were also significant differences at habitat level (comparing Ea, Cym/Th, Tc and unvegetated areas) (ANOVA; Corg content: F3,28 = 15.6, p < 0.001; Sediment density: F3,28 = 14.2, p < 0.001; Corg stocks: F3,28 = 3.5, p < 0.05; Fig. 2). Ea sediments contained higher Corg content (Tukey’s HSD, p < 0.01) while lower sediment density levels (Tukey’s HSD test, p < 0.05) compared to the other seagrass habitats and unvegetated areas. For the Corg stocks, there were no significant differences among the seagrass communities. Unvegetated sediment contained lower Corg content than Ea- and Tc-sediments (Tukey’s HSD, p < 0.01), but did not differ from Cym/Th sediments. The sedimentary Corg contents and stocks in the fringing mangrove habitats varied substantially, while the mean level was manifold higher than those in the seagrass- and unvegetated-habitats (Fig. 2a, c). The mean content of inorganic carbon in the sediment ranged from ~ 12% in mangroves to ~ 40% in the habitats dominated by Ea or Cym/Th and ~ 60% in Tc-dominated meadows and unvegetated areas.

Considering the upper meter of sediment (from long cores; Fig. 3), Corg content varied among the Ea sites and downcore the sediment of individual cores (Fig. 3a). However, on average the Corg content was relatively constant with depth (Fig. 3a). The average Corg content and stocks in the upper meter of Ea sediments (content: 0.89 ± 0.04%, Fig. 3a; stocks: 86 ± 11 Mg C ha−1, mean ± SE, Fig. 3b) were significantly higher than in unvegetated sites (content: 0.35 ± 0.03%, Fig. 3a; stocks: 46 ± 11 Mg C ha−1) (Kruskal–Wallis, p < 0.01).

Fig. 3
figure 3

Organic carbon (Corg) content (%) profiles in the upper meter of sediments in Enhalus acoroides (green) and unvegetated (yellow) areas (a), top meter Corg stocks (Mg Corg ha−1, mean ± SE) including the measures from each site (b) and the locations where the long cores were samples in Tsimipaika Bay (c)

Total 210Pb specific activities in the upper 40 cm of Ea sediment were constant with depth (range: 19.1 ± 1.2 to 85.5 ± 4.5, Bq kg−1) and agreed with 226Ra specific activities (range: 29 ± 2 to 106 ± 8, Bq kg−1) (Supplementary Fig. S2). The lack of excess 210Pb in the six long cores of Ea sediment suggests negligible recent net accumulation of sediment.

Influence of scale-dependent metrics on sediment Corg content

Overall PLS-modelling performance

The PLS regression models showed that all metric categories (i.e. distance-, seascape-, plant- and sediment-metrics) were represented among the foremost predictors of sedimentary Corg content (%) in the topmost 25 cm of sediment, both when including all seagrass seascapes and when specifically focusing on seascapes of each of the three seagrass communities (i.e. Ea, Cym/Th and Tc; Table 2). The Q2-statistics for the different seagrass seascape models showed a cross-validated variance ranging from 20 to 46% (i.e. higher than the defined significance limit of 5%) and hence the models displayed good predictability (Table 2). The cumulative portion of all predictors combined (R2y) explained 57–88% of the variance of the Corg content, which shows that the models displayed high level of determination (Table 2). In the PLS models, between 10 and 13 of the 27 predictors tested had a variable of importance (VIP) value over 1 and thus contributed significantly to the respective model (for ranking order, see Table 2).

Table 2 Summary of partial least squares (PLS) regression models ranking the most important variables predicting Corg content (%) in the topmost 25 cm of sediment for the seagrass seascapes in general and for community-specific seagrass seascapes.

Influence of plant- and sediment-metrics

Seagrass structural complexity measures in terms of canopy height (positively correlated) and shoot density (negatively correlated) were among the most important predictors explaining the sedimentary Corg content in the models for all seagrass seascapes and Ea-dominated seascapes (Table 2). In contrast to these two models, canopy height was negatively correlated to Corg in Tc-dominated sediment (Table 2). Among the plant metrics, the various seagrass cover predictors also contributed to the performance of some models, although to a lesser degree than canopy height and shoot density (Table 2). Besides the importance of plant structure measures as within-meadow predictors of Corg, sediment metrics in terms of sediment density (most important), infaunal abundance and infaunal biomass (all negatively correlated to sedimentary Corg content) also contributed (in various ways depending on seascape model) to the performance of the different models (Table 2).

Fig. 4
figure 4

Principal component analysis (PCA) biplot displaying correlative assessments of scale-dependent predictors (distance-, seascape-, plant- and sediment-metrics), sedimentary Corg content at different depth sections and the model observation metrics of each seagrass seascape represented by Enhalus acoroides, Cymodocea spp./Thalassia hemprichii and Thalassodendron ciliatum in Tsimipaika Bay. For abbreviations, see Table 1. PC1: R2x = 0.317, PC2: R2x = 0.133, Q2cum = 0.117 (i.e. > significance limit of 0.07)

Influence of distance- and seascape-metrics

Spatial pattern analysis based on distance- and seascape-metrics demonstrated significant relationships between several of the examined predictors and the level of Corg content in seagrass sediment (Table 2). A variety of seascape metrics measuring seascape proportion of different habitats, focal area size and patch heterogeneity (i.e. number, richness, mean size and edge of patches) were related to sedimentary Corg content (Table 2). Seascape proportion of Ea had a strong positive influence on sedimentary Corg levels in seagrass seascapes in general and in Ea-dominated seascapes (Table 2). Seascape proportion of the two other seagrass communities (i.e. Cym/Th and Tc) and focal area size of seagrass also seem to be important for the sedimentary Corg content, although the strength and nature of the relationships showed community-specific variation (Table 2). For instance, opposite to the positive influence of seascape proportion of Ea on sedimentary Corg, seagrass seascapes in general showed a negative relationship between seascape proportion of Cym/Th or Tc and Corg content in the sediment (Table 2). The significant patch-related seascape metrics showed mixed influences on sedimentary Corg content (Table 2). In particular, edge was a strong predictor, with an opposite influence on Corg content in sediment of Cym/Th-dominated seascapes (negative relationship) and seascapes dominated by Tc (positive relationship).

The distance metrics showed community-specific relationships to sedimentary Corg content in seagrass meadows of the bay (Table 2). For instance, distance to deforested mangrove (including the three distance metrics to deforestation, Table 1) were negatively related to Corg content in general seagrass seascapes and those dominated by Ea, while positively in Cym/Th-dominated seascapes and those dominated by Tc (Table 2).

Multivariate PCA model

The PCA-biplot visualized clear community-specific clusters grouping the different seascapes in relation to the scale-dependent predictors (i.e. distance-, seascape-, plant- and sediment-metrics) and all variables of sedimentary Corg content (Fig. 4). Generally, sedimentary Corg content was positively coupled to the model observation metrics of the Ea-dominated seascapes, while the seascapes dominated by Cym/Th or Tc showed opposite patterns, with negative relationships to the Corg content (Fig. 4). The community-specific clusters show that seagrass canopy height was closely and positively related to the Ea-dominated seascapes, while seagrass shoot density was clearly positively related to seascapes dominated by Cym/Th or Tc (Fig. 4). Particular distance metrics of importance (for PC1) were proximity to deforested mangrove and distance to open ocean (Fig. 4). Closer proximity to deforested mangrove and larger distance to open ocean seemed to influence the model observation metrics of the Ea-dominated seascapes, while larger distances from the deforested mangrove seemed to particularly influence the metrics of the Cym/Th-dominated sites (Fig. 4). Sediment depth represented by the different depth sections of organic carbon content (%Corg) and accumulated stocks (AccCorg) were highly clustered and therefore had minor influence on the model performance (Fig. 4).

Source and fate of organic material

The isotopic signals did not vary among the different seascape habitats (i.e. Ea, Cym/Th, Tc and Unveg) (Fig. 5a), while the landscape context showed an influence (Fig. 5b). Overall closeness to mangrove (deforested areas in particular) generated a strong mangrove/terrestrial organic material signal, while the seagrass contribution of Corg in the sediment seemed to be larger towards the open ocean. A significant PLS regression model (Q2-statistics = 24%, R2y = 53%) between distance metrics and organic material originating from the mangrove/terrestrial source in the sediment showed that distance to open ocean (ranked 1, positively correlated) and distance to channel (ranked 2, negatively correlated) were significantly (VIP values > 1) related to the variation in contribution of mangrove/terrestrial organic material to the sedimentary Corg pool (Fig. 5c).

Fig. 5
figure 5

Mean of relative contributions in percent (± SE) of organic material originating from mangrove/terrestrial environments, suspended particulate organic matter (SPOM) or seagrass in the different seagrass-dominated seascapes (a). Relative contributions (%) of organic material from seagrass plants and imported Corg, from mangrove/terrestrial vegetation and SPOM present in the bulk sediment of seagrass meadows determined by isotope signals of carbon (δ13C) and nitrogen (δ15N) (b). Summary of partial least squares (PLS) regression model ranking in a descending scale from left to right the most important distance metrics predicting the percent organic material originating from mangrove or other terrestrial source in the seagrass seascapes. Numbers above/below bars present variable influence on the projection (VIP) values. Brown bars represent predictor metrics with VIP values above 1, i.e. metrics that contribute significantly to model performance (c). For habitat abbreviations, see legend of Fig. 2

Discussion

We found that seascape configuration and mangrove deforestation had a strong influence on Corg content and source composition in coastal sediments of the studied land–sea continuum. Our findings suggest that the carbon sink capacity of seagrass meadows is driven by factors operating at multiple spatial scales (from m to 10 s km). Landscape configuration (spatial arrangement of habitat patches and geographic distances to surrounding seascape features), patch heterogeneity (degree of continuity) and within-patch attributes (seagrass structural complexity and sediment properties) were all metrics of importance for the models’ performance. The outcome revealed a strong land-to-sea gradient, possibly forced by hydrodynamically driven exchange of organic matter across the bay landscape which may be enhanced by recent (last decades) escalation in export of organic carbon from deforested mangroves, with higher Corg levels and a clear mangrove/terrestial signal found in surface sediment (0–25 cm) of seagrass seascapes located closer to deforested mangrove area. The substantial input of allochthonous carbon derived from adjacent habitats to the seagrass meadows emphasizes the importance of cross-habitat connectivity within the coastal seascape. Lateral transfer of organic carbon from mangrove areas (cleared and intact) and increased carbon content in the surface sediment may, however, not necessarily lead to accumulation of carbon stocks, as we observed no recent net accumulation of sediments (where the long cores were collected). This suggests well mixed sediments, which commonly occur in high-energy environments (Arias-Ortiz et al. 2018). Of particular importance for conservation of blue carbon sinks may be this study confirms that seascapes comprising large continuous seagrass meadows had higher Corg levels in comparison to more diverse or patchy seascapes.

Relative importance of landscape configuration, patch heterogeneity and within-patch attributes

A multitude of scale-dependent metrics were important predictors influencing spatial patterns and variability of Corg content in seagrass sediment. While spatial heterogeneity of landscapes is fundamentally scale-dependent (Wu 2004), determining the relative importance of drivers at multiple spatial scales can be very challenging to unravel, especially given the broad range of process-oriented patterns in nature. Our study showed that the sedimentary carbon content varied spatially from low to moderate levels across seascapes in the embayment. The Corg levels were lower than the global mean of ~ 2% (Fourqurean et al. 2012), but comparable to other tropical and subtropical seagrass areas (e.g. Fourqurean et al. 2012; Lavery et al. 2013; Phang et al. 2015), including studies from the western Indian Ocean (Dahl et al. 2016a; Githaiga et al. 2017; Belshe et al. 2018; Gullström et al. 2018, 2021). The top meter of sedimentary Corg stocks (i.e. from the Ea-dominated meadows) is comparable to estimates from other subtropical regions (72–90 Mg Corg ha−1; Miyajima et al. 2015; Serrano et al. 2019) and two times higher than Corg stocks observed in unvegetated sediments.

The analysis of spatial attributes demonstrated that landscape configuration and position within the coastal landscape strongly influence sedimentary Corg in seagrass meadows. In accordance, Ricart et al. (2020) found that sedimentary Corg stocks were related to landscape configuration (based on seagrass meadow type) in a coastal estuary. In general, a higher seascape proportion of Ea positively influenced the sedimentary Corg content. Thus, the larger the area covered with the focal seagrass, the more Corg was found in the sediment, especially in the Ea-dominated seascapes. In the Cym/Th-dominated seascapes, on the other hand, the proportion of unvegetated area was positively related to sedimentary Corg content, which could be explained by their location in vast sand flats and the patchy occurrence of this seagrass community in the bay. Any occurring organic material (autochthonously- or allochthonously-produced) will most likely be trapped and end up in the seagrass patches rather than in the surrounding unvegetated areas. In fact, high patch richness had generally a negative influence on the sedimentary Corg content, even though this picture was not evident in the community-specific models, where such relationships might be masked by the fact that many habitat types may be missing in the seascape units (i.e. the 1 km buffer zones). We observed that larger continuous seagrass areas may have a greater potential to contain higher amount of Corg than more diverse or patchy areas, as the autochthonous plant material is more likely to end up within a vast meadow (Ricart et al. 2017). Such a habitat–area relationship has previously been suggested in the western Indian Ocean region (Gullström et al. 2018). The edges of habitat patches measured as perimeter–area relationships also seem to have an important but contrasting (species- or community-specific) influence on the levels of Corg in seagrass sediment, as seen when comparing the models of the Cym/Th- (negative relation) and Tc- (positive relation) dominated seascapes. The model performance of the Cym/Th-dominated seascapes hence harmonizes with the continuous seagrass area assumption (emphasized by Gullström et al. 2018), as a lower perimeter–area relationship seems to generate a higher Corg content in the sediment. Patches with high perimeter–area relationship may indicate fragmentation into small patches, which may negatively influence the sedimentary Corg content in seagrass meadows (Ricart et al. 2015b, 2017; Oreska et al. 2017; Gullström et al. 2018). The opposite influence seen in the Tc-dominated seascapes cannot be explained by the same assumption, but rather be a result of the position and shape of Tc-meadows as narrow strips at the edge of the seagrass zonation towards the open sea and/or the large perimeter–area ratio of surrounding habitats in the seascape units.

The most important distances-based metrics were found along the land-to-sea gradient, especially those assessing the association to deforested mangrove. In general, the models revealed that the sedimentary Corg content in the upper 25 cm of sediment was higher in closer proximity to the deforested area and lower closer to the open ocean. Such land–sea linkages have recently been studied in estuarine environments (Ricart et al. 2020), with the findings showing a high proportion of land-based organic carbon in seagrass sediment. The land–sea coupling seen in our study can have several explanations that act in combination. First, outwelling of organic matter and Corg from cleared mangrove areas may have occurred. Second, the complex network of rivers, mangrove creeks and channels in the vicinity of the deforested mangrove may increase the transport of organic material from land and from the cleared mangroves. Third, the seagrass seascapes in the area adjacent to the cleared mangrove are dominated by Ea, where higher sedimentary Corg levels were observed. Fourth, the zone between the intact mangrove and the open ocean was slightly narrower compared to the zone outside the deforested mangrove, wherefore seagrass seascapes in this area may be influenced by closer proximity to both mangrove and open ocean. Recently, Arias-Ortiz et al. (2020) showed that 20% of the soil carbon stocks in the upper meter of cleared mangroves in Tsimipaika Bay have been lost in the last 10 years due to deforestation through several pathways, including CO2 emissions and DOC export, which may have implications for adjacent ecosystems such as seagrass meadows. We could, however, not observe any influence of proximity to cleared mangrove area on carbon stocks (down to one meter depth) of the seagrass sediment. One reason may be that mangrove deforestation is a relatively recent disturbance (occurring since the 1990s in the region), hence the integration of Corg stocks in the upper meter may hide any difference in Corg inputs during recent times, which are expected to have an effect in the upper layers only. Additionally, we did not observe recent net sediment accumulation at these sites as shown by the 210Pb data (Supplementary Fig. S2), suggesting that these long seagrass cores are not representative of recent sediment dynamics. The lack of excess 210Pb is consistent with intense sediment resuspension and transport, usually observed in high-energy environments generating well-mixed sediments (Arias-Ortiz et al. 2018). The sampling location of some of the long cores may be a reason for a low Corg content; we speculate that the factors “closeness to channels” (with higher water flow) and “focal area size” could downplay the factor “closeness to mangrove”. But most likely, the lower number of long cores together with the high variability in Corg content observed within them, might mask some of the patterns observed in the short cores, the number of which is 15-fold larger and cover a more widespread area across the bay. However, we did see clear differences in Corg content and stocks between the vegetated and unvegetated sediments in both long and short cores, suggesting that seagrasses may play an important role in preserving the Corg already present in the sediments.

Considering the influence of plant traits, our findings suggest that seagrass structural complexity and cover influence the plants’ efficiency to sequester or trap allochthonous organic material moving within the seascape. The PLS and PCA models showed that seagrass canopy height (mainly driven by the variation in shoot height of Ea) and cover positively influenced the sedimentary Corg content and that many seagrasses may have strong effects on carbon sequestration through their high primary productivity. This was most likely due to high seagrass plant biomasses and their effecient ability to trap material from surrounding areas (Samper-Villarreal et al. 2016). Sediment density seems also to have a significant role for Corg levels in seagrass habitats, as shown by the models in this study and by earlier studies elsewhere (e.g. Avnimelech et al. 2001; Dahl et al. 2016b, 2020a; Röhr et al. 2018). High Corg content is generally associated with low sediment density since organic material has a lower weight compared to minerogenic particles (Winterwerp and van Kesteren 2004) and a higher water content (Avnimelech et al. 2001). Other sediment characteristics, such as mud content and degree of grain size sorting, are known to positively influence the sedimentary Corg content and inversely affecting sediment density (Röhr et al. 2018). Both grain size distribution and sediment density are linked to the hydrodynamic exposure, where sheltered areas have lower sediment density and higher mud content (Dahl et al. 2020b) as these sediment bed types are more easily eroded (Jacobs et al. 2011). Sediment density in our study was negatively correlated to sedimentary Corg in the seagrass seascapes dominated by the higher growing climax species (i.e. Ea and Tc), which was potentially due to the higher organic matter content in the sediment of these larger seagrass species compared to the Cym/Th-dominated sediment, and that larger species may have a higher potential for reducing the water flow with the canopy. The negative influence of seagrass shoot density on the Corg content may seem contradictory since increasing shoot densities normally slows the water flow (Peterson et al. 2004) but is likely due to the inverse relationship between shoot density and canopy height (seen e.g. in Gullström et al. 2006). This means that the Corg content, especially in Ea-dominated meadows, was higher in areas with high seagrass canopies while low shoot density. In addition, the higher Ea canopy seen outside the deforested mangrove area, where the water often was more turbid, may be a plant response to limited light conditions. In contrast to Ea-dominated meadows, very dense Tc-dominated meadows may limit allochthonous carbon to reach the underlying sediment. Overall, our study confirms that both landscape processes and plant–sediment properties of seagrass habitats play an important role for carbon accumulation in the sediments, and thus for these areas’ carbon sink function.

Cross-habitat connectivity and influence of altered land-use

In coastal environments, organic matter is regularly transferred across habitat boundaries (Bouillon and Connolly 2009; Hyndes et al. 2014), which depends on different direct and indirect mechanisms of carbon transfer at seascape scale (Huxham et al. 2018). For instance, exchange of carbon (e.g. litter, particulate and dissolved organic and inorganic carbon) between vegetated coastal/terrestrial areas and seagrass meadows can be substantial (Maher et al. 2013). Global estimates show that 50–72% of the organic carbon pool of the surface sediment in seagrass habitats can originate from terrestrial vegetation or other allochthonous sources (Gacia et al. 2002; Kennedy et al. 2010). Our study showed a high mean proportion of allochthonous organic carbon originating either from mangrove/terrestrial areas or as SPOM that exceeds 70% in the seagrass sediment, although the mean levels varied among seascapes (from 10 to 59%) and across seagrass community types (71% for all seagrasses, 77% for Ea, 68% for Cym/Th and 58% for Tc). Such high contribution of allochthonous material to carbon accumulation and storage has been seen in other tropical seagrass areas (e.g. Chen et al. 2017; Ricart et al. 2020). Carbon subsidies from mangroves to other shallow-water habitats (such as seagrass meadows) can occur due to multiple interconnecting factors, such as rainfall creating runoff (Huxham et al. 2018), hydrodynamic and tidal forces (Slim et al. 1996), animal movement (Kneib 1997) and/or major plant systems’ efficiency in trapping carbon (Kennedy et al. 2010). The influence of seascape configuration may be changed with major alterations in land use, as has occurred in an extensive part of our studied embayment. The variability in compositions of Corg from the isotopic analysis shows that transfer of carbon from intact mangroves to seagrass meadows differed from the observations in seagrass meadows next to the deforested areas, where the major source of carbon was mangrove-originated. Patterns of carbon movement and exchange in this coastal environment seem to be driven by a gradient across the land–sea interface, with higher presence of mangrove/land-originated carbon in seagrass meadows at closer distance to mangrove (both intact- and deforested-areas) and farther from the open sea. Even though studies in areas where mangrove and seagrass occur next to each other have showed high carbon subsidies from mangrove to seagrass areas (e.g. Gonneea et al. 2004; Kennedy et al. 2010; Gullström et al 2018), our results indicate that clearing of mangroves increases such outwelling effect, inducing a shift in carbon transfer from mangrove areas to adjacent seagrass meadows. Further, our results also show that closeness to channels can reduce sedimentary Corg in seagrass sediment, which is most likely due to erosion by increased hydrodynamic forces (Dahl et al. 2018).

Concluding remarks

This study emphasized that the carbon sink function of tropical seagrass meadows is driven by multiscale spatial determinants, including landscape configuration, patch heterogeneity, within-patch attributes (plant–sediment properties), land-use change and land-to-sea gradients.

Such information regarding spatial heterogeneity is of fundamental importance to understand complex patterns and processes in ecology (Pickett and Cadenasso 1995). In blue carbon science, understanding spatial patterns of variability and what is driving the observed patterns is important to successfully identify carbon sink hotspots and to avoid double-counting in order to develop conservation prioritizations and blue carbon projects.