No-take marine reserves promote oligotrophic reef bacterioplankton communities across the Great Barrier Reef
Authors: Marko Terzin1,2,3,5*, Steven J. Robbins4, Kim-Anh Lê Cao5, Sara C. Bell1, Katherine E. Dougan4, Julian Zaugg4, Renee K. Gruber1, Michael J. Emslie1, Daniela M. Ceccarelli1, Samuel Chaffron6,7, Philip Hugenholtz4, Nicole S. Webster1,4,8, David G. Bourne1,2,3, Yun Kit Yeoh1,3, Patrick W. Laffy1,3,*
1Australian Institute of Marine Science, PMB no3 Townsville MC, Townsville QLD 4810
2College of Science and Engineering, James Cook University, Townsville, 4811
3AIMS@JCU, James Cook University, Townsville QLD 4811
4Australian Centre for Ecogenomics, School of Chemistry and Molecular Biosciences, The University of Queensland, St Lucia, QLD 4072
5Melbourne Integrative Genomics and School of Mathematics and Statistics, University of Melbourne, Melbourne, Parkville VIC 3052
6Nantes Université, École Centrale Nantes, CNRS, LS2N, UMR 6004, F-44000 Nantes, France
7Research Federation for the Study of Global Ocean Systems Ecology and Evolution, FR2022/Tara Oceans GOSEE, F-75016 Paris, France
8Institute for Marine and Antarctic Studies, University of Tasmania, TAS, 7001
*Corresponding authors
This study investigates the relationship between seawater microbiomes, environmental variables, and reef protection status (No-Take Marine Reserves vs. fished reefs) across the Great Barrier Reef (GBR). We integrated 876 prokaryotic metagenome-assembled genomes (pMAGs), physico-chemical water data, benthic cover, and fish survey data to identify microbial indicators of reef zoning and predict environmental conditions from microbial community composition.
The conceptual framework proposed in this study linking reef protection status, fish biomass, benthic cover, nutrient dynamics, and seawater microbial community structure.
The pipeline is entirely script-based and interpreted (R, Python, Julia); no compiled standalone binary is required. All software dependencies and version numbers are listed below, organised by analysis stage, and full session details are provided in sessionInfo.txt.
- Tested on: R 4.3.2 and Python 3.12
- Operating systems: the analysis scripts are platform-independent and have been run on Linux and macOS
- Hardware: no non-standard hardware is required; all figure scripts run on a standard desktop or laptop
See the Key R Packages and Versions and External Software & Tools tables below for the complete dependency list with versions.
git clone https://github.com/mterzin/fishy_microbes.git
cd fishy_microbesR packages install from CRAN and Bioconductor using the standard installers:
install.packages(c("vegan", "randomForest", "tidyverse", "ggplot2",
"patchwork", "igraph", "glmmTMB", "DHARMa"))
# Bioconductor packages
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(c("mixOmics", "ALDEx2", "phyloseq", "microbiome"))# Create conda environment (optional)
conda create -n modularity_env python=3.12
conda activate modularity_env
# Install required packages
pip install networkx pandas numpy scipy matplotlibTypical install time: on a normal desktop computer this may take up to ~2 hours, dominated by compiling the R/Bioconductor dependencies.
Each figure script (Figure_1 to Figure_6) is self-contained: it reads the processed input tables linked from this README and reproduces the corresponding manuscript figure and its underlying statistics. The processed input tables and the underlying primary data are openly available (see Data Availability).
Expected output: running the scripts reproduces the published main figures and their associated statistical results, e.g. the ~71% MINT sPLS-DA classification accuracy, the PERMANOVA/dbRDA variance partitioning, and the co-occurrence network metrics.
Alongside each .Rmd file we provide a rendered HTML report for all six scripts, so that anyone replicating the pipeline can see exactly what the expected output looks like at each step without re-running the code.
Expected run time: each figure script runs in minutes on a standard desktop; the permutation-based analyses in Figure_2 (999 permutations for MINT sPLS-DA, Random Forest, and ALDEx2 validation) are the most time-consuming and may take longer depending on hardware.
To reproduce a given figure and its associated analyses, open the corresponding R Markdown script in RStudio (or knit it from the command line) and run it against the provided input tables:
rmarkdown::render("Figure_2_Nature_Communications.Rmd")Each script maps directly to a manuscript figure and the analyses described in the Materials and Methods (see Analysis Workflow below). The code is released under the MIT License, permitting reuse and adaptation on the user's own data.
| Data type | Repository | Accession/Link |
|---|---|---|
| Metagenomic Sequences | EBI BioProject | PRJEB82623 |
| pMAGs (5,283 high-quality) | Zenodo | 10.5281/zenodo.17109887 |
| Processed abundance tables | Zenodo | 10.5281/zenodo.17109887 |
| Physico-chemical variables | IMOS-AODN | 10.25845/Q4XH-YN10 |
| Benthic cover & fish data | AIMS LTMP | AIMS Data Portal |
| Assembly & binning code | bioRxiv | Robbins et al. 2025 |
Metagenomes were assembled using the Aviary v0.3.3 pipeline, which generated:
- Hybrid assemblies (Illumina + Nanopore) for 27 sites
- Short-read-only assemblies for 21 sites
- Long-read processing: Guppy v5.0.16 (superaccuracy basecalling), Porechop (adapter/barcode trimming)
- Assembly: Aviary v0.3.3 with metaFlye (long-read assembly), racon and pilon (polishing), and metaSPAdes (short-read assembly)
- Binning: MetaBAT1, MetaBAT2, MaxBin2, CONCOCT, VAMB, and Rosella
- Refinement: DAS Tool v1.1.2
- Quality assessment: CheckM v1.2.2 and CheckM2 v1.0.2 (pMAGs retained at a quality threshold of ≥50, calculated as completeness − 3×contamination)
- Dereplication: CoverM v0.6 at 95% Average Nucleotide Identity (ANI)
- Result: 5,283 pMAGs → 876 "species-resolved" pMAGs95%ANI
- Read mapping: minimap2 v2.18 (via CoverM) for abundance estimation
- Taxonomy: Genome Taxonomy Database Toolkit (GTDB-Tk, release R214)
- Functional annotation: anvi'o v8 (with Prodigal v2.6.3 and HMMER v3.3.2) with KEGG Orthology (KO) database
- Metabolic pathways: KEGG module completeness (358 modules detected)
Full metagenomic processing methods (assembly, binning, taxonomy, abundance) are described in Robbins et al. 2025.
- Map of 48 offshore reefs across 7 GBR sectors
- Sampling timeline: 4 transects (Nov 2019 – Jul 2020) and 7 sectors
- Reef protection status and GBR zoning categories
This script implements a comprehensive analytical framework to identify and validate microbial indicators of reef protection status: 1. Community-level assessment:
- PCA: Principal Components Analysis to explore major sources of variation in microbial community composition (
mixOmicsv6.26.0) - PERMANOVA: Permutational Multivariate Analysis of Variance testing the effect of reef protection status on microbial communities while accounting for sampling trip, geographic sector, and reef identity as spatiotemporal covariates (
veganv2.6-4; 9,999 permutations) - dbRDA: Distance-based Redundancy Analysis for constrained ordination of protection effects 2. Indicator identification with spatiotemporal integration:
- MINT sPLS-DA: Multivariate INTegration Sparse Partial Least Squares Discriminant Analysis that identifies microbial indicators discriminating between NTMRs and fished reefs while accounting for sector-specific variation (
mixOmicsv6.26.0). Leave-One-Group-Out Cross-Validation (LOGOCV) with iterative training on six sectors and validation on the remaining sector was used to determine optimal number of MINT sPLS-DA components and features. - Permutation testing: Zone-label shuffling (999 permutations within sectors) to generate null distribution and calculate p-value and Cohen's d effect size. 3. Independent validation methods:
- Random Forest classification: Non-linear supervised learning with leave-one-sector-out cross-validation, permutation testing (999 iterations), feature importance assessment, and methodological concordance analysis (Jaccard similarity, Cohen's Kappa) (
randomForestv4.7-1.1) - ALDEx2 pairwise comparison: Welch's t-test and Wilcoxon rank-sum tests on CLR-transformed abundances without covariates
- ALDEx2 GLM: ANOVA-like differential expression with Generalized Linear Models incorporating spatiotemporal covariates (Open_or_Closed_to_fishing + Sampling_trip + SECTOR_N_S) using Dirichlet-multinomial modeling (128 Monte Carlo samples, FDR correction α = 0.05) (
ALDEx2v1.34.0) 4. Validation of indicator robustness: - Presence/absence analysis: Detection frequency quantification for indicator MAGs across 190 samples using raw count data to assess whether abundance differences reflect genuine ecological variation versus genome size bias
- Read-based validation: Assembly-independent verification using DIAMOND (v2.0.9) mapping against NCBI nr database, MEGAN (v6.23.0) taxonomic profiling, and MINT sPLS-DA on read-based profiles to confirm indicator selection across independent methods
- MINT sPLS: Integrating 876 MAGs with 54 environmental variables
- Biplots and clustered image maps (CIM)
- Co-occurrence networks: Connectedness and cohesion metrics (Herren & McMahon 2017)
- Network comparison: NTMR vs. fished reef networks (connectedness and cohesion)
- Regression: Genome size, GC content, and KEGG module completeness vs. network properties A custom-made Python script implements co-occurrence network analysis to compare microbial community structure (MODULARITY) between NTMR and fished reefs: 1. Network construction (FlashWeave):
- Co-occurrence networks: Constructed for each sector × zone combination (n = 14 networks) using FlashWeave86 (v0.19.2) via the
flashweaveJulia package - Data transformation: Raw count tables were CLR-transformed with pseudocounts added (1e-6)
- Network parameters: Sensitive mode (α = 0.01, max_k = 2) with a minimum of 4 samples per network
- Edge filtering: Only positive partial correlations (weight > 0) were retained, consistent with the predominance of positive associations in global plankton interactomes (98.5% positive edges; Chaffron et al., 2021)83 2. Modularity analysis:
- Binary network conversion: Positive correlation networks were converted to binary (unweighted) graphs by retaining edges with weight > 0
- Modularity calculation: Computed using the Clauset-Newman-Moore greedy algorithm87 as implemented in the
greedy_modularity_communitiesfunction fromNetworkXv3.6.188 - Community detection: Modularity measures the degree of network compartmentalization, where higher values indicate more compartmentalised community structure 3. Statistical comparison:
- Between-zone comparison: Mann-Whitney U tests comparing modularity values between NTMR (n = 7) and fished (n = 7) sector-specific networks
- Effect size: Cohen's d calculated to quantify magnitude of differences 4. Visualization:
- Boxplots comparing modularity distributions between protection statuses with sector-specific shapes and trip colors
- KEGG module completeness: Comparing completeness scores of 358 KEGG modules between fished reef-enriched and NTMR-enriched indicator pMAGs
- Top 45 modules: Heatmap of modules with greatest between-group differences, grouped by metabolic category (carbohydrate metabolism, energy generation, biosynthesis of cofactors, vitamins, amino acids, and lipids)
- Statistical testing: Wilcoxon rank-sum tests comparing mean completeness scores between groups
1. Random Forest environmental prediction:
- Model performance: RF models predicting continuous environmental variables, evaluated across 50 stratified permutation tests per variable (80/20 train/test split stratified by GBR sector)
- High accuracy predictions (R² > 0.6): seawater temperature (R² = 0.74), salinity, particulate nutrients (POC R² = 0.74, PN R² = 0.66), and dissolved inorganic phosphorus (PO₄³⁻ R² = 0.69, TDP R² = 0.79)
- Low-to-moderate accuracy (R² < 0.6): dissolved nitrogen species, silicate, benthic cover variables, and most fish groups (exception: corallivore fish biomass, median R² = 0.72) 2. Microbial niche inference:
- Niche modeling: Robust optimum (RO) method (Chaffron et al. 2021) applied to top 50 RF predictors per environmental variable
- Niche bounds: Q1 (lower bound), Q2 (optimum), Q3 (upper bound)
- Specialist vs. generalist taxa: Narrow niche ranges associated with high RF prediction accuracy; broad niches associated with weaker predictive power
- Example: Top 50 temperature predictors show narrow thermal niche (Q1–Q3: 27.38 ± 2.11°C to 28.38 ± 1.73°C; optimum Q2: 27.84 ± 1.88°C), predominantly Flavobacteriales (46%)
| Package | Version | Citation | Purpose |
|---|---|---|---|
| mixOmics | 6.26.0 | Rohart et al. 2017 | MINT sPLS-DA, MINT sPLS |
| vegan | 2.6-4 | Oksanen et al. 2022 | PERMANOVA, dbRDA |
| randomForest | 4.7-1.1 | Liaw & Wiener 2002 | Random Forest classification/regression |
| ALDEx2 | 1.34.0 | Fernandes et al. 2013 | Differential abundance testing |
| glmmTMB | 1.1.10 | Brooks et al. 2017 | GLMMs for environmental variables |
| DHARMa | 0.4.7 | Hartig 2022 | GLMM residual diagnostics |
| phyloseq | 1.46.0 | McMurdie & Holmes 2013 | Microbiome data handling |
| microbiome | 1.24.0 | Lahti & Shetty 2017 | CLR transformation |
| igraph | 1.5.1 | Csárdi & Nepusz 2006 | Network analysis (graph handling) |
| tidyverse | 2.0.0 | Wickham et al. 2019 | Data wrangling & visualization |
| ggplot2 | 3.5.1 | Wickham 2016 | Publication-quality graphics |
| patchwork | 1.2.0 | Pedersen 2024 | Plot composition |
| dataaimsr | 1.0.0 | Australian Institute of Marine Science | Spatial data access |
| gisaimsr | 1.0.0 | Australian Institute of Marine Science | Spatial data access |
| Software | Version | Citation | Purpose |
|---|---|---|---|
| Guppy | 5.0.16 | Oxford Nanopore Technologies | Nanopore basecalling (superaccuracy) |
| Porechop | (see repo) | Wick et al. 2017 | Adapter/barcode trimming |
| Aviary | 0.3.3 | Robbins et al. 2025 | Metagenomic assembly & binning |
| metaFlye | (see repo) | Kolmogorov et al. 2020 | Long-read metagenome assembly |
| metaSPAdes | (see repo) | Nurk et al. 2017 | Short-read metagenome assembly |
| DAS Tool | 1.1.2 | Sieber et al. 2018 | Bin refinement |
| CheckM | 1.2.2 | Parks et al. 2015 | Genome quality assessment |
| CheckM2 | 1.0.2 | Chklovski et al. 2023 | Genome quality assessment |
| GTDB-Tk | R214 | Chaumeil et al. 2019 | Taxonomic classification |
| CoverM | 0.6 | Aroney et al. 2025 | MAG dereplication & read mapping |
| minimap2 | 2.18 | Li 2018 | Read mapping (via CoverM) |
| anvi'o | 8 | Eren et al. 2015 | Functional annotation (KEGG Orthology) |
| FlashWeave | 0.19.2 | Tackmann et al. 2019 | Co-occurrence network inference |
| Python | 3.12 | Python Core Team | Modularity analysis scripting |
| networkx | 3.6.1 | Hagberg et al. 2008 | Co-occurrence network analysis & modularity (Clauset-Newman-Moore) |
| scipy | 1.12.0 | Virtanen et al. 2020 | Statistical tests in Python |
| DIAMOND | 2.0.9 | Buchfink et al. 2015 | Read-based validation |
| MEGAN | 6.23.0 | Huson et al. 2016 | Taxonomic profiling |
| Inkscape | 0.92.5 | Inkscape Project | Figure compilation |
Full session details available in sessionInfo.txt.
The modularity analysis requires a Python environment with the following dependencies:
# Create conda environment (optional)
conda create -n modularity_env python=3.12
conda activate modularity_env
# Install required packages
pip install networkx pandas numpy scipy matplotlibThis code is released under the MIT License (approved by the Open Source Initiative). The full license text is included in the repository (LICENSE), and the code may be reused and adapted on the user's own data.
If you use this code or data, please cite:
Terzin, M., Robbins, S. J., Lê Cao, K.-A., et al. No-take marine reserves promote oligotrophic reef bacterioplankton communities across the Great Barrier Reef. Nature Communications (under revision).