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 < 0.05 and ** for p < 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 < 0.05 and ** for p < 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

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.