Dataset: Impacts of submarine groundwater discharge on benthic community composition and functional diversity on coral reefs in Mo'orea, French Polynesia in 2023

Detailed methods are outlined in the results publication Barnas et al. (2025) and summarized here. 

Study site and characterization

Mo‘orea, French Polynesia, is a tropical volcanic island with coastal fringing coral reefs where SGD is distributed through fissures in the reef plate (Knee et al. 2016; Hagedorn et al. 2020). Local fishers’ knowledge of an SGD seep informed the location of our survey site, and the presence of SGD was confirmed through spatial and temporal radon (Hagedorn et al. 2020, 2024) and biogeochemical surveys (Silbiger et al. 2023). We identified a focal seepage point on the western shore of Mo‘orea and haphazardly chose 19 survey locations downstream of the seep to study the effects of SGD on taxonomic and functional diversity. All survey locations had hard substrate with an average depth of 0.6 m and were within 150 m of the SGD seep, experiencing a gradient of SGD influence. Our field site experiences consistent northwestward unidirectional flow averaging 0.15 m/s (Silbiger et al. 2023), distributing SGD in a predictable alongshore gradient.

Biogeochemical measurements associated with SGD influence were assessed through discrete water sampling from high and low tides during the day and nighttime in August 2021 (n = 4 measurements per survey location). See Silbiger et al. (2023) for detailed methods and descriptions of the SGD gradient. In brief, water samples were collected concurrently at each time point in acid-washed, triple-rinsed 1 L HDPE bottles. Salinity, temperature, and pH were immediately measured using portable sensors (salinity accuracy ± 1.0% psu and precision = 0.1 psu, temperature accuracy ± 0.3 °C and precision = 0.1 °C, YSI Pro2030, Xylem Inc., Washington D.C, USA.; pH [total scale] accuracy ± 0.002 and precision = 0.001, tris-calibrated ROSSTM double junction electrode, Orion Star A325, Thermo Fisher Scientific Inc., Waltham, MA, USA). The water samples were also filtered through a 0.22 μm Sterivex filter before being frozen at −20 °C for subsequent nutrient analysis for concentrations of silicate [SiO₃²⁻], phosphate [PO₄³⁻], and nitrate + nitrite [N+N]). The samples were brought to the S-LAB at the University of Hawai‘i, where they were analyzed on a Seal Analytical AA3 HR Nutrient Analyzer (N+N: detection limit [DL] = 0.009 and coefficient of variation [CV] = 0.3%; PO₄³⁻: DL = 0.011 and CV = 0.2%; SiO₃²⁻: DL = 0.03 and CV = 0.5%). We calculated the coefficient of variation (CV = 100 × standard deviation/mean) for each biogeochemical parameter to characterize the SGD gradient for this study. CV was selected because sites most affected by SGD experienced both more extreme mean values and higher variability as SGD is pulsed onto the reef in association with the tidal cycle—SGD fluxes are highest during low tide (Burnett et al. 2006).

Community surveys

Benthic communities were surveyed via snorkeling at each survey location and at the SGD seepage point in June–July 2022. Our survey methods captured the species composition of coral, macroalgae, sponges, corallimorphs, anemones, and cyanobacteria. Composition was assessed within 2 × 2 m plots using a uniform point-count method with 200 evenly distributed points. Organisms at each point were identified to the species level when possible, or to the lowest possible taxonomic unit (Payri et al. 2000; Bosserelle 2014). Of the 51 taxa identified in this study, only six of those taxa could not be identified to the species level. In these cases, broader taxonomic classifications were necessary when identifying organisms in the community (i.e., ‘Crustose Corallines’ [CCA], Cyanobacteria unknown, Porifera unknown, Dictyosphaeria sp., Verongida sp., and turf). Therefore, we use the term ‘taxa’ instead of ‘species’ for accuracy in this dataset. Importantly, given our understanding of the life history of these broader groups, the use of these broad taxonomic classifications did not hinder our trait-based identifications. Taxa unidentifiable in the field were photographed and fragmented or collected whole for later identification. Substrate types were also identified at each point as sand, rubble, dead coral, or live coral to give context for taxon presence and abundance at each location.

Rugosity (an in situ measurement of structural heterogeneity) was measured by laying a 2 m length chain (15 mm link size) over the benthos at three parallel locations within the survey area at each location. We then calculated the ratio of the transect length of the draped chain to the total linear chain length for each measurement (Risk 1972). Mean rugosity was calculated by the average of these three ratios and subtracted from one, such that higher values reflect greater structural heterogeneity. We use the term ‘structural complexity’ as a synonym for ‘rugosity’ throughout for ease of interpretation.

Classification of functional traits

Each identified taxon was categorized into functional groups, which were selected for their contribution to broader community ecosystem functioning: phyla, morphology, calcification type, and trophic group (McGill et al. 2006). The combination of these functional groups comprises each taxon’s functional entity (FE), which provides context for each taxon’s ecological role within its community (Villéger et al. 2011; Chao et al. 2014). For example, the morphology of stony corals has been linked to photosynthetic and calcification efficiency, such that weedy branching corals exhibit greater rates of calcification than digitate or encrusting species (Alvarez-Filip et al. 2013). Conversely, branching and encrusting corals with minimal self-shading exhibit higher rates of photosynthesis and respiration than dense digitate species with self-shading and reduced interstitial flow (Carlot et al. 2022; Gattuso et al. 1999; Dennison and Barnes 1988). Relative growth rates among scleractinian corals are also dependent on morphology, such that tabular and branching species exhibit faster growth compared to those with massive morphologies (Zawada et al. 2019; Madin et al. 2020). Calcification functional traits provide insights into rates of calcification as well as to the resilience of calcifiers under environmental stress (Pentecost 1991). The phyla and functional traits specified within each functional group encompass the possible phyla and traits available from the full surveyed community taxon pool. Functional identification of each taxon was accomplished using the World Register of Marine Species (WoRMS), CoralTrait Database, AlgaeTraits, species-identification guides, and primary literature. Notably, we were unable to identify all organisms to the species level. However, the functional entities ascribed to these broader classifications were consistent with characteristics of these taxa, both in the literature and according to our observations.

Taxonomic and functional diversity

We took a multi-framework approach to identifying taxonomic and functional diversity, using a combination of raw data, multidimensional space, and dissimilarity-based methods (Mammola et al. 2021). We calculated three diversity metrics to measure community shifts along the SGD gradient: proportional taxon richness (raw data), functional entity richness (raw data), and volume of functional entity trait space (multidimensional space). We also used Gower’s distance metric and Bray−Curtis dissimilarity matrices to characterize functional dispersion and taxa dissimilarity, respectively, as described in the statistical analyses section below (dissimilarity-based method). Taxon richness was determined as the total number of unique species or taxonomic units within each survey plot. Similarly, each taxon was represented by one functional entity (FE), where each FE encompassed the unique combination of functional traits from all functional groups—phyla, morphology, calcification type, and trophic group (Villéger et al. 2011). FE richness was determined by the total number of unique FEs within each survey plot. Relative taxon richness and FE richness were calculated as the total number of unique taxa or FEs present within each survey location relative to the total number of taxa (Taxon richness_T) or total number of functional entities (FE richness_T) observed across the full community, as follows:

% Taxon richness = 100 × (Taxon richness ÷ Taxon richness_T)
% FE richness = 100 × (FE richness ÷ FE richness_T)

The number of functional entities present at each site may have been equal to or less than the total number of taxa, and FE richness < taxon richness indicates functional redundancy, where more than one taxon shared the same functional entity and occupied a similar functional role in the community (Yachi and Loreau 1999). Functional entity volume, described as the volume of FE in multidimensional trait space, represents the dispersion of functional entities in multidimensional space through FE dissimilarity (Teixidó et al. 2018; Villéger et al. 2011). High FE volume indicates greater richness and dissimilarity across functional entities in a given surveyed community and therefore a wider range of functional roles, with less overlap in functionality. To calculate FE volume, a dissimilarity matrix of each survey location was calculated for FE using the daisy function with Gower’s distance metric (de Bello et al. 2013) in the cluster package in R, version 2.1.3., accessed 24 June 2023 (Teixidó et al. 2018; Maechler et al. 1999). Volumes of each survey site were calculated using the convhulln function in the geometry package in R, version 0.4.7, accessed on 24 June 2023 (Roussel et al. 2005).

Statistical analyses

We used multiple statistical approaches to test the effect of SGD on taxa and functional richness as well as community composition. We employed a regression approach to assess continuous changes in environment and communities along the SGD gradient. Indeed, recent reviews highlight the power of using regression-based experimental designs, which better characterize mechanisms compared to ANOVA designs (Idjadi and Edmunds 2006). We used individual general linear models (GLM) to determine the effect of SGD on the suite of functional and taxonomic diversity metrics while controlling for structural complexity, which could impact benthic taxonomic diversity by affecting settlement substrate (Idjadi and Edmunds 2006). To test the effect of structural complexity on %Taxon richness, %FE richness, and %FE volume in trait space, we used GLMs with mean structural complexity as the independent variable. We then calculated residuals of each diversity metric as a function of structural complexity. These residuals were used to test the impact of SGD on diversity above and beyond the effect of structural complexity. Due to the overall dominance of stony coral and fleshy macroalgae within the study site, as well as the ecological relevance of these functional groups to overall ecosystem health within a coral reef (Hatcher 1990; Hoegh-Guldberg et al. 2007), we additionally assessed the taxonomic and functional diversity of coral and fleshy macroalgae separately along the SGD gradient. All taxa used for the coral and macroalgae analyses were identified to the species level. Because there are several biogeochemical metrics commonly associated with SGD (i.e., variability in salinity, temperature, pH, and nutrients) (Taniguchi et al. 2019), we used a model-selection approach to determine the dominant SGD-related physicochemical variables and possible interactive effect of structure; selection involved comparing the AIC_C (Akaike information criterion, corrected for small sample size) of regression models. We tested both linear and polynomial regressions because communities exposed to various intensities of SGD may exhibit different relationships with diversity along the gradient in response to distinct biogeochemistry at each location.

We assessed functional-trait dispersion across surveyed species in multidimensional functional space using a principal coordinate analysis (PCoA) with the Gower metric. The functional space was created by calculating pairwise distances between taxa for four functional groups. To test the effect of SGD on community composition along the gradient, we used generalized additive models (GAM) to fit nonlinear relationships to the full suite of SGD parameters on community composition. Taxa and FE composition dissimilarities were visualized through an nMDS with a Bray−Curtis dissimilarity index, and we used the ordisurf function in the vegan package, versions 2.4.0-2.4.2 (Oksanen et al. 2003) to create a smooth fit of each parameter in ordination space. All analyses were completed in R version 4.3.2 (R Core Team 2023), and all visuals were produced with ggplot2 (Wickham 2016).