Testing biodiversity-ecosystem function relations in nearshore marine sediments
Abstract
Rapid biodiversity loss is raising concerns about potential declines in the functioning of Earth’s ecosystems. Although decades of research explore biodiversity-ecosystem function (BEF) relations, empirical BEF assessments have lagged behind theoretical advances, particularly in marine benthic systems. Here, we incubate intact sediment push cores to examine the relationship between macroinfaunal community composition and benthic nutrient cycling in three nearshore sub-Arctic sites. First, we quantitatively assess potential effects of taxonomic and functional diversity, as well as community-weighted trait means, on oxygen and nutrient fluxes. Second, we examine fluxes in relation to macrofaunal abundance, oxygen consumption (a proxy for total core metabolism), and abundance of key functional groups, to test fundamental expectations based on ecological theory. We report distinct macrofaunal communities and contrasting benthic fluxes among sites, with oxygen and ammonium largely driving multivariate inter-site flux differences. Diversity indices and community-weighted trait means collectively explained ~76% of the variation in multivariate fluxes but provided little insight into the mechanistic links between diversity and functioning. In fact, we find that total macrofaunal abundance was the key driver of benthic fluxes at our sites, rather than functional community composition, which may have important implications for coastal conservation planning. Overall, our findings emphasize the highly context-dependent nature of BEF relationships and highlight the need to improve empirical understanding of these patterns in complex, natural ecosystems.
Data were collected from three intertidal soft-sediment sites around the island of Newfoundland, Canada, using manual sediment push cores. Nutrient (ammonium, silicate, phosphate, and combined nitrate/nitrite) and oxygen fluxes were measured during 24-hour incubations of sealed sediment cores. Cores were then sectioned into depth layers (0-2 cm, 2-5 cm, and 5-10 cm), and macrofauna were sieved through a 500 micron mesh, enumerated, and identified to the lowest possible taxonomic level.
This dataset includes:
- Raw abundance data and functional trait information, as well as nutrient and oxygen concentrations measured in each core;
- Taxonomic and functional diversity indices (including Hill numbers) calculated using common methods in the field;
- Nutrient and oxygen flux rates calculated using common methods in the field;
- Code used to produce the figures, tables, and analyses in Clinton et al. (2025).
Description of the data and file structure
File structure
Code: This folder contains all scripts needed to produce the figures, tables, and analyses in Clinton et al. (2025).
Data: This folder contains raw and processed data (.csv files).
Description of raw data
A core x species matrix that contains raw abundances for each taxon identified in each sediment core. Column headers Alitta.sp to Tellinidae.indet are taxon names following binomial nomenclature (e.g., <Genus>.<species>). The column core_id follows the format: <site code><sampling date id>_<enrichment type>_<replicate id>. Note that enrichment type is unrelated to the analysis conducted in Clinton et al. (2025).
Contains community-weighted mean trait values for each core. For details on how these values were calculated, see Clinton et al. (2024). The column core_id follows the format: <site code><sampling date id>_<enrichment type>_<replicate id>. The 3 possible values for site are: NS = Newman Sound, NH = Neddie's Harbour, and SP = St. Paul's Bay. All remaining columns (diet_Hb to log_body_size_scaled) contain unitless values between 0 and 1. Note that the traits Diet, Feeding guild, and sediment reworking are each divided into several modalities, which are denoted in the column name using an underscore, as seen below:
Diet: _Hb = Herbivore, _C = Carnivore, _Dt = Detritivore.
Feeding guild: _Pr = Predator, _Gr = Grazer, _Sc = Scavenger, _Ff = Filter feeder, _SD = Surface deposit-feeder, _SSD = Sub-surface deposit-feeder.
Reworking: bd = Biodiffusor, sm = Surficial modifier, updown = Upward/downward conveyor, epi = epipelagic
log_body_size_scaled: Log transformed maximum body size (mm), scaled between 0 and 1
Contains taxonomic diversity indices for each core. For details on how these values were calculated, see Clinton et al. (2024). Column names can be interpreted as follows:
site: Name of the site at which that sediment core was collected.
core_id: Unique identifier for each sediment core. Follows the format: <site code><sampling date id>_<enrichment type>_<replicate id>.iet
shannon: Shannon-Weiner index
simpson: Simpson's diversity index
richness: Taxonomic richness (i.e., the number of unique taxa observed within a given core)
pielous_evenness: Pielou's evenness index. Value of 1 indicates perfect evenness (i.e., all taxa are equally abundant), value of 0 indicates maximum unevenness (one taxon dominates completely).
Contains functional diversity indices for each core. For details on how these values were calculated, see Clinton et al. (2024). Column names can be interpreted as follows:
Nb_sp: Number of unique taxa observed within a given core.
Tot_weight: Number of organisms observed within a given core.
FRic: Functional richness
FDiv: Functional divergence
FEve: Functional evenness
FDis: Functional dispersion
FSpe: Functional specialization
FOri: Functional originality
A species x trait matrix, which includes trait information for each taxon. Also outlines the functional group (FG) to which each taxon was assigned, based on these traits. Each functional group and trait is described in detail below. Note: The traits Diet and Feeding guild are fuzzy coded (values between 0 and 1) and divided into several modalities, which are denoted in the column name using an underscore. For example, the trait diet has three modalities: Hb (Herbivore), C (Carnivore), and Dt (Detritivore). The codes used to define each modality are also outlined below:
Diet: diet_Hb = Herbivore, diet_C = Carnivore, diet_Dt = Detritivore.
Feeding guild: _Pr = Predator, _Gr = Grazer, _Sc = Scavenger, _Ff = Filter feeder, _SD = Surface deposit-feeder, _SSD = Sub-surface deposit-feeder.
Reworking: bd = Biodiffusor, sm = Surficial modifier, updown = Upward/downward conveyor, epi = epipelagic
Body size: Values in the body size column are the maximum recorded size (in mm) for that taxon.
Mobility: Value between 0 and 1 indicating the degree of mobility exhibited by the taxon within the sediment matrix. None (value: 0) = Sedentary or only moving within a fixed tube structure. Low (value: 0.3) = Limited free movement (e.g., withdraws into sediment when disturbed). Moderate (value: 0.6) = Slow movement within the sediment matrix via non-permanent burrow formation. High (value: 1) = Freely mobile within sediment in a permanent, excavated burrow system. Outlines which taxa were assigned to which functional groups
FG: free_move_bd_Hb_Dt = free-moving, biodiffusing herbivores/detritivores, hunters = free-moving, biodiffusing carnivores, other = other, copepod = small epipelagic copepods (referred to as "swimmers" in Supplementary material of Clinton et al. 2025), burrowers = surficial modifiers with permanent burrow systems, low_mob_sm = low-mobility surficial modifiers, snails_and_co = free-moving surficial modifiers, heads_down = sub-surface deposit feeding upward conveyors, heads_up = surface deposit feeding downward conveyors
Contains site-level information, as well as variables collected during the sediment core incubations. Column headers are defined below:
site_name_full: Collection site. Note that Alien Cove refers to Newman Sound
site: Abbreviation for collection site. AC = Alien Cove (Newman Sound), SP = St. Paul's, NH = Neddie's Harbour
lat: Latitude of collection site.
long: Longitude of collection site.
substrate: Categorical description of substrate at collection site (sand, mud, or gravel)
bottom_temp: Water temperature (degrees Celsius) at the site during collection.
salinity: Salinity (ppm) at site during collection.
day: Day on which collection occurred.
month: Day on which collection occurred.
year: Year in which collection occurred.
settlement_time: Number of hours cores were left undisturbed between collection and incubation, to ensure all suspended sediment had settled out of the water column.
treatment: Category of core (none indicates regular sediment cores, none_bw indicates cores that were incubated and contained only bottom water, bottom_water indicates oxygen measures taken from the store of replacement bottom water collected at each site.
core_ID:
Unique identifier for each sediment core. Follows the format: <site code><sampling date id>_<enrichment type>_<replicate id>.
water_height: Height (cm) of the water overlying sediment in each core. Used to calculate total volume.
oxy_0h: Oxygen concentration (μmol/L) in each core, measured at the beginning of the incubation (hour zero).
oxy_4h: Oxygen concentration (μmol/L) in each core, measured 4 hours into the incubation.
oxy_8h: Oxygen concentration (μmol/L) in each core, measured 8 hours into the incubation.
oxy_12h: Oxygen concentration (μmol/L) in each core, measured 12 hours into the incubation.
oxy_16h: Oxygen concentration (μmol/L) in each core, measured 16 hours into the incubation.
oxy_20h: Oxygen concentration (μmol/L) in each core, measured 20 hours into the incubation.
oxy_24h: Oxygen concentration (μmol/L) in each core, measured 24 hours into the incubation.
Contains duplicate measurements of nutrient concentrations (mmol m -3) for each core_id, at each time point in the incubation (start: 0h, mid: 12h, end: 24h). Nutrient concentration columns follow the format: <replicatenumber>_<nutrient>.
Code
000-study-site-map.R
Generates components of the study site map (Fig. 1) in Clinton et al.
- Figure outputs: study_sites_map.tiff, NorthAmerica_inset_map.tiff
001-nutrients-and-oxygen-data-setup.R
Script used to clean nutrient and oxygen data, and join them into one data frame.
- Input:
IncubationData_UnenrichedCores.csv,NutrientData_UnenrichedCores.csv - Processed data outputs:
Oxygen_and_Nutrient_concentrations_raw.csv
002-nutrients-and-oxygen-apply-corrections.R
Applies corrections to nutrient and oxygen values, required to account for the nutrients and oxygen added via bottom water additions during incubations. Then, calculates flux rates in units of mmol/m2/day and umol/m2/day for each of the 4 nutrients and for oxygen.
- Input:
Oxygen_and_Nutrient_concentrations_raw.csv - Processed data outputs:
Oxygen_and_Nutrient_concentrations_corrected.csv,Oxygen_and_Nutrient_FluxRates_AllCores.csv,
Mean_Oxygen_and_Nutrient_FluxRates_AllCores.csv
003-nutrients-and-oxygen-plots.R
Creates a plot of oxygen consumption in each core and bar plots of nutrient flux rates.
- Input:
Oxygen_and_Nutrient_concentrations_corrected.csv,Oxygen_and_Nutrient_FluxRates_AllCores.csv,
Mean_Oxygen_and_Nutrient_FluxRates_AllCores.csv - Figure output:
OxygenConsumption_AllSites.tiff,Nutrient_FluxRates_AllSites_ByCore.tiff,Nutrient_FluxRates_AllSites.tiff
004-nutrients-and-oxygen-site-comparisons.R
Tests for differences in flux rates among sites using univariate (GLS) and multivariate (PERMANOVA, nMDS) analyses.
- Input:
Oxygen_and_Nutrient_FluxRates_AllCores.csv - Figure output:
spiderplot_nutrient_fluxes.tiff
005-community-diversity-data-setup.R
Script used to clean community composition data and join them into one data frame. Includes taxonomic and functional diversity indices, community-weighted trait means, and functional group abundances.
- Input:
functional_groups.csv,abundances_by_core.csv,Oxygen_and_Nutrient_FluxRates_AllCores.csv,
TD_indices.csv,FD_indices.csv,CWM_by_core.csv - Processed data output:
All_Explanatory_and_Response_Variables.csv
006-hill-numbers.R
Estimates the first three Hill numbers (q = 0, 1, 2) following Chao et al. (2014).
- Input:
abundances_by_core_for_calculation_of_indices.csv - Processed data output:
Hill_numbers.csv
007-diversity-and-fluxes-multivariate-analysis.R
Performs Redundancy Analyses (RDAs) on multivariate benthic fluxes as a function of diversity indices, CWMs, and Hill numbers. Variation partitioning describes the portion of variance explained by each set of explanatory variables.
- Input:
All_Explanatory_and_Response_Variables.csv,Hill_numbers.csv - Figure outputs:
RDA_Indices_ggvegan_plot.tiff,RDA_CWMs_ggvegan_plot.tiff
008-mechanistic-links-FGs-and-fluxes.R
Uses Generalized Additive Models (GAMs) to explore potential drivers of ammonium and oxygen flux. Produces supplementary plots of the relationships between NH4+ flux, O2 flux, and total macrofaunal abundance.
- Input:
All_Explanatory_and_Response_Variables.csv - Figure outputs:
Oxygen_flux_vs_macrofaunal_abundance.tiff,Ammonium_flux_vs_macrofaunal_abundance.tiff,
Ammonium_flux_vs_oxygen_uptake.tiff
References
Clinton ME, Snelgrove PVR, Bates AE (2025). Testing biodiversity-ecosystem function relations in nearshore marine sediments. Estuarine, Coastal and Shelf Science.
Clinton ME, Snelgrove PVR, Bates AE (2024). Macrofaunal diversity patterns in coastal marine sediments: Re-examining common metrics and methods. Mar Ecol Prog Ser 735:1-26.