A tale of two key species in a subtropical mudflat: four-fold density increases produce minimal ecological response in macrofauna

Keywords: benthic ecology, biological traits​, bioturbation​, ecological functioning​, ecosystem services​, environmental management​, experiment, invertebrate​.

Soft-sediment intertidal flats provide many valuable ecosystem services, such as supplying food, fuel and construction materials, protecting coastal communities from storm surges and flooding, providing habitats for fish and birds, or cultural services, such as tourism and recreation (Barbier et al. 2011; Himes-Cornell et al. 2018). Many of the ecosystem services provided by mudflats depend directly or indirectly on the activities of benthic macrofauna (Crowe and Frid 2015; Passarelli et al. 2018). For instance, in soft sediments, macrofaunal assemblages mediate ecosystem functions, including primary and secondary productivity, carbon cycling, nutrient cycling, sediment stability, sediment oxygenation and bioturbation (Norling et al. 2007; Snelgrove et al. 2014). Large amounts of organic matter are delivered to estuaries from adjacent terrestrial and coastal ecosystems (Middelburg et al. 1995; Herbert 1999). This organic matter is then broken down by microbes present in intertidal sediments and nutrients are regenerated that are essential for primary productivity in the mudflats by microphytobenthos, and adjacent coastal ecosystems and so support wider marine food webs (Herbert 1999; Van Colen 2018). Feeding and bioturbation by macrofauna influence the rates of these activities, with nutrient sediment effluxes being ~four-fold higher when macrofauna are present (Karlson et al. 2007). The specific contribution of a species to ecosystem functioning varies depending on the unique characteristics, or biological traits, they possess (e.g. their morphological, behavioural and life-history attributes). Changes in biodiversity, and hence the mix of traits expressed, can therefore result in changes in functioning (Griffin et al. 2009). However, many questions remain regarding how the biological composition of an assemblage relates to its functioning (i.e. the biodiversity–ecosystem functioning (BEF) relationship; Hooper et al. 2005; Bulling et al. 2010).

So as to extrapolate our understanding of BEF relationships to ‘real world’ scenarios, we need to better understand how a broad range of species contribute to functioning, how their contribution varies across time and space, and what drives these differences. This knowledge is needed to interpret changes in natural systems against an increasingly complex patchwork of human activities (Díaz et al. 2003; Srivastava and Vellend 2005; Thrush et al. 2017). A recent review of research on soft-sediment habitats showed that despite being critical for coastal and offshore ecosystem functioning (Snelgrove 1997; Snelgrove et al. 2014), the BEF relationships in benthic communities are under-studied in many geographic areas (Dissanayake et al. 2018). Growing pressures from environmental change, such as climate and many other synergistic human pressures, are affecting the structure and functioning of ecosystems and the services they deliver to humans. If we wish to anticipate how natural systems will function in the future, we need to achieve a better understanding of the different ways in which ecosystems respond to anthropogenic change, and the extent to which they can compensate for any changes (Holling 1973; Ghedini et al. 2015; Landi et al. 2018; Thomsen et al. 2019).

One approach to understanding the BEF relationships and the degree of functional redundancy in an ecosystem (Oliver et al. 2015; Dissanayake et al. 2018) and, hence, its ability to withstand perturbation, is to quantify biological traits as a proxy for functioning (Bremner et al. 2003, 2006). The potential contribution of individual taxa to the processes underpinning ecological functioning can be derived from their biological traits. Experimental studies in laboratory mesocosms (Biles et al. 2003; Waldbusser et al. 2004; Braeckman et al. 2010; Karlson et al. 2016) and in the field (O’Connor and Crowe 2005; Clare et al. 2016) have examined how the biological traits of macrofauna influence functioning. Manipulative in situ field experiments can provide more realistic insights because they integrate the species interactions and abiotic processes within natural systems (Paine 1966; Dayton et al. 1974; Menge 1976; Wilson 1991; Raffaelli and Moller 1999; Booty et al. 2020; Hawkins et al. 2020). However, rigorous design is required to minimise artefacts of the experimental set-up, such as those owing to the presence of artificial structures, which may bias the outcomes (Hulberg and Oliver 1980; Woodin 1981; Hall et al. 1990; Underwood 1990; Peterson and Black 1994; Benedetti-Cecchi and Cinelli 1997; Como et al. 2006; O’Malley and Hunt 2020).

The present study used in situ experiments to examine the response of an intertidal macrofaunal assemblages to changes in the density of two large epibenthic taxa. The species selected for use were the Hercules club mud whelk, Pyrazus ebeninus (Bruguière; adult body length >3 cm and mean biomass was 41 mg per individual), and the Australian sentinel crab, Macrophthalmus (Mareotis) setosus H. Milne Edwards (adult body length >5 cm and mean biomass was 494 mg per individual). These species are abundant at the study site (McCoy’s Creek) and throughout in south-eastern Queensland mudflats (Dissanayake et al. 2020a) and are likely to make a significant contribution to several ecological functions (Bishop et al. 2007; Tanaka et al. 2013). Large, active consumers (grazers or predators) often play pivotal roles in the dynamics of communities (key stone species (sensu Paine 1966), and those that make a major contribution to bioturbation can substantially change the sedimentary environment (e.g. Snelgrove et al. 2014; Morais et al. 2019). The mudflat macrofaunal assemblage and its functioning (e.g. the primary productivity, bioturbation and oxygenation, and carbon decomposition) could respond to increased densities of these potentially ecologically important manipulated taxa because of the corresponding increase in physical disturbance and consumption of labile organic material. A lack of response, in contrast, would suggest that the effects of the increased biomass of these two large taxa were (i) offset by compensatory changes in the functioning of the resident assemblage, or (ii) this subtropical mudflat ecosystem can accommodate a greater biomass of these taxa without affecting the remaining assemblage (Naeem et al. 1999).

The response of the resident macrofaunal community and its contribution to key ecological functions were examined by increasing the density of the mud whelk Pyrazus ebeninus and the crab Macrophthalmus (Mareotis) setosus within enclosed experimental plots on mudflats at McCoy’s Creek, near the mouth of the Pimpama River, south-eastern Queensland (27.8222°S, 153.3778°E; Fig. 1a–c). Situated within a Marine National Park, the mudflat is fringed to landward by mangrove forests that are dominated by Rhizophora sp., Bruguiera sp. and Avicennia marina. Local land-use comprises a mixture of nature reserves, cultivated forests, marshlands and waterways making up 63% of the subcatchment area (that totals 10 × 10⁶ km²), with the remaining 37% being residential (Dissanayake et al. 2020a). A walking trail and residences are present within 750 m of the site; however, overall, it was minimally affected by direct human activities. We thus considered it suitable for the determining ecological functioning in a ‘natural’ system. The mean spring tidal range at McCoy’s Creek mudflats is 1.7 m (0.2–1.9 m above datum). The cages were installed in the mid-intertidal (~0.8 m above port datum Brisbane Bar) where the flats were exposed for ~4 h per tide (N. G. Dissanayake, pers. obs.). The site is typical of mudflats in Moreton Bay, comprising poorly drained soft muds with a silt–clay content ranging from 9 to 18%, 1.8 to 2.4% organic matter (by loss on ignition, LOI) and redox discontinuity layer (RDL) depth range of 1.5–3.5 cm. Measures of the biotic (e.g. macrofaunal density, diversity and taxonomic composition) and abiotic parameters (e.g. silt and clay percentage, standing stock biomass (chlorophyll-a concentration), organic matter content on LOI and RDL depth) between winter and summer at McCoy’s Creek mudflats are fully described in Dissanayake et al. (2020a, 2020b).

Fig. 1.  (a) Study site within the context of Australia. ACT, Australian Capital Territory; NSW, New South Wales; NT, Northern Territory; SA, South Australia; Qld, Queensland; Tas., Tasmania; Vic., Victoria; WA, Western Australia. (b) The study-site area at McCoy’s Creek subcatchment in south-eastern Queensland, Australia (27.8222°S, 153.3778°E) is highlighted by a black box and the water ways are indicated in grey. (c) The study site at the mouth of McCoy’s Creek opening to Pimpama river, indicating the three experimental set-ups in white dots. Experiments 1 and 2 ran from 10 to 26 June 2018 and from 25 August to 9 September 2018 respectively, for 16 days. The intertidal mudflats are illustrated in black strips and the mangrove habitat area in light grey. (d) Plan view of the enclosure cage (diameter 0.51 m, area of sediment enclosed 0.2 m²) and the 0.0625 m² area of sample retrieved at the end of the experiment. (e) Enclosure cage as deployed in the field. (f) Experimental three layout (not to scale) showing the randomised block design for the five treatments with varying densities of M. setosus (Ms) and P. ebeninus (Pe), that ran for 30 days from 29 August to 29 September 2019 and 90 days from 29 August to 25 November 2019. Treatments were IC, initial control; FC, final control (both without cage), CC, cage control; 4Ms:Pe, 4× the natural population density of M. setosus; Ms:4Pe, 4× the natural population density of P. ebeninus; 3Ms:3Pe, 3× the natural density of both species.

Pyrazus ebeninus is a non-selective epifaunal deposit feeder, which attains a length of up to 11 cm, and consumes organic detritus and algae on the sediment surface (Bishop et al. 2007). The whelks have been consumed by aboriginal people since at least the early days of European settlement and, aside from sea turtles, have few natural predators (McPhee 2017). As it crawls across the flats, the broad foot disturbs the upper layers of the sediments leaving distinctive tracks. The Australian sentinel crab M. (M.) setosus, which reaches up to 4-cm carapace length, is a deposit feeder consuming mostly microphytobenthos in both surface and subsurface sediments (Kanaya et al. 2008; Barnes 2010) where it constructs and inhabits burrows (Poore 2004; Supplementary Table S1). The contrasting biological traits of the two taxa suggest that they might make different contributions to mudflat functioning. Although it would have been valuable to complement these additions with some ‘trait removals’, we did not do so because the physical disturbance to sediments and non-target species would have been so great as to invalidate or undermine any experimental effects.

To assess whether the animal additions led to any changes in functioning, the impacts on both the macrofaunal communities present and some physico-chemical properties of the sedimentary environment were determined at the beginning and end of the experiments. Both the numbers and biomass of the macrofaunal species present were determined because animals of different body size may make different contributions to functioning (Thrush et al. 2006). Sediment organic matter content was compared because carbon decomposition reflects community metabolism (Wantzen et al. 2008; Cesar and Frid 2009), changes in the benthic chlorophyll-a concentration (i.e. the standing stock of microphytobenthos) was used as a metric of primary production potential (MacIntyre et al. 1996; Huot et al. 2007). The depth of the RDL was measured as a proxy for the degree of bioturbation and associated sediment oxygenation. The RDL depth is mostly controlled by the dynamic equilibrium between the downward transport of oxygen from the surface by bioturbators and the consumption of oxygen by organisms within the sediment (Kristensen et al. 2012; Snelgrove et al. 2014).

Initial experiments were conducted during 10–26 June (Experiment 1) and 25 August–9 September 2018 (Experiment 2), with sampling at 0 and 16 days for both. For Experiment 1, the densities of the crab M. setosus (Ms) and whelk P. ebeninus (Pe) were manipulated to produce treatments with one species occurring at natural densities and the other at double (2Ms:Pe and Ms:2Pe treatments), and both species at 1.5× natural density (1.5Ms:1.5Pe; Table 1). In Experiment 2, higher densities were used, with one taxon at ambient and the other at quadruple (4Ms:Pe and Ms:4Pe), and both taxa at 3× ambient (3Ms:3Pe; Table 1). The results from these initial experiments suggested that no changes in community composition or ecological functioning had occurred, and so a longer-duration experiment (Experiment 3) was implemented using the same density treatments as for Experiment 2, but was sampled after 30 and 90 days, to validate these results. The 3-month experimental duration was selected to increase the chances of detecting an ecological effect while minimising the potential for experimental artefacts, which tend to increase over time (Hall et al. 1990; Raffaelli and Moller 1999; Gallucci et al. 2008). Experiment 3 was conducted during spring of 2019, from 29 August to 25 November 2019, with sampling after 0, 30 and 90 days.

Table 1.  Macrofaunal densities used in experiments at McCoy’s Creek (27.8222°S, 153.3778°E; Fig. 1a–c), opening to Pimpama River, south-eastern Queensland, Australia, from June to September in 2018 and 2019.

In total, five replicates of five treatments were used (open control plots, caged controls and three different animal densities) with durations of 16, 30 and 90 days across the three experiments (Table 2, Fig. 1d). Cage controls, i.e. a cage enclosing an unmanipulated area of sediment, was established to account for any experimental artefacts that could be introduced by the cages (Hulberg and Oliver 1980; Woodin 1981; Peterson and Black 1994; Raffaelli and Moller 1999). The effects of experimental artefacts vary with habitat, type of manipulation and the target species, therefore it is necessary to include controls that may estimate such effects (Harrison 2004).

Table 2.  The 11 biological traits and 45 trait modalities used to characterise the macrofaunal taxa observed at McCoy’s Creek mudflats (27.8222°S, 153.3778°E), Pimpama River, south-eastern Queensland, Australia (Fig. 1a–c).

In June 2018, the ambient densities of the crab and whelk at McCoy’s Creek mudflat were determined by sampling 10 box cores. The mean (±s.e.) ambient densities were 10 ± 2 M. setosus individuals per 0.2 m² and 3 ± 1 P. ebeninus individuals per 0.2 m² (0. 2 m² being the area of a cage). To achieve the desired experimental densities, animals were added to cages (Table 1), which remained in place throughout the experiment. Using this approach, we had no prior knowledge of the initial in situ density of taxa in each cage (as that would have required destructive sampling), and so, the total number of animals enclosed within a cage varied among replicates.

The P. ebeninus used for the experimental additions were collected from mudflats at Jacobs Well (27.7747°S, 153.3613°E) ~5 km to the north (Fig. 1b), and M. setosus was collected from McCoy’s Creek (>400 m from the experimental site). Animals were collected 1 week before initiation of the field experiments and were maintained in the laboratory at 10°C (to reduce metabolic activity, ambient mean water temperature was 22°C) in aerated seawater collected from the site and renewed daily.

The experimental plots were enclosed in cylindrical mesh cages (Fig. 1d, e) constructed of aluminium ‘gutter wire’ (with 0.5-cm mesh) to maintain the required densities of the manipulated species. Cages were secured with cable ties and covered with wire-mesh lids reinforced with bamboo canes (Fig. 1e). This was sufficient to contain the juvenile and adult crabs (carapace width 0.8–4.2 cm) and mud whelks (shell length 5.0–6.0 cm) added to the experimental cages. We made no attempt to restrict the movement of small fauna, including post-larvae of the manipulated taxa. All cages were installed 15 cm deep into the sediment, and, so, assuming that none of the animals burrowed sediment below this depth they would have remained within the cages. Escape of the animals via burrowing seems unlikely because the previously recorded mean RDL was 3.3 cm deep (ranging from 0.5 to 6.5 cm) at McCoy’s Creek and was 4.5 cm deep in mudflats across the region (Dissanayake et al. 2020a). The cages enclosed ~0.2 m² of sediment surface each, meaning that when the box cores were extracted (0.0625 m² × 15 cm deep) to sample macrofauna at the end of the experiment, an 8–13 cm wide ‘buffer zone’ remained unsampled to minimise any edge effects (Fig. 1d).

The experimental protocol was specifically designed to allow analysis using fixed-factor ANOVA (time and treatment) following Underwood (1981), with additional consideration of pseudoreplication issues (Hurlbert 1984), although we recognise that such approaches are seen by some as dated (Beninger et al. 2012). To control for any spatial variation across the site, 10 experimental blocks (one replicate treatment in each block) were established parallel to the shore, placed separated by at least 2 m. The experimental blocks for treatments with differing temporal sampling points (i.e. after 30 or 90 days of duration) were assigned to alternate blocks to minimise the influence of any local environmental variation, such as shading by trees or small-scale variations in organic content. Within each block, the experimental cages or plots were situated at least 1 m apart and their relative positions were randomised (Fig. 1f). The three different experiments were set up >100 m from each other because it was clear that the site(s) had not fully recovered from the disturbance caused during the previous experimentation.

On the initiation of the experiment, the cages were installed during low tide on two consecutive days (Fig. 1c, f, Table 1). Macrofauna were added to their respective treatment cages and the covers were secured; this marked the beginning of the experiment for that cage.

Following installation, five open-control samples for Day 0 were retrieved using a box-core from undisturbed sediments adjacent to each experimental block. The RDL was measured using a ruler (to the nearest 1 mm) as each core was extracted. Sediments were sieved in situ over 0.5-mm mesh and the macrofaunal residues were preserved in 90% ethanol with Rose Bengal stain. For each control, three samples for chlorophyll-a analyses were collected from the sediment surface by inverting a disposable Petri dish and scraping off the underside to collect a sample of ~9-cm diameter × 1.4-cm depth (the size of the dish; Grinham et al. 2007). Three PVC sediment cores (3-cm diameter, 15-cm depth) were extracted for determination of sediment particle-size distribution and organic matter content by loss on ignition. All sediment samples were kept on ice in an insulated cool box to prevent degradation, and samples for chlorophyll-a were immediately wrapped in foil to exclude light. Once returned to the laboratory, all samples were frozen at −20°C.

Visual checks were made on the cages every 7–14 days to ensure their integrity and remove any accumulated debris to minimise the risks of physical damage or additional shading of plots. Half the plots were sampled after 30 days and the remainder after 90 days (Fig. 1f). During the low tide on two consecutive days, the enclosure cages were carefully removed, and box-core samples were extracted from the centre of each cage and control plot and were sieved and stored as described above (Fig. 1b). The redox discontinuity layer depth was measured, chlorophyll-a and sediment cores were collected from every plot as per the Day 0 sampling.

The preserved sieved residues from each box-core were carefully examined under a stereoscopic microscope and the macrofauna were enumerated and identified to species level (as far as practicable) by using standard Australian taxonomic works (Beesley et al. 1998, 2000; Glasby and Fauchald 2003; Australian Faunal Directory, see https://biodiversity.org.au/afd/home; Atlas of Living Australia, see https://www.ala.org.au/). Dissections were made to aid taxonomic identification where necessary (i.e. molluscs were dissected from their shells, polychaete chaetae were removed). Identifications were verified with taxonomic experts at the Museum and Art Gallery of the Northern Territory, the Queensland Museum and Museum Victoria, Melbourne. Scientific names were verified using, and follow, the World Register of Marine Species (https://marinespecies.org/).

The mean biomass was determined, by using a total of 10 individual animals (randomly selected) per taxon oven-dried at 60°C (until they reached a constant weight) and the dry weight was measured to 0.01 g. The mean dry weights were used to convert the species-abundance data to species-biomass data for each taxon, before the biological-trait analysis (see ‘Data analyses’ section).

Samples for analyses of physicochemical variables were thawed overnight at 4°C and were manually homogenised for 5 min. Chlorophyll-a was extracted by adding 40 mL of 90% reagent grade acetone to 20 g of wet sediment (maintaining a 1:2 ratio of sediment to solvent). The samples were wrapped in foil, mechanically shaken for 2 h, refrigerated (4°C) for 24 h (all in the dark). Samples were then centrifuged at 26 000g for 10 min and filtered through Watmann 20-μm papers (Grinham et al. 2007), and the absorbance was measured using spectrophotometry at 630-, 647-, 664-, 691- and 750-nm wavelengths and chlorophyll-a concentration was calculated (Ritchie 2008).

Organic matter content of sediment was determined using LOI, where 15 g of wet soil from each core was oven-dried at 105°C until a constant weight was obtained, and samples were combusted at 550°C in a muffle furnace. Dry weights were measured at both temperatures and the equation of Heiri et al. (2001) was used to calculate organic matter content as LOI at 550°C.

For sediment particle-size analyses, 50 g of oven-dried sediment (60°C for 48 h) was sieved across the full particle-size spectrum (2, 1, 0.710, 0.5, 0.25, 0.180, 0.125, 0.063, <0.063 mm; Percival and Lindsay 1997).

The macrofaunal assemblage was characterised by 11 biological traits representing morphological, life-history and behavioural aspects of each taxon’s biology. Each trait was characterised by 3–5 trait modalities (e.g. maximum adult body size was divided into the following four size groups: <10, 10–30, 30–50 and >50 mm, Table 2), giving a total of 45 modalities. The affinity of each taxon to each modality was assigned such that the ‘total’ for a trait summed to one. This ‘fuzzy coding’ (splitting the total among the modalities of a single trait) approach allowed both diversity and plasticity in the biology of the organisms, and any uncertainties to be captured (Chevene et al. 1994). Information on biological traits was obtained from the literature, in particular the Biological Traits Information Catalogue (MarLIN 2006), the Marine Species Identification Portal (http://species-identification.org, accessed 4 February 2020), BOLD systems (Ratnasingham and Hebert 2007), Polytraits (Faulwetter et al. 2014) and selected papers (e.g. Macdonald et al. 2010; Kristensen et al. 2012; Queirós et al. 2013). The species–traits matrix was multiplied by the species-abundance matrix to give an abundance-weighted score for each trait modality in each sample. Similarly, the species–traits matrix was multiplied by the species–biomass matrix to give trait modality (dry weight) biomasses for each sample.

The benthic macrofaunal assemblages from the different treatments were compared using the total number of individuals, species richness, Shannon–Wiener diversity and Pielou’s evenness, and the multivariate taxonomic and trait composition (both density and biomass weighted). The abundances of selected taxa and traits were also compared, with the latter providing estimates of potential ecological functioning on the basis of the distribution of selected trait modalities. Changes in the biomass of selected trait modalities (represented by the manipulated taxa) in the remaining community were compared among treatments (Table S1). Functioning was further considered by exploring differences in the physical environment by comparing sediment chlorophyll-a concentration, RDL depth, organic matter content and the proportion of silt and clay among treatments.

Comparisons were initially based on two-way ANOVA approaches (including ANOSIM) with experimental duration and treatment factors. To avoid pseudoreplication, the mean values were calculated for each cage, and replication was provided by the cages (Hurlbert 1984). Given their violation of the assumptions of parametric ANOVA, univariate biotic indices were compared using Kruskal–Wallis or Mann–Whitney U tests. Community composition was compared using non-metric multidimensional scaling (nMDS) on log(X + 1) transformed Bray–Curtis similarity, analysis of similarity (ANOSIM), and the similarity percentage routine (SIMPER) identified pairwise differences that contributed to the dissimilarities. Multivariate analyses were completed using PRIMER (ver. 6 Beta, Plymouth Marine Laboratory, UK), and univariate statistics were performed in IBM SPSS Statistics (ver. 28.0, IBM Corporation, see https://www.ibm.com/support/pages/downloading-ibm-spss-statistics-280).

Additionally, the sediment parameters in the open and cage controls were compared to determine any ‘cage effects’, and because none were detected the results from the open and cage controls were combined into a single ‘combined control’ treatment. Although there is a diversity of opinions on this approach (e.g. Queen et al. 2002; Colegrave and Ruxton 2017), this pooling will have increased power in any subsequent analyses.  The density treatments were compared with the combined controls to ascertain the impacts of macrofaunal manipulations on the measured biotic and abiotic variables. Analyses of changes in the macrofaunal assemblage were completed with both the manipulated taxa included and those excluded from the analysis to detect whether the individual additions affected the remaining community.

Given that the absolute number of the manipulated taxa in a plot was unknown, it follows that the actual density associated with each treatment was also unknown a priori; the treatments (Table 1) represent the number of animals added. It is therefore potentially informative to also examine the influence of the actual density, as determined at final sampling, of the manipulated taxa on the macrofaunal assemblage. The generalised linear models were used to examine the effects of the experimental duration (30 or 90 days; fixed factor), the number of recovered crabs, and the number of recovered whelks (random factors) on the macrofaunal community structure, such as number of individuals, species richness, Shannon–Wiener diversity, Pielou’s evenness, and taxonomic and trait composition (PCA1 scores for the latter two, both indexed by density and biomass).

During winter (June–August) to spring (September–November), the McCoy’s Creek mudflat macrobenthic assemblages were dominated by annelids, molluscs, and arthropods. In total, 4888 individuals and 54 taxa were recorded from the 7.5 m² of mudflat sampled across the three experiments. Of the 54 taxa (Supplementary Table S2), 96% were described to species level, with nemerteans and nematodes being recorded to phyla level only. Over 54% of the individuals observed were accounted for by three species, including the predatory polychaete Aglaophamus australiensis and the deposit-feeding polychaete Barantolla lepte (24 and 15% of all the individuals respectively), and a suspension-feeding bivalve Hiatula alba (15% of all individuals). Six species, the gastropod Pyrazus ebeninus, the crabs Mictyris longicarpus, Gelasimus vomeris and Macrophthalmus setosus, the bivalve H. alba, and the polychaete A. australiensis comprised 90% of the biomass of the assemblage. Sediments comprised of 9.8–19.5% silt- and clay-sized particles, 1.8–5.7% organic matter and 8.3–12.6 mg chlorophyll-a m⁻², with a RDL depth of 2.0–3.8 cm (Table 3), none of which differed among seasons (Mann–Whitney U tests, P > 0.05).

Table 3.  Median macrofaunal biodiversity (±inter quartile range, IQR) and mean sediment parameters (±s.e.) across the open controls sampled in Experiment 1 (June 2018), Experiment 2 (August–September 2018) and Experiment 3 (August–September and November 2019) at McCoy’s Creek mudflat (27.8222°S, 153.3778°E, Fig. 1a–c), south-eastern Queensland, Australia.

The macrofauna were sampled between mid- and late winter (July–August) and early spring (September) in two consecutive years (2018–2019), and so the presence of any inter-annual or seasonal variations in assemblage composition could have confounded the results. To determine the likelihood of this, the unmanipulated open controls were compared among years and seasons.

Considering only the unmanipulated open controls, no significant interannual variations were found for macrofaunal density, species richness, Shannon Wiener diversity or Pielou’s evenness (Mann–Whitney U tests, P > 0.05) among experiments running through the same season. However, macrofaunal taxonomic composition indexed by density (ANOSIM, global R = 0.269, P < 0.05; Supplementary Fig. S1a) and biomass (ANOSIM, global R = 0.366, P < 0.05; Supplementary Table S3, Fig. S1b) significantly differed between years, with an average SIMPER dissimilarity of 45–48%. The median densities of four species significantly differed between 2018 and 2019, with the bivalve Laternula anatina (Mann–Whitney U tests, Z = −2.562, P < 0.05), the polychaete Owenia australis (Z = −1.993, P < 0.05) and the decapod Trypaea australiensis (Z = −2.509, P < 0.05) having lower densities in 2019, whereas the decapod Mictyris longicarpus (Z = −3.375, P < 0.05) had higher median density in 2019.

No significant seasonal shifts (Kruskal–Wallis test, P > 0.05) were observed in taxonomic diversity (species richness, Shannon–Wiener or Pielou’s evenness), including the manipulated taxa, among the open controls from across the seasons (late winter to late spring). Even though the taxonomic composition did not differ when indexed by biomass (ANOSIM, P > 0.05, Fig. S1c), a significant seasonal shift was observed when indexed by density (manipulated taxa included; ANOSIM, global R = 0.216, P < 0.05, Fig. 2a). After 90 days, the open controls differed from both Day 0 (pairwise ANOSIM, R = 0.384, P < 0.05) and Day 30 (pairwise ANOSIM, R = 0.256, P < 0.05), with an average SIMPER dissimilarity of 44–47%. The polychaete Phyllodoce novaehollandia, and decapods Australoplax tridentata and M. setosus, increased in abundance between late winter and early spring, whereas the abundance of bivalve H. alba decreased (Fig. 2c). A significant seasonal shift in macrofaunal trait composition (weighted by density, manipulated taxa included) was observed for Experiment 3 between 0 and 30 days only (ANOSIM, global R = 0.193, P < 0.05, Fig. 2b–d, see the ‘Temporal shifts in the macrofaunal assemblage at McCoy’s Creek’ section of the Supplementary material). No seasonal shifts in biological-trait composition were observed for the biomass-weighted trait composition (one-way ANOSIM; P > 0.05, Fig. S1d).

Fig. 2.  nMDS ordination of the Bray–Curtis similarity of the macrofaunal (a) taxonomic composition (log(X + 1) transformed) and (b) trait composition, indexed by density, in the open controls in late winter (Day 0), early spring (after 30 days) and late spring (after 90 days) of the experiment at McCoy’s Creek mudflat, Pimpama River, south-eastern Queensland, Australia (27.8222°S, 153.3778°E). Median (±inter-quartile range), (c) densities (per 0.06² 5 m²) of the 14 species, and (d) the 16 biological trait modalities that contributed >70% to the SIMPER dissimilarity between the open controls in late winter (Day 0), early spring (after 30 days) and late spring (after 90 days). Asterisks indicate significant differences (at P = 0.05) for species and trait modality densities: *between 0 and 30 days; **between 0 and 30 days and 0 and 90 days; and any names highlighted in bold differed between 30 and 90 days (Mann–Whitney U test).

None of the physico-chemical parameters measured significantly differed between the cage controls and the open controls on any sampling date (Supplementary Table S4), and, so, the cages do not seem to have introduced any physical artefacts on the enclosed sediments. Similarly, neither the taxonomic (two-way nested ANOSIM, for time_{density/biomass}, P > 0.05; for cage effect_{density/biomass}, P > 0.05; Fig. 3a, b) nor the trait composition (two-way nested ANOSIM, for time_{density/biomass}, P > 0.05; for treatment_{density/biomass}, P > 0.05; Fig. 3c, d) of the macrofaunal assemblage significantly differed between the cage and the open controls after 16, 30 or 90 days. Given the absence of cage effects, data from the open and cage controls were combined for each sampling date and are hereafter referred to as the ‘combined controls’.

Fig. 3.  nMDS ordination of the Bray–Curtis similarity of the macrofaunal taxonomic composition (log(X + 1) transformed) for the open controls and cage controls after 16, 30 and 90 days of the experiment when indexed by animal (a) density, and (b) biomass. Macrofaunal trait composition for the open controls and the cage controls after 30 and 90 days when (c) density and (d) biomass weighted. Mean (±s.e.) number of crabs M. setosus and whelks P. ebeninus per 0.0625 m² area observed in the combined controls and the three density treatments (used to determine the effectiveness of the treatments) in (e) Experiment 1 and (f) Experiment 2 after 16 days; and, (g) Experiment 3 after 30, and (h) 90 days at McCoy’s Creek mudflat, south-eastern Queensland, Australia (27.8222°S, 153.3778°E).

The mean densities of the manipulated taxa recovered from the sampled area of each experimental plot (the central 0.0625 m²) at the end of Experiment 1 showed that neither M. setosus (ANOVA, P > 0.05) nor P. ebeninus (P > 0.05) densities significantly differed between the combined controls and the three density treatments (Fig. 3e, Supplementary Table S5). This suggested that the desired densities were not maintained in the central, sampled, area for either species. Casual observation showed that the whelks were aggregating around the cage perimeter and on the cage mesh (N. G. Dissanayake, pers. obs.).

For Experiment 2, the densities of M. setosus and P. ebeninus were increased up to 4× their ambient density (200 M. setosus individuals m⁻², 60 P. ebeninus individuals m⁻²), and the resulting mean densities recovered from the centre of the cage after 16 days were higher in the treatments than in the combined controls (Fig. 3f), although whelks continued to aggregate around the perimeter and on the mesh. These differences were significant for both species (ANOVA, F_(4,25) = 6.604, P < 0.05 and F_(4,25) = 32.126, P < 0.05 respectively; see the ‘Effectiveness of the density manipulations’ section and Table S5 of the Supplementary material). Thus, these densities were applied in Experiment 3, which had an extended duration. In Experiment 3, the recovered mean densities of the crabs and whelks were verified as being higher in the treatments than in the combined controls after 30 (crabs: ANOVA, F_(4,25) = 7.54, P < 0.05; whelks: F_(4,25) = 11.22, P < 0.05, Fig. 3g, Table S5) and 90 days (crabs: F_(4,25) = 2.80, P < 0.05; whelks: F_(4,25) = 25.32, P < 0.05, Fig. 3h; see the ‘Effectiveness of the density manipulations’ and Table S5).

When M. setosus and P. ebeninus were added to experimental plots at 3–4× their ambient density, neither the median macrofaunal density, species richness, Shannon–Wiener diversity, Pielou’s evenness (Kruskal–Wallis test, P > 0.05, Fig. 4a–d) nor the taxonomic composition of the remaining macrofaunal community (i.e. that which excluded the manipulated individuals; ANOSIM; density and biomass P > 0.05, Fig. 4e–j, Supplementary Table S7) significantly differed between the combined controls and density treatments. Nor did these parameters differ among the three density treatments after 16 days in Experiment 2, or after 30 or 90 days in Experiment 3 (Table S7). These comparisons were made after excluding the abundances of the manipulated taxa from the assemblage data, so that only the remaining, and not manipulated, fauna was compared.

Fig. 4.  Median (±interquartile range, IQR) density and taxonomic diversity of the macrofaunal assemblages retrieved from the centre of the four treatment cages after 16, 30 and 90 days at McCoy’s Creek mudflat (27.8222°S, 153.3778°E), south-eastern Queensland, Australia, in 2018–2019. Treatments were: the combined controls; 4× the mean ambient density of M. setosus (4Ms:Pe); 4× the mean ambient density of P. ebeninus (Ms:4Pe); and, 3× the mean ambient density of both M. setosus and P. ebeninus (3Ms:3Pe). (a) The median number of individuals (N) per 0.0625 m², (b) median species richness (SR), (c) median Shannon–Wiener diversity (H′), and (d) median Pielou’s evenness (J) of the macrofaunal assemblages for each treatment (based on five replicate cores). nMDS ordination of the Bray–Curtis similarity (log(x + 1) transformed) of the density weighted macrofaunal taxonomic composition among the four treatments after (e) 16 days in Experiment 2, (f) 30 days and (g) 90 days in Experiment 3. The biomass weighted macrofaunal taxonomic composition (excluding manipulated taxa) among four treatments after (h) 16 days in Experiment 2, (i) 30 and (j) 90 days in Experiment 3.

When the manipulated taxa were included, the trait composition weighted by density did not differ among the four treatments after 16, 30 or 90 days (ANOSIM, P > 0.05, Supplementary Fig. S4a–c, Table S7). The trait composition when indexed by biomass significantly differed between the combined controls and the three density treatments after 16 days (ANOSIM, Global R = 0.537, P < 0.05) of Experiment 2, after 30 (ANOSIM, Global R = 0.291, P < 0.05) and 90 (ANOSIM, Global R = 0.418, P < 0.05, Fig. 5a–c, Table S7) days of Experiment 3. Pairwise ANOSIM showed the trait composition of the combined controls and 4Ms:Pe density treatment did not differ after 16 and 30 days (ANOSIM, P > 0.05, Fig. 5a, b, Supplementary Table S6), but both differed from the treatments that included elevated densities of the whelks (ANOSIM, P < 0.05, Fig. 5a, b, Table S6). After 90 days, macrofaunal trait composition (Fig. 5c, Supplementary Table S8) of the combined control differed from the Ms:4Pe treatment only (P < 0.05). The high whelk-density treatment had a trait composition significantly different from 4Ms:Pe (P < 0.05, Fig. 5c, Table S8) and 3Ms:3Pe (P < 0.05, Fig. 5c, Table S8), and the high crab treatment (4Ms:Pe) also differed from 3Ms:3Pe (P < 0.05, Fig. 5c, Table S8). Fourteen trait modalities contributed >70% of the cumulative SIMPER dissimilarity among the treatments after 16, 30 and 90 days (Fig. S4d–f, see the ‘The influence of the manipulated species on the resident macrofaunal assemblage’ section of the Supplementary material).

Fig. 5.  nMDS ordination of the Bray–Curtis similarity of the biomass-weighted macrofaunal trait composition (manipulated taxa included) among the four treatments (the combined controls; 4× the mean ambient density of M. setosus (4Ms:Pe); 4× the mean ambient density of P. ebeninus (Ms:4Pe); and 3× the mean ambient density of both M. setosus and P. ebeninus (3Ms:3Pe) at McCoy’s Creek mudflat (27.8222°S, 153.3778°E), south-eastern Queensland, Australia, after (a) 16 days in Experiment 2, (b) 30 and (c) 90 days in Experiment 3. Mean (±s.e.) trait modality biomass (log(X + 1) transformed) including and excluding manipulated taxa for (d) feeding, (e) sediment dwelling depth and (f) bioturbation traits in the four treatments after 90 days (Experiment 3). Asterisks (*) denote trait modalities that differ significantly (Tukey’s post hoc multiple comparisons showed that, in each case, only the MS:4Pe treatment differed from all others).

The mean biomasses of deposit feeders, surface modifiers and those with sediment dwelling depths of 0–2 mm were higher in the high whelk-density treatments than in all of the others after 90 days (ANOVA, P < 0.05, Fig. 5d–f), reflecting the greater biomass of the traits of the whelk (Table 2). The macrofaunal additions did not significantly alter the trait composition of the remainder of the assemblage whether weighted by density or biomass (when the manipulated taxa were excluded) for Experiment 2 (ANOSIM, P > 0.05, Supplementary Fig. S3a, d and Table S7) or Experiment 3 after 30 or 90 days (ANOSIM, P > 0.05, Fig. S3b, c, e, f). None of the other biological traits differed among treatments (Fig. 5d–f). When the biomasses of the manipulated taxa were excluded, none of the biomasses of traits for the remaining community differed between the combined controls and the density treatments (ANOVA, P > 0.05, Fig. 5d–f).

Generalised linear model performed on the macrofaunal assemblages (manipulated taxa excluded) after removing a single outlier (Supplementary Fig. S5). The analyses confirmed that the experimental duration and number of recovered crabs and whelks did not influence the structure of the unmanipulated component of the macrofaunal assemblage after 30 or 90 days (Wald chi-square test, P > 0.05, Supplementary Table S9).

The mean sediment chlorophyll-a concentration and the percentage of silt-and clay-sized particles did not significantly differ between the combined controls and the three density treatments at any time (ANOVA, P > 0.05, Fig. 6a, d). The mean organic matter content did not significantly differ between the combined controls and density treatments after 16 and 30 days (ANOVA, P > 0.05), but after 90 days, it was 78–109% greater in all three treatments than in the combined controls (ANOVA, Tukey’s post hoc test, P < 0.05, Fig. 6b). The redox discontinuity layer was 30% deeper than the combined controls in the two most extreme treatments (4Ms:Pe and Ms:4Pe, Fig. 6c) after 30 days (Tukey’s post hoc, P < 0.05), and 70–95% deeper in all treatments after 90 days (Tukey’s post hoc test, P < 0.05), suggesting that more bioturbation may have been occurring, and more organic matter was accumulating in the treatments towards the latter stages of the experiment.

Fig. 6.  Mean (±s.e.m.) sediment parameters in the presence of increased densities of M. setosus (crab) and P. ebeninus (whelk) after 16, 30 and 90 days of the experiments at McCoy’s Creek (27.8222°S, 153.3778°E), south-eastern Queensland, Australia. (a) The benthic chlorophyll-a (Chl-a) concentration (mg m⁻²), (b) organic matter (OM) content (percentage on loss on ignition at 550°C), (c) redox discontinuity layer (RDL) depth (cm), and (d) the percentage of silt- and clay-sized particles (<0.063 mm). Treatments: 4Ms:Pe, 4× the mean ambient density of M. setosus; Ms:4Pe, 4× the mean ambient density of P. ebeninus; 3Ms:3Pe, 3× the mean ambient density of both M. setosus and P. ebeninus.

Coastal ecosystems are under intense pressure from anthropogenic activities (Williams et al. 1990; MEA 2005; IPCC 2014; McPhee 2017; Newton et al. 2020), and an improved understanding of their dynamics is required to anticipate how marine ecosystems, and the services they deliver, might change in the future. In situ experimental density manipulations of organism densities have shown how communities and the associated ecosystem functions respond to change (Wilson 1991; Bertness et al. 2014; Ling et al. 2015; Aguilera et al. 2018; Booty et al. 2020). Therefore, we hypothesised that the macrofaunal assemblage structure, their collective biological traits (and, by extension, mudflat functioning), and the abiotic environment would respond to increased densities of the crab Macrophthalmus setosus and the mud whelk Pyrazus ebeninus. Within the subtropical mudflats in Queensland, large increases in densities of these two taxa did not produce changes in the structure or biological-trait composition of the remaining macrofaunal assemblage, leading us to accept the null hypothesis of no detectable effect on the assemblage. The community did not resist (sensu Dayton et al. 1974) the experimental animal additions. However, we were able to reject the null hypothesis for some sediment properties; RDL was 70–93% deeper, and organic matter content was 75–100% higher after 90 days in treatments with elevated densities of crabs or whelks, than in the combined controls.

Several previous studies of both present-day and past benthic assemblages have shown that, contrary to simplistic expectations, loss of species or individuals from an assemblage often results in patterns of species change that conserve trait composition and, so, ‘potential’ ecological functioning (Cesar and Frid 2012; Frid and Caswell 2016; Caswell and Frid 2017; Meyer and Kröncke 2019; Gammal et al. 2020). In the present study, the abundances of the unmanipulated taxa did not measurably respond to the animal additions. However, as the assemblage gained individuals (from the manipulations), the functioning the McCoy’s Creek mudflats differed from the pre-manipulation configuration. Therefore, in this case, functions such as organic matter decomposition, bioturbation and oxygenation were not conserved. Rather, the absolute amount of some functions increased following our experimental additions and this altered state persisted for at least 90 days. However, because species composition, and hence traits, did not change, neither would the range of potential functions.

Direct measurement of ecological functioning is challenging (Bolam et al. 2002; Gammal et al. 2017) because it involves many factors related to the chemical, physical and biological components of the system (Hooper et al. 2005). Changes in ecological functioning may be stimulated by experimentally manipulating taxa with presumed important functional roles (based on their biological traits; Duffy et al. 2003; Thrush et al. 2006; Cardinale et al. 2011; Norkko et al. 2013; Gammal et al. 2017), which was the basis of the approach taken in the present study. For instance, a decrease in microphytobenthic biomass has been observed when large deep-burrowing deposit feeders are removed from some systems, whereas removal of large surface-dwelling suspension feeders can increase sediment organic matter content, and, so, its decomposition (Thrush et al. 2006). In the present study, increases of up to 4× the natural densities of the surface deposit-feeding whelk (up to 60 individuals m⁻²) and burrowing subsurface deposit-feeding crab M. setosus (up to 200 individuals m⁻²) did not significantly alter the remaining macrofaunal assemblage. The manipulated macrofaunal densities were maintained (as far as possible for the purpose of this study) throughout the experiment and the lack of any difference between the caged and open controls suggests that there were no experimental artefacts, and predators (of a size that could be excluded by the cages) do not seem to have been a major factor at this time of year. As a minimum, this suggests that the additional food requirements of these extra individuals were met without measurably affecting the abundances of other macrofaunal taxa. Whereas we did not measure the availability of food for these deposit feeders directly (i.e. we did not distinguish the labile organic matter from the total), we do know that the total organic matter in McCoys Creek sediments were within the range of other mudflats on the eastern coast of Australia where these species occur (Morelli and Gasparon 2014; Dissanayake et al. 2020a).

The dietary breadth and plasticity of ‘deposit feeders’ on mudflats makes unravelling the trophic pathways difficult (Melville and Connolly 2005; Ehrnsten et al. 2019, 2020), and reduces the likelihood of strong competitive interactions. The magnitude of manipulation at McCoy’s Creek may not have been sufficient to induce strong interspecific competition. However, a similar experiment at Quibray Bay, Sydney, (Bishop et al. 2007), which manipulated P. ebeninus at lower densities than those used in the present study, found that after 2 months, the macrofaunal abundance and species richness of the resident assemblage was halved for plots with 44 whelks m⁻² compared with those with just 4 whelks m⁻². Microphytobenthos standing stock was also 20% lower when whelk density was high. The natural whelk densities at Quibray Bay were ~17% higher, and the surface sediments were mesotrophic, compared with McCoy’s Creek (70–90 mg chlorophyll-a m⁻² and 3–14 mg m⁻² at the two sites respectively). Although the microphytobenthic standing stock at McCoy’s Creek was at the lower end for intertidal flats composed of fine cohesive sediments, which can range from 1 to 560 mg chlorophyll-a m⁻² globally (Underwood 2001), it was typical for other muddy-sand shores across Moreton Bay (3–20 mg m⁻² in winter and 4–22 mg m⁻² in summer; Grinham et al. 2007; Dissanayake et al. 2020a, 2020b). The lower autochthonous production at McCoy’s Creek might have made the macrofauna more dependent on non-living organic matter (including the contribution from the mucus and faecal material produced by the greater numbers of whelks).

Although there was no evidence for changes in the remaining macrofauna at McCoy’s Creek, the physical environment did change in a way similar to that observed in other systems where bioturbator density is altered (e.g. Otani et al. 2010; Morais et al. 2019). Epifauna, such as foraging mud snails, e.g. P. ebeninus (Edwards and Welsh 1982; DeWitt and Levinton 1985) and large vertebrates that disturb surface sediments (e.g. shorebirds; Booty et al. 2020) or burrowers can affect surface deposit feeders by decreasing microphytobenthic production (Botto and Iribarne 1999; Webb and Eyre 2004; Lomovasky et al. 2006; Volkenborn and Reise 2007; Volkenborn et al. 2007). For example, the complete removal of Macrophthalmus japonicus from temperate mudflats of the Tama Estuary, Japan, during autumn increased mean macrofaunal density by 53% in the subsequent summer compared with controls. Densities of the surface deposit-feeding polychaete Hediste sp. and the filter-feeding bivalve Corbicula japonica increased in the crab’s absence (Tanaka et al. 2013). Frequent burrowing by M. japonicus destabilised surface sediments and reduced microphytobenthic production, which inhibited Hediste sp. and C. japonica feeding. The consequences of removing M. setosus were not explored at McCoy’s Creek, but increased M. setosus density and the corresponding ~2 cm deeper sediment redox-potential discontinuity layer (RDL) did not drive any shifts in macrofaunal taxonomic or trait composition. It is possible that decreased bioturbation would have stimulated a response from the other taxa, as occurred with M. japonicus, once it dropped below a particular threshold. More detailed research into these mudflats is needed to establish the trophic relationships and pathways and to understand how they can support such large increases in consumer biomass for >90 days with no detectable impacts on the remaining assemblage.

More detailed research into these mudflats is needed to establish the trophic relationships and pathways and to understand how they can support such large increases in consumer biomass for >90 days, with no detectable impacts on the remaining assemblage. The impacts of the density manipulations on chlorophyll-a might be confounded by small grazers <0.5 cm (the cage mesh size) that were able to move into and out of the cages.

A biological-traits approach describes ‘potential’ ecosystem functioning (Bremner 2008; Cesar and Frid 2012) and previous studies (Thrush et al. 1997; Cesar and Frid 2009; Clare et al. 2016) have identified mismatches between potential and actual functioning because the traits of a taxon used in the coding are generalisations (Bremner 2008; Tyler et al. 2012). In practice, organisms show considerable plasticity in their behaviours and, so, only a portion of the potential functioning may be expressed at any one time. The magnitude of potential functioning increased when taxa were added to the experimental cages at McCoy’s Creek, and the increased sediment organic matter content and deeper RDL in cages with elevated densities of the whelk and crab suggest that some changes in actual functioning occurred. Specifically, sediments accumulated more carbon, perhaps as a result of increased whelk faecal pellets and mucus production, and were bioturbated to a greater depth.

Given the lack of significant changes in the abundance and composition of the remaining macrofaunal assemblage, it seems likely that the altered sediment OM load and deeper RDL were a direct result of the foraging or burrowing behaviour of the two large taxa added, although the effect did not differ between the high-crab and high-whelk treatments. However, it is possible that individuals across the community responded to the additions by altering, within their envelope of plasticity, the amount and types of ecological functions expressed (e.g. Cesar and Frid 2009). This could have limited the scale of the changes observed by partially masking or mitigating the effects of the added taxa. The strength of this effect might be expected to change over time, as different community level processes began to manifest, such as reproduction and recruitment, or vary with seasonal changes in the community and resource availability (Raffaelli and Moller 1999; Underwood 2000; Jenkins and Uyà 2016; Meyer and Kröncke 2019).

So as to anticipate how mudflat ecological functioning will change under anthropogenic pressure and apply this knowledge to management, we need a better understanding of trophic interactions and BEF relationships across spatial and temporal scales. Experimental manipulations of animals in situ can help us better understand these relationships. Changes in the abundances of experimental animals in the McCoy’s Creek mudflats increased the prevalence of some biological traits (at a minimum those provided by the added individuals) and will presumably have altered competitive (and predator–prey) interactions. The assemblage accommodated the enhanced densities of the crab and whelk for the duration of the experiment, although we suspect that over time the system would ultimately revert to the pre-manipulated state as the introduced individuals were lost by natural mortality. This suggests that these systems may have the capacity, at least at certain times of the year, to support greater macrofaunal biomass and, so, deliver higher levels of ecosystem services.

The data used to generate the results in the paper will be available in PANGEA (https://www.pangaea.de), data publisher for earth and environmental science. Supplementary results for this paper are also available in the Supplementary material.

The authors declare that they have no conflicts of interest.

The work presented here was funded by Griffith University International Postgraduate Scholarship to N. G. Dissanayake from the School of Environment and Science, Griffith University, and the Environmental Futures Research Institute, Queensland.

Supplementary material is available online.