Congruent life stage hotspots identify essential fish habitats to fill biodiversity conservation gaps

Open-science code and data to model, map, and analyze biodiversity-based Essential Fish Habitats (EFHs) across the Western Mediterranean Sea using hierarchical joint species distribution models (HMSC) and hurdle models, followed by spatio-temporal prediction, hotspot analysis, and MPA coverage assessment.

All figures can be regenerated directly from the provided outputs (once the hurdle stack is downloaded from Zenodo), without rerunning scripts 1–5.

Purpose

• Quantify species distributions and abundances with HMSC (presence/absence and abundance components).
• Generate spatio-temporal predictions and hurdle-model products for multiple life stages.
• Identify emerging hotspots/coldspots and evaluate marine protected area (MPA) coverage of EFH patterns.

Methods and Workflow

1. Model specification: Assemble inputs and define HMSC models (1-Define HMSC models.r).
2. Model fitting: Fit presence/absence and abundance models (supports HPC workflows) (2_1-Fit presence absence model.r, 2_2-Fit abundance model.r); Slurm job scripts for cluster runs are available in Models/.
3. Diagnostics: Assess convergence, effective sample sizes, and PSRF (3-Models convergence and fit.r).
4. Spatial predictions: Generate posterior samples and means over the spatial grid (4-Spatial predictions with fitted HMSC.r).
5. Hurdle outputs: Combine presence/absence and abundance to produce hurdle predictions and derived rasters (5-Hurdle model predictions.r).
6. Hotspot analysis: Emerging Hotspot Analysis (EHSA) and visualization scripts (6-Emerging Hotspot Analysis.r, Figure_3_4.r, Figure_5_6.r).
7. MPA coverage: Quantify overlap of hotspot classes with marine protected areas (Figure_7.r).

Inputs and Data

• Primary inputs: RData objects and spatial layers in Inputs_HMSC/.
  • Species/community matrices, traits, and covariates (e.g., allData.RData, Traits.RData).
  • Prediction grids and protection layers (shapefiles in Grid/ and protectionPLactivities_niv1&2&3&4/).
• Model checkpoints: Fitted models and MCMC outputs in Models/.
• Hurdle prediction stack: Spatio-temporal raster stack available on Zenodo . Use download_zenodo_spatiotemporal.R to place r_stack_Hurdle_0.05_0_1000_19992021.Rdata into Outputs/Spatio_temporal_prediction/; download this file before running scripts that depend on the hurdle stack (e.g., 6-Emerging Hotspot Analysis.r, Figure_2.r, Figure_3_4.r, Figure_5_6.r).

Outputs

• Posterior samples and means for presence/absence and abundance: Outputs/Presence_Absence_posteriors/, Outputs/Abundance_posteriors/, Outputs/Presence_Absence_mean_posterior/, Outputs/Abundance_mean_posterior/. Full posterior outputs are several hundred GB and not hosted on Zenodo; they can be requested from [AUTHOR_NAME].
• Hurdle predictions (combined): Outputs/Hurdle_prediction/. These combined hurdle outputs are large and not mirrored on Zenodo; contact [AUTHOR_NAME] if you need access.
• Spatio-temporal hurdle stack: Outputs/Spatio_temporal_prediction/. Source data are hosted on Zenodo (see the hurdle prediction stack entry above).
• EHSA results: Outputs/EHSA/.
• Figures and supplementary plots: Figures/.

Installation and Environment

• R version: Tested with R (>= 4.2.x).
• Core packages: Listed in 0-Load libraries.r; major dependencies include Hmsc, raster, sf, ggplot2, dplyr, zenodor, and others.
• Install dependencies: The loader script installs missing CRAN packages; zenodor is installed from GitHub. From the repo root:

source("0-Load libraries.r")

If zenodor is absent, it is pulled via remotes::install_github("FRBCesab/zenodor") inside download_zenodo_spatiotemporal.R.

Basic Usage

1. Load packages and data paths:

   source("0-Load libraries.r")

2. Define and prepare models:

   source("1-Define HMSC models.r")

3. Fit models (presence/absence and abundance) — adjust settings for local vs HPC:

   source("2_1-Fit presence absence model.r")
   source("2_2-Fit abundance model.r")

4. Run diagnostics and convergence checks:

   source("3-Models convergence and fit.r")

5. Produce spatial predictions and hurdle outputs:

   source("4-Spatial predictions with fitted HMSC.r")
   source("5-Hurdle model predictions.r")

6. Retrieve the published hurdle stack (optional if not already present):

   source("download_zenodo_spatiotemporal.R")

7. Run Emerging Hotspot Analysis (EHSA):

   source("6-Emerging Hotspot Analysis.r")

8. Generate figures as needed (e.g., Figure_1.r, Figure_2.r, Figure_3_4.r, Figure_5_6.r, Figure_7.r or Figure_S*.r).

Repository Structure

• 0-Load libraries.r: Dependency management and loading.
• 1-Define HMSC models.r: Model setup and data preparation.
• 2_1-Fit presence absence model.r, 2_2-Fit abundance model.r: Model fitting workflows.
• 3-Models convergence and fit.r: Diagnostics.
• 4-Spatial predictions with fitted HMSC.r: Spatial posterior predictions.
• 5-Hurdle model predictions.r: Combine components into hurdle outputs.
• 6-Emerging Hotspot Analysis.r: Emerging Hotspot Analysis workflow.
• Figure-generation scripts: all Figure_*.r scripts in the repository root (including main and supplementary figures; e.g., Figure_7.r, Figure_S8_S9.r).
• 📁 Inputs_HMSC: Input data and spatial layers.
• 📁 Models: Fitted model objects and initialization files.
• 📁 Outputs: Generated predictions, EHSA results, and hurdle outputs (git-ignored intermediates where appropriate).

Citation

If you use this code or data, please cite:

• Dataset: [AUTHOR_NAME] (2026). Community-based Essential Fish Habitats in Mediterranean Sea (Version V1.0.0). Zenodo.
• Code: This GitHub repository (XXXXX/Biodiversity_based_EFHs, main branch).
• Associated publication: XXXXXXXX.

License

This project is released under the GNU General Public License v3.0. See LICENSE for details.