Ray bioturbation rates suggest they shape estuary processes

Abstract

Bioturbation of sediments is a key ecosystem service in estuarine and marine ecosystems, and rays (superorder Batoidea: skates, stingrays, electric rays and shovelnose rays) are among the largest bioturbators, modifying their habitat through foraging and predation. Ray activities cycle nutrients, increase oxygen penetration and re‐stratify sediments. However, given rays are globally threatened, it is unclear to what extent the loss of rays, and the ecosystem services they provide, would affect ecosystem processes. This study assessed the likely amount of sediment displaced annually by rays during foraging activities at an estuary scale. To achieve this, an aerial drone was used to map daily ray bioturbation activity, as evident from the presence of feeding pits. High‐resolution, 2.6 cm px⁻¹, Digital Elevation Models (DEM) were created and used to measure the volume of sediment displaced by feeding pits. We found rays within the Brisbane Water estuary in NSW, Australia excavated 1.20 (±0.68) tonnes of sediment per day within this 1443 m² intertidal area, or a rate of 575.2 cm³ m⁻² per day. This bioturbation rate is relatively high compared to bioturbation rates documented in other ray species. Spatial autocorrelation analysis indicated that the ray feeding pits were significantly clustered by location as well as size (P < 0.01), suggesting size segregation of ray foraging. When bioturbation rates were conservatively extrapolated across measured feeding area in the estuary, we calculated rays displace 57.6 (±32.4) tonnes of sediment per day, or 21.0 (±11.4) kilotonnes per year. This underlines how ray bioturbation likely shape estuary processes, and how the loss of rays and their ecosystem services would likely have considerable impacts on estuarine sedimentary ecosystems. These findings emphasize the ecological importance of rays in an ecosystem and a need to better understand the consequences of anthropogenic pressures to these services.

Introduction

Bioturbation refers to an animal directly or indirectly affecting the physical, chemical and biological properties of their habitat (Kristensen et al., 2012). Occurring in terrestrial (Eldridge & Mensinga, 2007; Ernst et al., 2009; Valentine et al., 2013) and aquatic (Oliver & Slattery, 1985; Ray et al., 2006; Williamson et al., 2021) ecosystems, species alter their habitat through feeding, digging and burrowing. Due to these processes, bioturbating species are classified as ecosystem engineers, modifying their habitat physically, changing the dynamics of ecosystems, and regulating the availability of resources for other species (Jones et al., 1994; Kristensen et al., 2012; Wright & Jones, 2006). As a result, many of these species are keystone species that disproportionately contribute to ecosystem services within their environment.

Rays (superorder Batoidea: skates, stingrays, electric rays and shovelnose rays) are key bioturbators and ecosystem engineers. Rays forage for buried prey by excavating sediment via hydraulic action, expelling water from their mouth and gills while simultaneously moving their pectoral fins, creating feeding pits that can remain on the substrate for several days (Gregory et al., 1979; Sasko et al., 2006; Wilga et al., 2012). Rays play a significant role in foraging facilitation (Boaden & Kingsford, 2012), stabilizing local prey populations (Ajemian et al., 2012; Hines et al., 1997) and influencing prey metapopulation source-sink dynamics (Peterson et al., 2001). In addition, ray bioturbation contributes to many ecological services within their ecosystem through the structuring of sediments, oxygen penetration and nutrient cycling (Harris et al., 2016; Lohrer et al., 2004; Wright & Jones, 2006). Rates of ray bioturbation can vary spatially and temporally and are influenced by physical characteristics of the seafloor, ray size, shape and behaviour (Crook et al., 2022; Myrick & Flessa, 1996). The distribution and density of prey are also influencing factors (Takeuchi & Tamaki, 2014) as animals forage in energetically rewarding areas with high prey densities (Charnov, 1976) resulting in spatially clustered patterns, which has been documented in other studies (Hines et al., 1997; Takeuchi & Tamaki, 2014). These patterns can also be a result of ray size, as larger rays can forage deeper into the sediment (Takeuchi & Tamaki, 2014), theoretically resulting in clusters of size-dependent feeding pits. These considerations highlight that turnover rates and volume of sediment displaced over time are likely to be habitat and species-specific and prey density dependent, varying the ecosystem services contributed by ray bioturbation.

Assessing the amount of sediment turned over by ray foraging activities provides vital information for quantifying the ecosystem services delivered by these marine organisms. Previous studies have quantified the abundance, volume and turnover of ray feeding pits for different species at different sites. For example, on an intertidal sandflat in Japan, Red Stingrays (Hemitrygon akajei) created 18103 feeding pits daily in a 1 km² area (Takeuchi & Tamaki, 2014), in contrast to a sandflat in New Zealand where Eagle Rays (Myliobatis tenuicaudatus) created just 21 new feeding pits daily over a similar area (Hines et al., 1997). Rays from the Dasyatidae family were estimated to turnover 2.42% of the area of 1.5 km² intertidal flats on Ningaloo reef in 7 days (O'Shea et al., 2011). A single Australian whipray (Himantura australis) is capable of excavating 200 tonnes of sediment in a single year (Crook et al., 2022), with other species completely turning over a 60 m² plot in 102 days (Atlantic Stingrays Dasyatis sabina) (Reidenauer & Thistle, 1981). While previous studies have mainly concentrated on quantifying ray-driven sediment turnover, they have not definitively established its significance to the ecosystem (Flowers et al., 2021). Understanding the ecological importance of rays is crucial to comprehend what population declines and loss of their ecosystem services mean for the wider ecosystem (Heithaus et al., 2010).

Understanding how ray bioturbation shapes aquatic ecosystems requires a quantitative spatial analysis of their bioturbation rates. This study assessed the volume and mass of sediment displaced over a 7-day period by rays during foraging activities in a known intertidal feeding area and extrapolated those patterns to an estuary scale in intertidal areas. We expected rays to displace high amounts of sediment on a daily and weekly basis, contributing to the bioturbation in the system. We predict to see spatial clustering of the feeding pits which will hypothetically reflect prey densities in the area. The methods developed here, and the results from this study provide insight into the importance of ray bioturbation in a temperate estuarine system.

Materials and Methods

Study site

Pelican's Spot is a small (9223 m²) low-lying island in the Brisbane Water estuary east of Woy Woy, NSW, Australia (Fig. 1). With a tidal range of 1.26–1.42 m, this area has abundant seagrass meadows (Zostera spp.) and sandy intertidal flats that Estuary Stingrays (Hemitrygon fluviorum) use for feeding and resting all year round (Bourke et al., 2023). There is limited data on what other ray species are present within the estuary. However, species from the families Trygonorrhinidae, Rhinobatidae, Dasyatidae, Myliobatidae, Rhinopteridae and Urolophidae are abundant in similar temperate estuaries along the coast of NSW and are likely to also inhabit this estuary. Despite H. fluviorum being the only confirmed species residing in the estuary, this species represents other ray species with similar morphologies and ecological roles as bioturbators implying that our findings may be applicable elsewhere (Oleksyn et al., 2020).

[IMAGE OMITTED. SEE PDF]

Drone surveys

Drone missions were conducted daily for seven consecutive days in May 2023. Surveys were flown at sunrise, during low tide and with minimal wind (<10 km h⁻¹) to reduce sun glint, water column effects and improve image quality (Joyce et al., 2018). The drone used for this study was an Autel Robotics Evo 2 Pro equipped with a 1-inch sensor camera that can capture 20-megapixel still photos. Controlled through the Autel Robotics V2 app, an autonomous oblique mapping mission was conducted at 20 m altitude which resulted in a ground sampling distance (GSD) of 0.45 cm px⁻¹. The drone was flown with the camera at NADIR (pointing straight down) across the mission area as well as with the camera at 45° angles to collect the oblique images, both with a front and side overlap of 80% to effectively capture the depth of ray feeding pits. This resulted in a polygon mission with numerous transects that had overlapping imagery. Missions took 80 min of flight time to complete, producing over 1400 images, and covering an area of 1443 m².

Field calibration

Since drone-derived digital elevation models can have relatively high measurement error in the z axis (Koci et al., 2017), prior to the week-long survey, a calibration survey was conducted to determine the measurement error between the feeding pit depth collected from drone imagery (hereby known as the drone-sourced depth) and the in situ, manually measured depths (hereby known as ground-measured depths). Ground control points consisting of yellow set squares were positioned next to 20 feeding pits where depth and width at the widest point were measured manually with a tape measure to the nearest centimetre. The drone then surveyed the area, and the ground control points were used to match the ground-measured depths to the correct pit in the drone imagery and extract the drone-sourced depth measurements. The measurement error applied to the drone-sourced depths was calculated by dividing the ground-measured depth by the drone-sourced depth, producing the mean and standard deviation of these measurement errors.

During the calibration survey, ten 100 mL syringe sediment cores were taken from the sediment in and around the feeding pits within the survey area. These sediment cores were used to approximate overall sediment densities in the area, which could then be used to relate calculated feeding pit volumes to sediment mass.

Data processing

Agisoft Metashape Professional version 1.8.4 was used to process the drone images on a desktop PC with 11th Gen Intel®Core™ [email protected], a GeForce RTX 3080 and 64GB RAM. This process is summarized in Table 1. Orthomosaics and Digital Elevation Models (DEM) from each day were imported into ArcGIS Pro v3.0.0 for further analysis. Models from days 2 to 7 were georeferenced in ArcGIS Pro by matching identifiable features in each model to the same feature in the day 1 model. A point was placed in the centre of all visible feeding pits in the day 1 model. The georeferenced orthomosaics from days 2 to 7 were then overlaid onto one another to identify historical and new feeding pits over a week period. Historical and new feeding pits were marked with a point shapefile placed in the centre of each pit.

Table 1 Data processing steps in Agisoft Metashape Professional.

Step	Process	Parameter	Processing Time
1	Align photos	Highest	50 minutes
2	Trim photos to remove outliers		N/A
3	Dense Cloud and Depth Maps	Medium quality and aggressive filtering	70 minutes
4	Digital Elevation Model	Created from dense cloud	5 seconds
5	Orthomosaic Model	Created from DEM	20 minutes

To calculate the volume of each pit, we delineated the area of each pit and obtained its depth. The outline of each pit was defined in the daily orthomosaic using unsupervised object-based image classification (Table 2). This created a raster layer distinguishing the pits from the surrounding substrate (Fig. 2). The layer was converted to a polygon and projected to WGS84 UTM56S. To extract only the relevant daily pit features (rather than existing pits or rocks/other features), the polygon layer was intersected with the previously digitized daily pits point layer. Zonal statistics was applied to the DEM within the area of each pit polygon, giving the maximum and minimum depths (depth range) within each pit. See workflow of image classification in Table 2 and Figure 2.

Table 2 Image classification steps in ArcGIS Pro.

Step	Process	Parameters	ArcGIS Pro tool
1	Image classification	Unsupervised classification Object-based Default classification schema Spectral detail = 10 Minimum segment size = 40 px	Classification wizard
2	Convert raster to polygon		Raster to polygon
3	Project polygon	WGS84 UTM56S	Project tool
4	Intersect projected polygon layer with the layer of existing pits	Select and delete background
5	Zonal statistics	Input = projected polygon Zone field = ObjectID Input value raster = DEM	Zonal statistics

[IMAGE OMITTED. SEE PDF]

Spatial autocorrelation analysis

All spatial autocorrelation assessment was conducted within ArcGIS Pro (version 3.1.3). To determine if the pits tended to be clustered together (i.e. not randomly distributed over the study area), we ran an Average Nearest Neighbour Analysis of the entire dataset of pits (days 2–7 merged) within the study area (1443 m²). In addition, we tested if the size of the pits (area) was spatially autocorrelated using a Global Moran's I test. To identify where clusters of pits with higher or lower areas than randomly expected were located a Cluster and Outlier Analysis (Local Moran's I) and Hot Spot Analysis (Getis Ord Gi*) were used with false discovery rates applied. Local Moran's I and Getis Ord Gi* tests were both run using inverse distance weighting (Euclidean distance) with a threshold distance of 7.5 m to conceptualize the spatial relationships (i.e. define the neighbourhood). This threshold distance was selected by determining the distance required for each pit, in each day's dataset, to have at least one neighbour (a requirement of the tests). The Calculate Distance Band from Neighbour Count tool was used for this. Multiple other conceptualizations of spatial relationships were also tried, and all produced very similar results lending confidence to our findings.

Field calibration analysis

To determine dry sediment density, the wet and dry weight of sediment cores (n = 10) was measured to obtain the ratio of dry sediment weight. This was calculated to give a mean ± SD (1.45 ± 0.12 g cm⁻³).

To obtain a calibrated depth measurement, the measurement error (3.61 ± 2.28 cm, mean ± SD) obtained from the calibration survey was used to convert drone-sourced depth measurements to calibrated depths. To account for the uncertainty in the measurement error, the drone-sourced depths were multiplied by the measurement error uncertainty (2.28), giving a calibrated uncertainty depth. Pits were standardized by removing depths outside of 2–41 cm following pit depth observations in Crook et al. (2022), to ensure outliers were removed. This removed approximately 80 feeding pits from the 1170 detected.

Data analysis

To measure the volume of the pits, we applied the formula for a cone, assuming feeding pits will follow a cone shape:1$\text{Cone volume}\left({\mathrm{m}}^3\right)=\mathrm{SA}\times \left(\frac{D}{3}\right)$where SA = surface area of the polygon (m²), D = the calibrated depth (m).

This formula was also applied to the correction uncertainty depth with D referring to the calibrated uncertainty depth (m).

To obtain the mass of the sediment in each pit we applied the formula:2$\begin{array}{c}\text{Mass}\left(\mathrm{g}\right)=\text{cone volume}\times \text{mean}\mathrm{dry}\text{density}\end{array}$where Cone volume = volume of each feeding pit (cm³), Mean dry density = the mean dry density of the dried sediment (g cm^-3).

To account for the uncertainties found within the depth correction factor (3.61 ± 2.28) and the mean dry density (1.45 ± 0.12), a propagation uncertainty formula was applied to each individual mass to obtain the standard deviation of both uncertainties.

If:3$\begin{array}{c}Z=X\times Y\end{array}$Then:4$\begin{array}{c}\Delta Z=Z\times \sqrt{{\left(\frac{\Delta X}{X}\right)}^2+{\left(\frac{\Delta Y}{\mathit{XY}}\right)}^2}\end{array}$where Z = mass (g), ΔZ = standard deviation of the mass (g), X = mean dry sediment density (g cm^-3), ΔX = standard deviation of the mean dry sediment density (g cm^-3), Y = individual cone volume (cm³), ΔY = cone volume with the correction uncertainty depth applied (cm³).

We now have the mass of sediment displaced from each pit (g) and the standard deviation mass (g) to account for both uncertainties in the depth correction factor and mean sediment density. The total mass of sediment displaced for pits each day was calculated. To calculate the total standard deviation mass we used the formula:5$\begin{array}{l}\text{Total standard deviation mass}\left(\mathrm{g}\right)\\ =\sqrt{\left(\mathrm{SD}{1}^2+\mathrm{SD}{2}^2+\cdots +{\text{SDth}}^2\right)}\end{array}$where SD = individual standard deviation mass of each pit (g).

To then find the average standard deviation mass per day, the following formula was used:6$\begin{array}{l}\text{Average standard deviation}\\ =\sqrt{\frac{\mathrm{SD}{2}^2+\mathrm{SD}{3}^2+\dots +\mathrm{SD}{7}^2}{n}}\end{array}$where SD = the total standard deviation mass for Days 2–7 (g), n = total number of days.

To determine any significant differences between the mass of sediment displaced each day, an Analysis of Variance (ANOVA) with a log transformation and Tukey Post Hoc Test was conducted in R Studio version 4.2.3 (RStudio Team, 2023).

Data extrapolation

To extrapolate measured bioturbation rates to the estuary scale, high-resolution aerial imagery (GSD; 5.5–7.5 cm px⁻¹) was obtained from Nearmap imagery of the Brisbane Water estuary and imported into ArcGIS Pro. Aerial imagery from the 5 October 2021 was used due to high water clarity. Intertidal areas with visible ray feeding pits in bare sediment were identified and digitized, and the total ray bioturbation area within the visible intertidal areas of the estuary was calculated. To standardize this process, ‘feeding areas’ met expertly determined criteria: feeding pits had to be visible (not under jetties or covered by trees), seagrasses were excluded due to difficulty identifying feeding pits, and feeding pits were distinct from other depressions in the sediment. Since these criteria exclude large areas likely to contain ray feeding and bioturbation (e.g. seagrasses), they would lead to a conservative measurement of feeding area across the estuary.

To extrapolate the mean mass and the mean standard deviation mass per day to the estuary scale, the following formula was used:7$\begin{array}{l}\text{Mass}\mathrm{per}\mathrm{day}\text{across the estuary feeding area}\\ =\frac{K}{\text{estuary area}}\end{array}$where K (tonnes m⁻²) = mass per day (average mass or average standard deviation mass) × study area, Study area = 1443 m², Estuary feeding area = 69287.9 m².

This daily rate was then extrapolated further over a year scale.

Results

A total of 1090 feeding pits were measured over a 7-day period in the 1443 m² study area (Fig. 3). The number of newly formed pits each day ranged from 74 to 329 with no recognizable trend (Table 3). The cumulative surface area covered by these pits was 90.41 m² (6.2% of the study area) with an average feeding pit surface area of 0.095 m² (±0.01 m²) and an average calibrated depth of 0.11 m (±0.03 m). The total volume of sediment displaced during ray feeding activities over the week was 4.95 m³ (Table 3), with a mean volume of 0.83 m³ (±0.52 m³) sediment displaced each day. This is approximately 575.2 cm³ per square metre per day. Incorporating the measured sediment density this equates to 7.19 (±3.45) tonnes of sediment displaced during feeding activities in the study area during the 7 days (Fig. 3; Table 3). On average, this is 1.20 (±0.68) tonnes of sediment each day. There was a significant difference in the log mass of sediment excavated each day (ANOVA, d.f. = 5, F = 10.72, P < 0.05; Fig. 4). Tukey's HSD pairwise tests suggest the log mass of sediment excavated on Day 2 was significantly different to Days 6 and 7, and Day 4 was significantly different to Days 3, 5, 6 and 7.

[IMAGE OMITTED. SEE PDF]

Table 3 Summary of ray feeding pits.

Day	Number of new pits	Mean pit area ± SD (m2)	Mean pit depth ± SD (m)	Total volume ± SD (m3)	Total mass ± SD (tonnes)
2	170	0.10 ± 0.12	0.13 ± 0.07	0.99 ± 0.62	1.44 ± 0.69
3	89	0.08 ± 0.04	0.08 ± 0.04	0.18 ± 0.11	0.26 ± 0.13
4	257	0.10 ± 0.09	0.16 ± 0.08	1.37 ± 0.68	1.98 ± 0.95
5	329	0.11 ± 0.13	0.10 ± 0.08	1.52 ± 0.96	2.21 ± 1.06
6	74	0.08 ± 0.09	0.09 ± 0.06	0.29 ± 0.18	0.42 ± 0.20
7	171	0.10 ± 0.11	0.08 ± 0.07	0.61 ± 0.38	0.88 ± 0.42
Total	1090			4.95 ± 3.13	7.19 ± 3.45

[IMAGE OMITTED. SEE PDF]

Total ray feeding habitat determined through visual identification of feeding pits from aerial imagery was 69288 m² in Brisbane Water estuary (Fig. 5), with Pelican's Spot covering 2% of the total area where ray feeding activity is observed. When extrapolating bioturbation rates collected in this study to an estuary scale, rays likely displace 57.57 (±32.44) tonnes of sediment per day, resulting in a conservative estimate of 21012.15 (±11840.28) tonnes of sediment displaced per year assuming constant rates between seasons.

[IMAGE OMITTED. SEE PDF]

The Nearest Neighbour Ratio of 0.8 (z = −11.63, P < 0.001) showed the location of pits to be significantly clustered. The Global Moran's I test revealed significant but weak positive spatial autocorrelation (P < 0.00) of feeding pit size for all days except day three (Table 4). The strongest positive spatial autocorrelation occurred on day two (Moran's I = 0.33, z = 11.18, P = 0.00). Hot spot analysis showed statistically significantly larger feeding pits were clustered in three areas (Fig. 6). A cluster of statistically significantly smaller pits was found at the eastern edge of the study area. Cluster, Outlier and Hot Spot Analysis showed similar results (refer to Appendix S1 to see the hot spot analysis results for each day's pits).

Table 4 Moran's I values assessing spatial autocorrelation in ray pit area throughout the study.

Days	Moran's I	P-value	z-value
2	0.33	0.00	11.18
3	0.07	0.15	1.43
4	0.17	0.00	7.12
5	0.24	0.00	12.62
6	0.10	0.05	1.98
7	0.11	0.00	3.57
All combined	0.15	0.00	27.20

[IMAGE OMITTED. SEE PDF]

Discussion

This study provides an account of ray bioturbation in a temperate, wave-dominated estuary. We calculate that rays in Brisbane Water estuary can displace approximately 20000 tonnes per year. The Estuary Stingray is the most dominant species in the estuary and likely responsible for a significant portion of bioturbation activity. It is also near-threatened with declining populations (Rigby & Derrick, 2021), making it likely that observed bioturbation rates are already lower than those prior to human colonization. This highlights the important role of rays and they have likely shaped many estuarine sedimentary habitats.

The considerable mass of sediment being displaced outlines how critical these rays are to the ecosystem services within their estuary. Numerous studies have assessed the significant effect of ray bioturbation on coral reefs (O'Shea et al., 2011) and other soft bottom environments worldwide (Hines et al., 1997; Reidenauer & Thistle, 1981; Takeuchi & Tamaki, 2014). Studies have suggested that rays may have a significant effect on infaunal community structure and stabilization (Hines et al., 1997; Thrush et al., 1991; VanBlaricom, 1982), and nutrient dynamics (D'Andrea et al., 2004; VanBlaricom, 1982), but further research is required on their role in nutrient dynamics (Flowers et al., 2021). A combination of foraging and high ray abundance may not only be important for bioturbation but may physically alter the hydrodynamics and topography of the seafloor, with some species seen to facilitate the erosion of sandbanks (H. australis and P. ater) (Crook et al., 2022). Rays also form mounds around feeding pits, creating rugosity within the seafloor (Crook et al., 2022). By altering intertidal topography and increasing suspended sediment rays can affect the geomorphology of an area influencing the movement of tidal currents (Grant & Madsen, 1979).

Rays in temperate estuaries face many anthropogenic pressures such as climate change, habitat degradation, commercial and artisanal fishing (Cavanagh et al., 2003; Dulvy et al., 2014; Stevens et al., 2005). These activities have resulted in population declines with as many as 19.9% of ray species listed in a threatened IUCN red list category (Dulvy et al., 2014; Last et al., 2016), which will have consequences for the ecological roles rays provide (Heithaus et al., 2010). In the context of rays as ecosystem engineers, a decline in population will have flow-on effects to ecosystem services. Despite finding 1090 feeding pits within the study area over 7 days, it is difficult to estimate ray population size in the area from this. Previously, the number and size of feeding pits was not linked to the abundance and demographics of rays in the Brisbane Water estuary, with feeding pits underestimating the size of rays found (Bourke, 2021). However, it is expected that bioturbation rates will decline linearly with ray abundance, as well as a decrease in the connected ecosystem services. If ray numbers continue to decline globally, this large bioturbation service within this ecosystem would be lost, underlining why additional efforts must be made to protect these animals.

Observed ray bioturbation rates were comparable to other studies. However, some studies have found low sediment turnover <150 cm³ m⁻² per day (D'Andrea et al., 2004; Hines et al., 1997; O'Shea et al., 2011), while others found that rays turn over much larger volumes of sediment per day >1000 cm³ m⁻² per day (Myrick & Flessa, 1996; Takeuchi & Tamaki, 2014; Thrush et al., 1991). Our study fits in the higher end of this range, recording a turnover of 575.2 cm³ m⁻² per day (Table 5). Grant (1983) found rays reworked 6.8–24.2 cm³ m⁻² (h ebb flow)⁻¹ of an intertidal sandflat in South Carolina, USA, significantly less than physical processes which reworked 6100–12000 cm³ m⁻² (h ebb flow)⁻¹. As such, the impact of biological processes to an ecosystem is thought to be background noise compared to the impact of the physical processes. Grant (1983)'s study showed that rays displaced just 10 kg of sediment over two summers in a 2000 m² area. In contrast, our study calculated that rays moved 1.2 tonnes of sediment in a day in a similar area. This variability in volume of sediment displacement is likely attributable to differences in species' foraging ecology, morphology, environmental conditions and prey population densities. Crook et al. (2022) discovered that two co-occurring species of rays, the Australian Whipray (Himantura australis) and the Cowtail Stingray (Pastinachus ater), turned over different volumes of sediment due to differences in feeding behaviours, despite being co-occurring and predating on the same prey. H. australis used more excavation type methods to forage prey, displacing more sediment than P. ater which used nondisruptive foraging methods. Differences within studies can also be attributed to the methodologies used in each study such as manual pit measurements that can disturb the structure of pits, aerial photography and direct-infilling (Flowers et al., 2021). While it was not possible to determine which species drove the bioturbation in our study site, our results corroborate evidence that their contribution to bioturbation is substantial in most cases (Table 5).

Table 5 Previous studies on ray sediment turnover on intertidal and subtidal habitats.

Study	Sediment turnover (cm3 m−2 per day)	Study location	Ray species
Myrick and Flessa (1996)	2883.3	Gulf of California, Mexico	Bat ray (Myliobatis californica) Round Stingray (Urolophus helleri)
Takeuchi and Tamaki (2014)	1429.3	Kyushu, Japan	Red Stingray (Dasyatis akajei)
Thrush et al. (1991)	1108.2	Manukau Harbour, New Zealand	New Zealand Eagle ray (Myliobatis tenuicaudatus)
Reidenauer and Thistle (1981)	611.1	Florida, U.S.A	Atlantic Stingray (Hypanus sabinus)
Our study	575.2	New South Wales, Australia	Most likely Estuary Stingray (Hemitrygon fluviorum)
VanBlaricom (1982)	420	California, USA	Round ray (Urolophus helleri) Bat ray (Myliobatis californica)
Hines et al. (1997)	180.8	Manukau Harbour, New Zealand	New Zealand Eagle ray (Myliobatis tenuicaudatus)
D'Andrea et al. (2004)	145	South Carolina, U.S.A	Stingray (Dasyatidae)
O'Shea et al. (2011)	34.5	Ningaloo Reef, Western Australia	Stingrays (Dasyatidae)

Other marine taxa change the biological and physical processes within their ecosystem through bioturbation. For example, many species displace sediment when foraging for benthic prey (Nelson et al., 1987; Shimek, 1977; Virnstein, 1977), with other taxa such as invertebrates displacing sediment through burrowing. The most notable invertebrate are callianassid shrimp, estimated to displace an average of 0.56 m³ m⁻² per year on a sand flat in Mexico (Myrick & Flessa, 1996), almost double the sediment displaced in the present study. For a small invertebrate, this is a significant displacement of sediment, with this species being documented to erase the physical effects of cyclones on sediment habitat within 6 weeks (Riddle, 1988). This indicates that the contribution of bioturbating marine species to ecosystem services can be substantial.

Spatial clustering of feeding pits is hypothesized to be attributable to the density-dependent foraging behaviours of rays (Hines et al., 1997). A reflection of foraging theory where animals will maximize the cost versus reward of foraging for example time to forage over energy obtained from the prey item (average-rate-maximizing theory) (Stephens & Krebs, 1986). This suggests that predators will distribute themselves across varying patches of prey densities consistent with prey availability (Fretwell & Calver, 1969; Sutherland, 1983). Eagle rays in New Zealand displayed a nonlinear response to the densities of their prey the marine clam Macoma liliana (Family Tellinidae), with low ray foraging activity seen in areas with low M. liliana densities (Hines et al., 1997). The study indicated that ray foraging rates increased significantly with prey density, and Red Stingrays showed similar patterns of prey density selectivity (Takeuchi & Tamaki, 2014). Despite not examining prey densities at our site, we can infer prey densities from the spatial clustering of the ray feeding pits in conjunction with foraging theory models. The spatial autocorrelation analysis indicated a positive relationship between feeding pit size and clustering. This could be driven by differences in prey depths, where larger rays are able to access deeper, more rewarding prey than smaller rays (Ebert & Cowley, 2003; Tillett et al., 2008). Future studies should aim to determine whether prey depth affects bioturbation depth, and how that relates to foraging ray size.

Bioturbation rates presented in this study are based on measurements with uncertainties and assumptions, providing a conservative analysis of the sediment displaced during ray feeding activities. Although our drone missions covered a large area, they only represent 2% of the total feeding area identified within the estuary, and it is possible that bioturbation rates differ in the other identified bioturbation zones. While we incorporated multiple measurement errors to account for the uncertainty in the drone derived measurements, we can only provide an estimation of the feeding pit depths. An incorporation of an RTK (real-time kinematics) GPS equipped drone and base station would reduce z measurement errors associated with the depth measurements due to its real-time processing of GNSS signals and ability to accurately and precisely capture elevation (Nota et al., 2022). However, this study provides a conservative estimate of ray bioturbation in one area and extrapolated that to an estuary scale using aerial imagery. It is also difficult for aerial imagery to capture all areas where feeding pits would be (e.g. under jetties and in deeper water), and bioturbation rates were only assessed in one estuary in one season. Different estuaries will have different physical characteristics, prey populations and different ray species with different feeding behaviours and morphologies. While we give an estimate of bioturbation within the Brisbane Water estuary, we can only imply similar results in Australian estuaries. Daily variability in the number of feeding pits could be attributed to the extent of the feeding areas across the estuary. Foraging activity was measured over the 7 days in a small area of the estuary; however, rays can potentially move over larger spatial scales than the feeding area examined here. There is limited data on the movement ecology and site fidelity of H. fluviorum (Constance et al., 2023); however, we assume rays move throughout the estuary to use a range of feeding areas as evident by the detected feeding areas in Figure 5. The movement in and out of these areas could contribute to the variation in the daily feeding activity in the area. If this species does display site fidelity, abiotic factors such as temperature and tides could also explain the variation (Oleksyn et al. 2020; Elston et al., 2022). Annual variability in temperature could influence the annual feeding rates, as rays are ectotherms whose metabolisms and feeding rates will correlate with ambient temperatures (Elston et al., 2022). An assessment of seasonal variation in feeding rates is an important avenue for further future research. Despite these technical considerations, the focus of this study was to conservatively quantify the sediment displaced by rays during feeding and provide an estimate on an estuary scale, instead of assessing the processes that influence ray bioturbation.

Overall, our findings indicate that the rays in the Brisbane Water estuary displace a large volume of sediment, contributing to the biological and physical ecosystem services within their estuary. This highlights the ecological importance of rays and a need to better understand the consequences of anthropogenic pressures to these services. These findings underline the need for more in depth research on how rays use Australian temperate estuaries, providing further knowledge on how anthropogenic impacts alter ray diet, movement and habitat use, aiding their conservation within estuarine systems.

Acknowledgements

The authors would like to thank Isaac Hampson for his assistance in the calibration survey. Open access publishing facilitated by The University of Newcastle, as part of the Wiley - The University of Newcastle agreement via the Council of Australian University Librarians.

Author contributions

Molly Grew: Conceptualization, Methodology, Formal analysis, Investigation, Writing – Original Draft, Review & Editing Troy F. Gaston: Funding acquisition, Conceptualization, Methodology, Validation, Resources, Writing – Review & Editing, Supervision Andrea S. Griffin: Conceptualization, Methodology, Validation, Formal analysis, Resources, Writing – Review & Editing, Supervision Stephanie J. Duce: Conceptualization, Methodology, Validation, Formal analysis, Resources, Writing – Review & Editing Vincent Raoult: Conceptualization, Methodology, Validation, Formal analysis, Writing – Review & Editing, Supervision.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Statement on inclusion

Our study brings together authors from a range of different disciplines from spatial science, animal behaviour and marine ecology, broadening the perspective of the manuscript. The author list includes scientists from the local area as well as outside of the state and includes a diverse gender author list with three of the five authors being women in STEM.

Data availability statement

The data will be made available via Supplementary Information.

References

Ajemian, M.J., Powers, S.P. & Murdoch, T.J. (2012) Estimating the potential impacts of large mesopredators on benthic resources: integrative assessment of spotted eagle ray foraging ecology in Bermuda. PLoS One, 7, [eLocator: e40227]. Available from: [DOI: https://dx.doi.org/10.1371/journal.pone.0040227]