DOI: 10.5061/dryad.w0vt4b93n
title: Bivalve mortality and burrowing behavior under compound stress of salinity and temperature: Data from a mesocosm experiment

description: The frequency of extreme weather events, such as heatwaves and flash
 floods, is increasing due to global climate change. While many studies
 have investigated the potential impact of high temperature and salinity
 fluctuations on benthic animals in tidal flats, few have explored how
 these stressors affect the relationship between native and introduced
 species, including the potential establishment of introduced species in
 response to extreme short-term events. This study investigates the
 compound effects of heatwaves and salinity changes on the survival and
 burrowing behaviors of native and introduced bivalve species.
 Specifically, the introduced species Ruditapes philippinarum and the
 native species Cerastoderma edule were cultivated in four salinity
 settings under five temperature blocks using a mesocosm experiment. The
 results show that R. philippinarum exhibits lower mortality than C. edule
 under temperatures exceeding 35 °C, especially those in lower salinity
 levels of 10 and 20 PSU. Interestingly, the burrowing behaviors of the two
 species did not differ significantly, indicating that they share similar
 ecological functions under future climate change. However, if heatwaves
 coincide with other stressors, like strong freshwater inputs, as shown in
 this study, a potential widening of the establishment window for
 introduced species may occur in future climate change scenarios. These
 findings provide valuable insights into how compound extreme weather
 events may drive species shift on tidal flats, with important implications
 for understanding population dynamics under future climate change.
description: # Bivalve mortality and burrowing behavior under compound stress of
 salinity and temperature: Data from a mesocosm experiment
 [https://doi.org/10.5061/dryad.w0vt4b93n](https://doi.org/10.5061/dryad.w0vt4b93n) ## Description of the data and file structure This repository archives the data and code associated with the manuscript titled, “Compound extreme events reshuffle the stacked odds in the gamble between native and introduced bivalves”. The materials are organized into two main components to support transparency and reproducibility: * **dataScript_figsTbs.zip**: A compressed folder containing the raw data and R scripts used to generate the figures and tables in the manuscript. * **neutral_process_model.R**: A standalone R script that runs the neutral process model described in the paper. Below is a breakdown of the folder contents and subdirectory structure. ### Files and directories #### File: dataScript\_figsTbs.zip This compressed archive contains experimental data and R scripts used to generate the figures and tables in the manuscript. The folder structure is organized to mirror the order of figures and tables for easy navigation. ##### Folder Structure Overview: ``` dataScript_figsTbs/ ├── README.txt ├── Fg3_survivalPercnt/ ├── Fg4_movt_intensity/ ├── Fg5_establish_window/ ├── Tb1_stats_survivalPercnt/ └── Tb2_stats_movtIntensity/ ``` ``` dataScript_figsTbs/ ``` This is the main directory containing all subdirectories related to data processing, analysis, and visualization for the manuscript. It includes: ``` ./README.txt ``` You are currently reading this file. It provides an overview and documentation for all files and folders included in the archive. ``` ./Fg3_survivalPercnt/ ``` This folder contains data and the R script used to generate Figure 3, which illustrates time-series variation in bivalve survival percentages during the experiment. ##### Contents: ``` ../sr_2sp.csv ``` A CSV file with survival data for two bivalve species under different salinity and temperature treatments. Columns include: * **date**: Date of data collection (YYYYMMDD format) * **exp_day**: Experimental day (elapsed days) * **amount**: Number of dead individuals * **species**: Species code (“RP” = *Ruditapes philippinarum*, “CE” = *Cerastoderma edule*) * **salinity**: Salinity level (5, 10, 20, or 30 PSU) * **treatment**: Temperature treatment (“C” = control, “H” = heated) * **room**: Climate room ID (Room 310 or 319) * **sediment**: Sediment type (“muddy” = D50: 120.06 µm; “sandy” = D50: 265.02 µm) * **csum**: Cumulative number of dead individuals - **sr**: Survival rate calculated as: sr = (total_number - csum) / total_number, where total_number = 12 per treatment combination. ``` ../survivalPercent_Fig3.R ``` This R script processes the survival data and generates a time-series plot. It calculates daily survival percentages across four salinity levels and two temperature treatments (control and heatwave). Plot features include: * facets (panels): One for each salinity level (in PSU) * x-axis: Experimental days grouped in ~5-day intervals * y-axis: Mortality percentages * red lines: Heated (heatwave) treatment * blue lines: Control (ambient temperature) * solid lines: Native species (*C. edule*) * dashed lines: Introduced species (*R. philippinarum*) ``` ./Fg4_movt_intensity/ ``` This folder contains data and an R script related to the analysis and visualization of bivalve movement intensity under varying salinity and temperature conditions, corresponding to Figure 4 of the manuscript. ##### Contents: ``` ../bdepth_mortality.csv ``` This file includes measurements of daily burrowing depth changes and associated metadata. Columns are as follows: * **sp_code**: Species code (“rp” = *R. philippinarum*, “ce” = *Cerastoderma edule*) * **date**: Date of measurement (DD/MM/YYYY format) * **sediment**: Sediment type — “muddy” (D50 = 120.06 µm) or “sandy” (D50 = 265.02 µm) * **treatment**: Heating treatment — “control” (ambient) or “heating” (heatwave simulated) * **salinity**: Salinity level — “S5”, “S10”, “S20”, “S30” (5, 10, 20, 30 PSU) * **b.sd**: Standard deviation of daily burrowing depth (in mm) * **max.change**: Maximum daily burrowing depth change (in mm) * **temp.m**: Mean low tide temperature in the mesocosm tank (°C) * **match_label**: Unique ID combining species, date, sediment, treatment, and salinity * **mr**: Mortality rate of the individual ``` ../movt_intensity_Fig4.R ``` This R script processes the burrowing depth data and generates a plot depicting changes in movement intensity over time. Daily burrowing depth was tracked using a thread-based method, and movement variability was quantified using the standard deviation of daily depth change (b.sd). To assess movement intensity over time, a time-series linear regression model was applied: log(Depth_variance) = log(a) + log(b) × Time, where the slope b served as the movement intensity index. To evaluate how salinity and temperature influenced this index, a polynomial regression model was used: lm(MovementIntensity ~ poly(Salinity, 2, raw = TRUE) × Temperature), with salinity treated as a four-level factor (5, 10, 20, and 30 PSU) and temperature categorized as either above or below 25 °C. Plot features are as below:  * **Solid lines**: Native species (*C. edule*) * **Dashed lines**: Introduced species (*R. philippinarum*) * Each data point reflects the slope from a time-series regression (movement intensity index). * No error bars are shown, as values are derived from individual-level slopes (n = 2,542 observations). ``` ./Fg5_establish_window/ ``` This folder contains survival data and an R script used to generate a conceptual figure (Figure 5) illustrating how compound stress from salinity and temperature may create a disturbance-driven “establishing window” that favors introduced species over native ones. ##### Contents: ``` ../sr_2sp.csv ``` This data file contains survival percentages for two bivalve species across various experimental treatments. Columns include: * **date**: Data collection date (YYYYMMDD format) * **exp_day**: Number of days since the experiment began * **amount**: Number of dead individuals recorded * **species**: Species identifier (“RP” = *R. philippinarum,* “CE” = *C. edule*) * **salinity**: Salinity level (5, 10, 20, or 30 PSU) * **treatment**: Temperature treatment (“C” = control, “H” = heated) * **room**: Climate room ID (Room 310 or 319) * **sediment**: Sediment type — “muddy” (D50 = 120.06 µm) or “sandy” (D50 = 265.02 µm) * **csum**: Cumulative number of dead individuals * **sr**: Survival rate calculated as: sr = (total_number - csum) / total_number, with total_number = 12 per salinity × treatment combination ``` ../establish_window_Fig5.R ``` This R script processes the survival data to create a conceptual figure illustrating how a disturbance-driven “establishing window” may widen the mortality gap between native and introduced bivalves under compound stress. To estimate cumulative thermal stress, the script calculates the daily mean temperature difference between heated and ambient (control) conditions during low tide periods, then computes cumulative sums in five-day intervals (aligned with each heatwave block). Establishing conditions for the introduced species (*R. philippinarum*) are inferred by comparing survival percentages between species under different salinity and cumulative temperature conditions, using the difference: Survival Difference = Survival_*R. philippinarum* – Survival_*C. edule* ``` ./Tb1_stats_survivalPercnt/ ``` This folder contains data and an R script used to generate Table 1, which summarizes the results of linear regression models assessing survival percentages of two bivalve species under heatwave conditions across varying salinity levels. ##### Contents: ``` ../bdepth_mortality.csv ``` This dataset includes burrowing depth measurements and associated environmental variables used for survival analysis. Column descriptions: * **sp_code**: Species code (“rp” = *R. philippinarum*, “ce” = *C. edule*) * **date**: Data recording date (DD/MM/YYYY format) * **sediment**: Sediment type — “muddy” (D50 = 120.06 µm) or “sandy” (D50 = 265.02 µm) * **treatment**: Temperature treatment — “control” (ambient) or “heating” (heatwave simulation) * **salinity**: Salinity level — “S5”, “S10”, “S20”, or “S30” (5–30 PSU) * **bDepth.sd**: Standard deviation of burrowing depth change (mm) * **max.change**: Maximum daily burrowing depth change (mm) * **temp.mean**: Mean simulated low-tide temperature (°C) * **mortalityRate**: Mortality rate of individual bivalves ``` ../stats_survivalPercnt_Table1.R ``` This R script fits linear regression models to examine how survival percentages are influenced by experimental duration, temperature treatment (control vs. heatwave), and species identity across different salinity levels. The model used is: log(Survival Percentage) ~ Experimental Duration × Heatwave Treatment × Species where the interaction term Experimental Duration × Heatwave Treatment captures the cumulative impact of increasing thermal stress over time. Model results are summarized by salinity treatment and species, with statistical significance indicated as: * for *p* &lt; 0.05 and ** for *p* &lt; 0.01.  ``` ./Tb2_stats_movtIntensity/ ``` This folder contains the data and R script used to generate Table 2, which summarizes a polynomial regression analysis of movement intensity indices across different salinity levels and heatwave treatments. The same raw data used for Figure 4 is used here. ##### Contents: ``` ../bdepth_mortality.csv ``` This dataset includes measurements of burrowing behavior and experimental conditions for individual bivalves. Column descriptions: * **sp_code**: Species identifier (“rp” = *R. philippinarum*, “ce” = *C. edule*) * **date**: Date of data collection (DD/MM/YYYY format) * **sediment**: Sediment type — “muddy” (D50 = 120.06 µm) or “sandy” (D50 = 265.02 µm) * **treatment**: Temperature treatment — “control” (ambient) or “heating” (heatwave simulation) * **salinity**: Salinity level — “S5”, “S10”, “S20”, “S30” (5–30 PSU) * **b.sd**: Standard deviation of daily burrowing depth change (mm) * **max.change**: Maximum daily depth change (mm) * **temp.m**: Mean low tide temperature in mesocosm tank (°C) * **match_label**: Unique identifier for each individual (composite of species, date, sediment, treatment, salinity) * **mr**: Mortality rate ``` ../stats_movtIntensity_Table2.R ``` This R script fits a polynomial regression model to evaluate how bivalve movement intensity responds to varying salinity and temperature conditions, using the formula: lm(Movement Intensity ~ poly(Salinity, 2, raw = TRUE) × Temperature). Movement intensity is quantified as the slope of time-series regressions representing the daily variance in burrowing depth (n = 2,542 observations). The analysis includes four salinity levels (5, 10, 20, and 30 PSU), with temperature categorized as either above or below 25 °C to distinguish heatwave effects. The model compares responses between two species: *C. edule* (native) and *R. philippinarum* (introduced). Statistical significance is indicated as * for *p* &lt; 0.05 and ** for *p* &lt; 0.01. #### File: neutral\_process\_model.R This R script implements a numerical model based on neutral process theory to simulate interactions between native and introduced species within a shared ecological niche. In this framework, species-specific competitive advantages are excluded, meaning that any shifts in species composition result solely from neutral, stochastic processes. The ecological niche is modeled as a collection of discrete sites, each of which may be occupied by a native species, an introduced species, or remain vacant. Transitions between these states occur through species-neutral processes such as natural turnover and random mortality events. At each time step, the model calculates the average state across all sites, providing a continuous measure of species composition over time. The model explores three mortality scenarios: 1. **Equal mortality**: Both species experience the same mortality rate (0.5). 2. **Increased native mortality**: Native species have a higher mortality rate (0.7) compared to introduced species (0.3). 3. **Increased introduced mortality**: Introduced species have a higher mortality rate (0.7), while native species have a lower rate (0.3). These scenarios demonstrate how even neutral dynamics, when paired with differing mortality rates, can drive shifts in species dominance. This conceptual model is self-contained and does not require any external data to run. ## Code/software All analyses in this repository were performed using **R 4.1.1** (R Core Team, 2024). Each script is thoroughly annotated to guide through the workflow. The following R packages were loaded before runing the scripts: * tidyr * zoo * ggplot2 * raster * viridis **Reference** R Core Team. (2024). *R: a language and environment for statistical computing* [Computer software]. R Foundation for Statistical Computing. [https://www.R-project.org](https://www.R-project.org) ## Access information Data was derived from the following sources: * A mesocosm experiment conducted at the Department of Estuarine and Delta Systems, NIOZ Royal Netherlands Institute for Sea Research, Yerseke, The Netherlands. Usage Restrictions: * The data are cross-referenced and linked to this study. For any further academic/commercial use of the raw data prior to the manuscripts publication, a formal request for permission must be submitted to the coauthors.