Study area and date
In the north-west of France, the macrotidal Bourgneuf Bay (1°-2° W, 46°-47° N; total area ~340 km2; Fig. 1) has an intertidal zone largely dominated by mudflats (exposed surface area ~100 km2). Bourgneuf Bay is situated south of the Loire estuary and is open to the sea along 12 km from the west to the north-west. C. gigas aquaculture here is of national importance and wild C. gigas reefs can account for over double the biomass of their farmed conspecifics65. Analysis of satellite observations covering 30 years of MPB biomass in the bay confirmed the co-occurrence of high MPB biomass with wild oyster reefs and cultivated stocks16 (Supplementary Methods: Wider Situation of the Reefs). Two small (each > 750 m2) wild C. gigas reefs and their immediate surroundings (10–100 m) in the north of Bourgneuf Bay were deemed suitable for experimental manipulation (yellow and orange regions in Fig. 1). Méléder et al.18 described MPB biomass as mostly concentrating around the 2 m isobath, the Falleron river channel (closest point ~400 m NNE from the eastern reef), and oyster farms. Covering this isobath, we superimposed a 350 * 350 m grid (12.25 hectares) to cover the two wild oyster reefs, orientated so that the ‘Y’ axis runs parallel to the slope of bathymetry (Fig. 1). The grid was split regularly into 49 ‘grand-cells’ of 50 * 50 m (i.e., n = 49) and each of those split into 25 cells of 10 * 10 m (i.e., n = 1225; Fig. 1). Four cells per grand-cell were chosen randomly for the sampling of meiofauna, granulometry, OM (see Table 3 for specific methodology), and macrofauna. Of these cells, only every second cell was processed for meiofauna because of time constraints in assessing their abundance.
Although there were only two oyster reef complexes (‘reefs’, hereon) in this study, multiple sampling cells fell on, or in close proximity to, each reef, so that each reef had many potential (though not independent) distance decay transects running from it capturing natural variation in spatial structure66. Comparing the ecological change following the experimental burning of oyster reefs (described below) against ecological change occurring at these two reefs over the previous 25 years16 also allowed us greater confidence to disentangle the treatment effects from typical variation. Through the centres of five grand-cells to the south of the extent, a transect forming an ‘L’ shape (Fig. 1) was sampled every 10 m for in situ MPB pigment composition and biomass. We used these data to complete the remote sensing approach for MPB biomass estimation (see below, Microphytobenthos). The western reef was slightly larger than the eastern reef and contained a large rock, ‘Roche Bonnet’, rising 0.5–1 m from the sediment. Outside the grid, another larger (200 * 80 m) wild oyster reef lies WSW at ~260 m distance from the western reef. The grid was sampled for the variables listed in Table 3 during the winter MPB low and early autumn peak seasons (see also ground-truthing in16), on the dates 18-19th September 2013 and 17-18th March 2014, before treatment, and on 7-8th October 2014 after treatment.
Microphytobenthos
We mapped MPB biomass by satellite remote sensing, following the method described in detail in Echappé et al. (2018). We used the same long-term record of high-resolution satellite images to analyse the spatial distribution of the normalised difference vegetation index (NDVI), a proxy of MPB chlorophyll a concentration at the sediment’s surface18,67, before and after treatment (individual image details in captions of Fig. 2 and Supplementary Figs. S4–S7). After atmospheric correction (FLASH and US40 aerosol model), the satellite-derived NDVI was validated against associated field measurements (r2 = 0.85, root-mean-square deviation, RMSE = 0.04, n = 57, P < 0.05; also Echappé et al., 2018). The availability of satellite images at both the desired spatial resolution (≤10 * 10 m) and meeting stringent conditions (cloud-free sky, study area emerged from tide, sun elevation >20°) was limited (Supplementary Results: Additional MPB Images). An optimal image was chosen as representative of MPB biomass patterns per season16. The study area would ideally be tidally uncovered for ~2 h before the image was taken, whereupon MPB biomass is concentrated at the sediment surface. The optimal images met this condition (i.e., Fig. 2).
To complete NDVI maps, in situ MPB pigment composition and biomass were retrieved by HPLC analysis from the 25 triplicates of sediment. These had been sampled using contact-cores to freeze the top 2 mm of sediment in situ with liquid nitrogen, with a metal surface 56 mm in diameter. Biomass was expressed by Chl a concentration (mg m−2), and dominance of MPB taxa was broadly assessed by ratio of pigment sources to Chl a: Fucoxanthin (Fuco), Diadinoxanthin (DD), Diatoxanthin (DT) and Chl c for diatoms. The ratio of unknown carotenoids (interpreted as by-products due to the low resolution of their absorption spectra) to Chl a was also analysed for ecological purposes (dominant taxa), whereas grazing pressure was investigated using the ratio of pheophorbid a to Chl a (methodological discussion in28).
Sediment variables
For laser granulometry, we sampled two depths, 0–5 cm and 5–10 cm, in triplicate at each cell. Each of the triplicate samples was put in a vial with water and sonicated. The particle size distribution was determined on a Mastersizer 3000 with a reporting range 50 nm to 3 mm. We also determined sediment percentage OM at two depths by mass loss on ignition in comparison to the oven-dried original (procedure as described for Macrofauna, also Table 3).
Macrofauna
We sampled macrofauna by a single 200 * 200 mm (depth * diameter) core per cell. Contents were placed into labelled buckets and sieved onshore (1 mm mesh). Soft-bodied polychaetes were picked out with forceps and preserved in buffered formalin during sieving. All material left on the sieve was bagged and preserved in formalin at the laboratory. Individuals were counted and measured by the longest axis (accuracy 0.1 mm, calipers); the deep-burrowing polychaete Diopatra biscayensis was counted by the presence of visible tubes above the sediment. Calibration curves from length to mass per species per season were built by identifying size classes by Sturges rule. Multiple individuals per size class (ideally n = 100) were measured to estimate mean organic mass per individual of each size class. Shell matter was physically separated from tissue, before both being dried in aluminium foil cups for 48 h at 60 °C and weighed (g) for tissue dry mass using a mass balance. Dry mass was then incinerated for four hours at 450 °C and reweighed (g; decrease in mass of the aluminium cup was also accounted for), the difference giving the organic matter mass (including residue in the shell matter), or ash free dry weight (AFDW). This number was divided by number of individuals. Calibration curves per species used first order polynomial curves for bivalves, unless numbers of size classes and individuals were small (<5 and <20 respectively), in which case the more restricted power curves were used. High variation in the mass ~size relationship, for gastropods because of issues separating shell from tissue via crushing, and for polychaetes because of their stretchable bodies, meant that power curves were again used. Calibration curves allowed the biomass of each species per sample to be estimated.
Total AFDW biomass estimates for the dominant species (≥100 individuals per season) over the study extent (0.1225 km2) were extrapolated from sample means. Biomass estimates for C. gigas were calculated in Le Bris et al.17, by a combination of remote sensing and field observations, before upper and lower range estimates were converted from wet weight to AFDW using the conversion equation for Crassostrea virginica68.
To estimate the distance-decay of predation from oyster reefs, the control reef was revisited over 24-25th July 2017. The original sampled cells were binned into distance classes from the reef i.e., 0–10, 10–20, 20–30 m, etc up to 60 m. Maintaining an even spread by distance class, 30 cells were randomly selected that were visited to deploy standardised fyke nets (nylon multifilament, mesh size ~3 mm, approx. 550 × 220 mm L × W, mouth diameter ~60 mm) over two low tides. Nets were then collected and their contents taken ashore for identification and measurement. Linear regression was used to estimate mean crab carapace width from distance from the reefs, with fishing failures (collapsed or lost nets) omitted. Mean carapace width, integrating both abundance and body size, gives an indicator of crab foraging intensity and depth, since deeper, larger infaunal prey are only available to larger individual crabs48. Our distance-decay model assumes that crabs were captured while foraging and must return to the reefs during low tide to avoid becoming prey themselves to birds4. High variance and a low sample size precludes the use of curvi-linear methods e.g., a polynomial regression term.
During the 2017 visit, we recorded reef-dwelling fauna by twenty haphazard 0.5 * 0.5 m quadrat samples over the control reef. Samples recorded the percent cover of stationary organisms, including C. gigas, counted individuals of mobile organisms, and measured carapace widths of crab species.
Meiofauna
We sampled two replicate 50 * 25 mm (depth * diameter) cores per cell for meiofauna, which were preserved in formalin onshore. In the laboratory, samples were sieved on a 40 µm mesh, re-suspended, thrice centrifuged, and extracted in a Ludox solution (density 1.31 g cm−3). The fraction was then sorted and counted by microscope down to large taxonomic groups (Nematoda, Cnidaria, Kinorhycha, Polychaeta, Oligocheta). We exclude Foraminifera and Ostracoda because they have external shells (a.k.a. tests) that accumulate in sediment, which make estimates of abundances less reliable. The huge amount of work counting meiofauna precluded the collection of any more than autumn 2013 data only,
Before-After Control-Impact (BACI) setup
To sacrifice the non-native oysters over the >750 m2 reef expanse, the most eastward of the two oyster reefs was entirely covered with a ~30 cm thick layer of straw and burned over two consecutive low tides (16-17th July 2014; termed the ‘treatment’). The wind was allowed to spread the fire rather than using any other additional substance (Supplementary Fig. S18 shows photos of the burning procedure and the effects on the oyster reef). A revisit to the reefs on the 19th July 2014 visually confirmed high-to-complete oyster mortality on the treatment reef only and that the reef shell base remained largely unmodified by the burning (e.g., Supplementary Fig. S18). The larger, westward reef remained as a control. Post-treatment sampling was 81 days afterwards to allow other natural processes that were disturbed by the burning to recover, although lingering artefacts of the burning cannot be ruled out (see discussion; also, Supplementary Results: The timescale of the impact).
Statistics and reproducibility
All our analyses and mapping techniques assumed at least second order stationarity (i.e., mean, variance and covariances constant over space). However, a bathymetric trend was expected in biological and sediment variables along the y-axis (Fig. 1), which could introduce spurious correlations, so linear trends were tested for by simple linear regressions. The residuals of significant trends were retained as the detrended variable (trends listed in Results), unless the trend was directly incorporated in the spatial prediction model (e.g., simple Kriging). Normality was also assessed for each variable by histograms and Shapiro–Wilk tests prior to parametric tests, leading to natural log or square root transformations where necessary (e.g., NDVI was LN transformed for early autumn 2013 but not for winter 2014; transformations listed in Results). If transformations failed to give an approximately normal distribution, non-parametric tests were used on the raw data (e.g., Wilcoxon’s test for pigment ratio). Distance measures were square root transformed.
To test for correlations between variables and MPB, satellite NDVI images were sampled to the 196 sample cells for consistency with other variables. P-values for correlation coefficients were calculated using the adjusted degrees of freedom from the Dutilleul et al.69 method, to counter the elevated risk of Type I statistical error with spatially autocorrelated data. Global (i.e., over the entire study grid) changes in variables from 2013 to 2014 early autumn sampling were tested by paired Wilcoxon tests. All tests were two-sided.
We used structural equation modelling (SEM) to assess which hypothesised ecological process(es) best explained the observed data70. SEM has previously been used to explain the impacts of oyster reefs on bird occurrences4, with a limited assessment of macrofauna. We contrast our main hypotheses based on pathway coefficients of four similarly-structured models, with 2013, 2014, and meiofauna data modelled separately to accommodate differences in the availability of variables and sites. The meiofauna model was dedicated to hypothesis (4), which expected a direct negative link between MPB and grazing, the latter inferred from nematode abundance, which could be mediated by copepod predation (inferred from abundance) and MGS. This was expected regardless of treatment impact, so hypothesis (4) was tested in the pre-treatment 2013 data only, which also halved the large workload involved in counting meiofauna. The other three models all tested hypotheses (1–3). Macrofaunal grazing was modelled as a latent variable approximated by the biomasses of M. balthica and S. plana. Variables are ‘latent’ when they are not directly observed but rather inferred from measured variables. The square root of distance from the reefs, reversed so that highest values were closest to the reefs, represented the non-linear distance decay of different processes according to each hypothesis. Hypothesis (1), local enrichment via oyster biodeposition, was not directly evidenced but emphasised a direct positive effect of distance from reefs on NDVI and OM. The other hypotheses emphasised indirect effects of distance from reefs on NDVI, via MGS (and OM) for the abiotic hypothesis (2), via macrofaunal grazers for the predation hypothesis (3), or via meiofaunal grazers for the meiofaunal (4) hypothesis. Hypothesis (3) expected grazer biomass to have a negative effect on MPB, while MGS, OM, and bathymetry were also allowed to affect distributions of both grazer groups.
Following the experimental treatment, we expected the link between distance from treatment reef (only) to NDVI to become insignificant under (1) or (3), visualised as a localised decrease in MPB relative to that surrounding the control reef. This was hypothesised because (1) killing the oysters halts biodeposition, or (3) because grazer abundance is no longer checked by predators, which were removed in the burning, although a recolonisation of the reef by predators may have occurred by the time of sampling (not recorded). Under (2), the experimental treatment is hypothesised to have no effect on MPB or sediment characteristics at the resolution of this experiment i.e., tens of metres (assuming the reef shell base remains intact; see Supplementary Fig. S18). Further details and SEM equations in Supplementary Methods: SEM. We use R package lavaan (v0.6-7)71. Parametric assumptions were tested for as above. As with our correlation analyses, we account for spatial autocorrelation by correcting SEM effective sample sizes for Moran’s I values using function ‘lavSpatialCorrect’ by Jarrett Byrnes (https://github.com/jebyrnes/spatial_correction_lavaan). This function corrects SEM standard error and significance values for spatial dependence. We assessed model fit by the comparative fit index (CFI), which ideally equals one, and the root mean square error of approximation (RMSEA), which ideally equals zero.
Mapping of ecological and environmental variables
The spatial structure in abundance, biomass, diversity, and environmental variables was mapped by geostatistics with the R package gstat (v2.0-6)72. We used Simpson’s 1 – D for diversity because it gives more weight to the common species we are interested in. Correlation strength between pairs of geographical samples (i.e., autocorrelation, characterised by semivariance) typically weakens with increasing distance between sample pairs (lag distance), a relationship plotted as the semivariogram. These are linear methods, so variables were detrended and normalised as explained above, except for the bivalve Scrobicularia plana biomass (Supplementary Methods: Geostatistical Details). Variograms used spherical and nugget models, which are well supported in ecology, firstly fitted by eye before secondarily using a weighted least-squares fitting algorithm to optimise sill and range estimates. The semivariograms were then used for simple Kriging interpolation, incorporating any significant trend coefficients back into the Kriging models as independent variables. Informing the Kriging model of these known trends (e.g., y-axis was expected to model bathymetry) allows prediction maps to not only reflect local variation but also natural gradients across the study extent. Analysis and mapping of biomass patterns were concentrated on macrofaunal species that contributed >1% of the total abundance. All mapping and analyses were performed in the statistical computing environment R (v4.0.2)73.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
