Experimental design
Reef structural complexity metrics—fractal dimension and height range—were manipulated using artificial habitat units, of standardized 15 × 15 cm planar area, but of varying three-dimensional geometry and surface area. There were 16 distinct designs (Fig. 1c), produced by crossing four levels of each fractal dimension and height range to produce 11 resultant levels of surface area. The range of levels of fractal dimension (2.02−2.52), height range (3−12 cm) and surface area (244.75−1,467.79 cm2) were centred around, but extended beyond, values observed for S. glomerata oyster reefs in Sydney, New South Wales, Australia (mean ± s.e.m., fractal dimension = 2.31 ± 0.009, height range = 7.67 ± 0.40 cm, surface area = 571.01 ± 17.22; Fig. 3).
Values of fractal dimension, height range and surface area were calculated from digital elevation models (DEMs) using the habtools package of R (ref. 55); DEMs for all artificial and natural surfaces were generated from three-dimensional models using the ‘mesh_to_dem’ function of habtools. All DEMs had the same resolution (1 mm), to allow for comparisons between surfaces. Fractal dimension was estimated using the fd() function and the variation method, based on measurements taken across six spatial scales obtained by subdividing each DEM into squares with side lengths of 15, 7.5, 3.75, 1.875, 0.9375 and 0.46875 cm. The variation method was used because it allows estimating the fractal dimension of real surfaces56, which do not necessarily behave like pure fractals (for example, they are rarely self-similar across scales). Height range was calculated using the hr() function, and surface area was calculated using the surface_area() function. Three-dimensional models of oyster reefs (15 × 15 cm, n = 28), were obtained using photogrammetry together with structure from motion52 (Supplementary Methods) from four naturally occurring reefs located at Towra Point Nature Reserve, New South Wales, Australia (34°0′56.21″ S, 151°10′ 44.57″ E). Reefs from Towra Point were selected because, at the time of the study, this was the only known location in the Greater Sydney region that contained multiple large undisturbed S. glomerata reefs39.
Three-dimensional models of the 16 artificial habitat units were designed using the three-dimensional graphics software Blender. Each of the 16 artificial habitat designs was three-dimensionally printed in polylactic acid filament and used as a master to create silicon moulds that allowed casting replicate units using concrete. A total of 500 concrete habitat units were produced (480 habitat units for hypothesis testing and 20 habitat units for assessing caging artefacts; see the next section). Concrete releases positive settlement cues for oysters, is easily moulded into a range of geometries and is used as a substrate in restoration projects44,57.
Field experiment
To assess the effects of habitat complexity descriptors on oyster recruitment, we used habitat units (n = 10 for each of the 16 designs) at three estuarine sites in the greater Sydney region, New South Wales, Australia, in October and November 2022: Brisbane Waters (33°30′40.29″ S, 151°10′18.75″ E), the Hawkesbury River (33°33′46.16″ S, 151°13′46.70″ E) and Port Hacking (34°04′41″ S, 151°06′38.03″ E). These sites are adjacent to S. glomerata oyster reefs, indicative of larval supply and support a diverse range of oyster predators58 (Supplementary Table 7). At each of the sites, five units of each design were caged (1 cm mesh; to exclude finfish predators) and five were uncaged (so they could be accessed by predators). At each of the three sites, we deployed units at a mid-intertidal elevation (0.6–0.9 m above the lowest astronomical tide; daily inundation period of about 56%), with units haphazardly interspersed with respect to treatment along the shoreline and separated by at least 15 cm from neighbouring units or natural habitat features (for example, rocks). After twelve months—a period long enough to allow for both settlement and post-settlement mortality to occur—we counted the number of oyster recruits on the surface of each experimental unit. This required scraping oysters out of habitat features to accurately enumerate them, preventing repeated sampling. To test for caging artefacts, partial cages that allow predator access were established for a subset of four designs with five replicates each at one randomly selected site (Brisbane Waters; Supplementary Table 8). We found no effects of caging artefacts in oyster recruitment (Supplementary Table 8). Research activities were performed under the New South Wales Department of Primary Industries scientific collection permit P07/0047-7.0 and approved by Macquarie University Animal Ethics Authority (2019/036).
Statistical analyses
Surface area and height range were base-10-log-transformed to improve model residuals. Oyster densities were square-root-transformed to improve model residuals. Oyster densities were calculated by dividing oyster counts found at each experimental unit by the surface area of the experimental units. We quantified habitat structural complexity–abundance relationships using a GLMM (for count data) and LMMs (for density data), with second-order polynomials to allow for non-linear relationships between predictor and response variables. Study sites were included as random effects. To test our hypothesis that increasing habitat structural complexity would confer benefits to oysters that extend beyond increasing surface area for larval settlement, and that most of these additional benefits would arise from predator mitigation, we first compared the effects of surface area on oyster counts (Fig. 2a) between habitat units with and without cages using a GLMM with a negative binomial distribution. We then compared the effects of both fractal dimension and height range on oyster densities at habitat units with and without cages (Fig. 2b,c) using LMMs to highlight the influence of predation in driving habitat structural complexity–abundance relationships. To test our hypothesis that ecosystem engineering maximizes recruit survival per unit area, we quantified the interactive effect of fractal dimension and height range (for uncaged units only) on oyster densities using a LMM (Fig. 3 and Supplementary Table 5). Residuals for all these models were approximately normal and were homogeneous when plotted against predictor variables. To test for caging artefacts, a zero-inflated generalized linear model was run using a negative binomial distribution to assess treatment effects (n = 2; uncaged and partially caged) on oyster counts. All analyses were conducted in R (ref. 59 ; v.3.6.1) and can be accessed at https://github.com/juan-muelbert/Oyster-complexity.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.