Specimens of Diopatra and Arenicola complexes, as well as sediment, were collected in Ria de Aveiro, a shallow, mesotidal costal lagoon located on the northwestern coast of Portugal, between 40° 30`N and 40° 51`N (Vaz et al., 2005), with a longitudinal gradient of salinity from about 0 to 36. Extensive intertidal mud and sand flats, salt marshes and islands can be found along its main channels (Costa et al., 2022). Given the large number of specimens needed and the effort associated with their capture, Diopatra and Arenicola were obtained from a fishing store that receives the specimens collected by local fishermen in Ria de Aveiro. Sediment was collected in Mira Channel, one of the southern channels, where the presence of Arenicola has been reported (Pires et al., 2015; Costa et al., 2022). Sediment was homogenized and macrofauna removed by hand-piking for use in housing and microcosms experiments. Specimens were transported to the laboratory and placed in containers (370 x 270 x 230 mm, L x W x H) with 120 mm of sediment overlain with 100 mm of artificial seawater (salinity 28, prepared by mixing Red Sea Coral Pro Salt (Germany) and reverse osmosis water). The salinity level was established based on the distribution areas of both species in Ria de Aveiro (Rodrigues et al., 2011; Vargas et al., 2017). The specimens passed by a housing period of 12 days to ensure the complete regeneration of the posterior part and to ensure that all the individuals were healthy and feeding (Pires et al., 2012). During the housing period, continuous aeration was provided using air-stones and the containers were supplied with fish food (Supreme 50, Skretting) every day by placing a small grain next to the entrance of each tube. If the food was not dragged into the tube within 5 minutes, it was removed. Food was not provided on the day before the polychaete specimens were transferred to the microcosms. During the housing period, the average water temperature in the containers was of 17.9 ± 1.0°C, pH 8.2 ± 0.3, oxygen concentration 7.3 ± 0.8 mg L−1, and salinity 28.0 ± 0.1 (mean ± s.d.).

A microcosm set-up was designed to simulate different degrees of Diopatra (Dp) habitat invasion by Arenicola (Ar). Sediment was distributed along 30 vertical glass aquaria (135 × 135 × 500 mm, internal dimensions), creating a 150 mm column (~ 1/3 of the height of the aquarium, 2.7 L of sediment) and filled with artificial seawater (~ 2/3 of the height of the aquarium, 5.5 L of water). The detailed timeline of the microcosms’ experiment, related to water replacement, housing period, acclimation, luminophores introduction, and monitoring of bioturbation and water column nutrients’ concentration is provided in Supplementary Table S1 . Each aquarium was gently aerated by releasing air bubbles through two capillary tubes (ø = 0.84 mm), connected by a silicon tube (ø = 4.00 mm) to an air pump (RESUN LP-100, 140 L min^-1). Aquaria were randomly distributed across the bench and exposed to natural light conditions with a 10 h light:14 h dark photoperiod. Five different ratio of Diopatra: Arenicola, with an equivalent constant density of ~180 ind.m^-2, plus Control (no organisms), were placed in the aquaria to be tested: 4 Diopatra (coded as 4Dp:0Ar), 3 Diopatra and 1 Arenicola (3Dp:1Ar), 2 Diopatra and 2 Arenicola (2Dp:2Ar), 1 Diopatra and 3 Arenicola (1Dp:3Ar), and 4 Arenicola (0Dp:4Ar) (5 replicates each, totaling 30 aquaria). Adult Diopatra individuals were introduced into the experimental units without their original tubes. However, by Day 0, prior to the introduction of luminophores and the start of the experiment, all individuals had already built new tubes (confirmed by direct observation). These newly constructed tubes were functional and appeared structurally stable throughout the experimental period. To minimize evaporation and prevent pronounced salinity shifts, all aquaria were covered with Parafilm^® Sealing Film. The water in each aquarium was not renewed and neither sediment nor food were added during the acclimation and experimental period.

To avoid additional manipulative stress before implementing the experiment, the polychaetes were measured (total length) and weighted only at the end of the experiment. Organisms passed by a depuration period of 24 h to clear their gut before weighing. The average length and weight were 74.08 ± 10.63 mm and 2.12 ± 0.46 g for Diopatra, respectively; and 86.62 ± 22.23 mm and 3.78 ± 0.75 g for Arenicola (mean ± s.d.), respectively.

Particle reworking was measured non-invasively using fluorescent sediment profile imaging (f-SPI; Solan et al., 2004; Lopes et al., 2018; Crespo et al., 2021), based on the visual detection of fluorescent luminophores on the sides of the aquaria. This profile imaging provides descriptors for faunal-mediated sediment reworking, which are surrogate measures for different aspects of infaunal organisms’ behavior (Teal et al., 2010; Hale et al., 2014). To avoid the strong initial activity associated with polychaetes establishment and burrowing, luminophores (40 g, 125 - 250 µm diameter, yellow color) were added on the top of the sediment, near the side of each aquarium which would be photographed later, at least 24 h after the introduction of organisms (D₀). This method allows the deposition of a 1–2 cm layer of luminophores on the side that will be photographed and ensures the supply of luminophores throughout the entire experimental period.

The aquaria were photographed within a dark box, under Ultraviolet light, using a GoPro HERO6 camera (12 megapixels sensor, set for 1/8 seconds exposure, diaphragm aperture of f/2.8, and film speed (light sensitivity) equivalent to ISO 3199). Each image obtained was subsequently cropped according to the internal width of the aquarium (135 mm = 2050 pixels, effective resolution = 63.4 µm px^-1), converted to a red-green-blue (RGB) stack and saved with JPEG compression (Joint Photographic Experts Group). Images were then analyzed using a custom-made plugin that runs within ImageJ (Version 1.48c), a java-based public domain program developed at the US National Institutes of Health (available at http://rsb.info.nih.gov/ij/index.html). The plugin output provides details on the mean (^f-SPIL_mean, time-dependent indication of mixing), median (^f-SPIL_med, typical short-term depth of mixing), and maximum (^f-SPIL_max, maximum extent of mixing over the long-term) mixed depth of particle redistribution (Hale et al., 2014), for each aquarium, after the semi-automatic conversion of the images to a binary data matrix of the distribution luminophore pixels (0 = background sediment, 1 = luminophore). For each image, the vertical distance between the highest and lowest points of the water sediment boundary (surface boundary roughness – SBR) was manually assessed. SBR is a proxy for surficial activity and is particularly useful in the case of species with high accretion abilities, such as Arenicola, which creates noticeable faecal casts on the surface of the sediment (Volkenborn et al., 2007). As bioturbation processes are cumulative, only the final day (day 21) image data were considered, and the values found for the tested species ratio were considered against those found for the procedural control.

Before the experiment inception, five sub-samples were analyzed to determine grain-size and organic matter content on the homogenized sediment. Sediment grain-size was analyzed by dry sieving standard method (Sutherland, 1998). Median and percentage of fines were used to classify the sediment, according to the Wentworth scale (Blott and Pye, 2001): very fine sand (median from 3 - 4 Φ; 63 - 125 µm); fine sand (2 - 3 Φ; 125 - 250 µm); medium sand (1 - 2 Φ; 250 - 500 µm); coarse sand (0 - 1 Φ; 500 - 1000 µm); very coarse sand (-1 - 0 Φ; 1000 - 2000 µm). The final classification adopted the description ‘clean’, ‘silty’ or ‘very silty’ when the % fines fraction of the total sediment, by dry weight, ranged from 0% to 5%, 5% to 25% and 25% to 50%, respectively (Blott and Pye, 2001). Samples with >50% fines content were classified as mud. Total organic matter content was obtained by loss on ignition as described by Kristensen and Andersen (1987).

At days 0, 2, 5, 8, 12, 16, and 21, 10 ml water aliquots (filtered through Whatman GF/C glass-fiber filters) were sampled on the top of the water column (0–5 cm), from each aquarium: this low volume ensured a loss lower than 10% on the microcosm total volume. Samples were stored at −20°C until analysis of dissolved inorganic nutrients (ammonium, NH₄-N; oxidized form of dissolved inorganic nitrogen, NO_x-N; phosphate, PO₄-P). The determination of the concentrations of NH₄-N and PO₄-P was performed following standard methods described in Limnologisk Metodik (Københavns Universitet, 1992). The concentrations of NO_x-N were determined using a flow injection system (FIAstar 5000 Analyzer, Höganäs, Sweden), following Strickland and Parsons (1972). Analytical quality control was ensured by calibration curves derived from standard solutions, run at the beginning of the analysis and in parallel with blanks and samples. Water temperature, pH, concentration of dissolved oxygen and salinity were measured in the mesocosm on every sampling day using a WTW – pH 330i/set equipped with SenTix^® 41; a WTW – cond 3110/set 1 equipped with TetraCon^® 325 and a WTW – Oxi 3210/set 2 equipped with CellOx^® 325-3.

The dependent variables related to sediment reworking descriptors (SBR, ^f-SPIL_mean, ^f-SPIL_med, and ^f-SPIL_max) and dissolved nutrient concentrations ([NH₄-N], [NO_x-N], and [PO₄-P]), on day 21, were analyzed using independent regression models with fauna ratios as the fixed term. The dependent variables were analyzed under the null hypothesis of no significant differences among the levels of the response variables across the different ratios of Diopatra (Dp) and Arenicola (Ar), and the Control (defaunated microcosms), and the statistical significance of their terms (p < 0.05) were tested. To deal with heteroscedasticity and avoid data transformation, the models were extended to include the appropriate variance covariate structure using a generalized least squares (GLS) estimation procedure (Zuur et al., 2009; Pinheiro et al., 2014). The variance covariate structure was determined using a restricted maximum likelihood (REML) estimation, supported by the Akaike Information Criteria (AIC) and visual comparisons of residuals plots, in contrast to the initial model. Pairwise comparisons of the treatments with different ratios of Dp and Ar were derived from each model’s t-table, re-expressed successively against each level as baseline.

Additionally, water nutrient concentrations of NH₄-N, NO_x-N and PO₄-P were represented in ordination analyses, using a Principal Component Analysis (PCA) (Anderson et al., 2008). All statistical analyses were performed using R statistical and programming environment (R Development Core Team, 2020) and the packages nlme (Pinheiro et al., 2014) for the regression analysis, and vegan (Oksanen et al., 2020) with functions metaMDS and envfit, to overlay vectors in the ordination analysis. Graphical outputs were produced with R’s package ggplot2 (Wickham, 2016).

The sediment sampled in Mira channel was classified as medium sand (125- 250 µm; φ = 2) with 8.00 ± 0.86% of fines (< 63 µm) and the total organic matter content was of 3.0 ± 0.5% (mean ± s.d.).

During the experimental period, the average water temperature in the aquaria was of 18.4 ± 1.2°C, pH 8.1 ± 0.2, oxygen concentration was 7.8 ± 0.7 mg L⁻¹ and salinity was 28.0 ± 0.2 (mean ± s.d.).

In the Control microcosms, luminophores were concentrated in the superficial layer (<2.0 cm) ( Figure 2A ), corresponding to the initial layer deposited on the top of the sediment. In the different treatments, the number of luminophores in the superficial layer tend to decrease, increasing in deeper layers with the increasing number of Arenicola ( Figures 2A–F ). All bioturbation parameters responded to the increasing density of Arenicola ( Figures 3A–D ; Table 1 ).

Surface boundary roughness (SBR) ranged between 0.26 cm (Control) and 6.48 cm (1Dp:3Ar). The highest mean SBR value was recorded in the 1Dp:3Ar treatment (4.18 ± 1.93 cm, mean ± s.d., Figure 3A ), and no significant differences in the mean SBR were found whenever Arenicola was present. Nevertheless, the monospecific treatment 4Dp:0Ar was significantly different from all the other treatments ( Figure 3A ; Table 1 ; Supplementary Table S2 ).

The mean mixed depth of particle redistribution (^f-SPIL_mean) for each aquarium ranged between 0.81 cm (4Dp:0Ar) and 5.86 cm (1Dp:3Ar), with the highest value of 3.87 ± 0.97 cm (mean ± s.d.), for the 0Dp:4Ar treatment ( Figure 3B ). There were no statistical differences between the Control and the 4Dp:0Ar ratio, between the 3Dp:1Ar, 2Dp:2Ar and 1Dp:3Ar ratios, and finally when comparing 2Dp:2Ar, 1Dp:3Ar and 0Dp:4Ar ratios ( Figure 3B ; Table 1 ; Supplementary Table S3 ).

The median mixed depth of particle redistribution (^f-SPIL_med) had its lowest value in the 4Dp:0Ar treatment, 0.83 cm, while the highest value was recorded in the 0Dp:4Ar treatment, with a ^f-SPIL_med of 6.01 cm. Collectively, the ratio 0Dp:4Ar had the highest mean ^f-SPIL_med, with 4.99 ± 1.18 cm (mean ± s.d.) ( Figure 3C ). Pairwise comparisons revealed no significant differences among Control, 4Dp:0Ar, and 3Dp:1Ar, among 3Dp:1Ar and 2Dp:2Ar, and among 2Dp:2Ar, 1Dp:3Ar and 0Dp:4Ar ( Figure 3C ; Table 1 ; Supplementary Table S4 ).

Finally, the maximum mixed depth of particle redistribution (^f-SPIL_max) ranged from 1.36 cm (4Dp:0Ar) to 11.26 cm (0Dp:4Ar). The lowest mean ^f-SPIL_max was found in the Control (1.59 ± 0.15 cm, mean ± s.d.), followed by the 4Dp:0Ar treatment (2.65 ± 1.18 cm, mean ± s.d.) while the highest mean of 8.03 ± 2.66 cm (mean ± s.d.), was recorded in the 1Dp:3Ar treatment ( Figure 3D ). When comparing the different species ratios, no significant differences were detected between the Control and the 4Dp:0Ar ratio, between the 4Dp:0Ar, 3Dp:1Ar and 2Dp:2Ar, and between 2Dp:2Ar, 1Dp:3Ar 0Dp:4Ar ratios ( Figure 3D ; Table 1 ; Supplementary Table S5 ).

The statistical model indicates that the concentration of all analyzed nutrients in the water column in the end of the experiment (21 days later) was influenced by species ratio ( Table 1 ). During the experimental period, the concentration of NH₄-N increased with the increasing number of Arenicola in the aquaria (2Dp:2Ar < 1Dp:3Ar < 0Dp:4Ar) ( Figure 4A ). However, the concentration of NH₄-N in the lower Arenicola densities (Control, 4Dp:0Ar, 3Dp:1Ar) have peaked around day 5, while in the intermediate levels (2Dp:2Ar and 1Dp:3Ar) the concentration of NH₄-N has peaked around day 16. The monospecific Arenicola treatment (0Dp:4Ar) showed an increase in the concentration of NH₄-N until day 21 ( Figure 4A ). At the end of the experiment, in the presence of fauna, the highest and lowest mean concentrations of NH4-N were registered in the monospecific scenarios ([NH4-N]_0Dp:4Ar = 8.07 ± 0.6, [NH4-N]_4Dp:0Ar = 1.49 ± 0.04 mg L⁻¹, mean ± s.d.). In the Control, the mean concentration of NH₄-N was 1.42 ± 0.03 mg L⁻¹ (mean ± s.d.). The statistical analysis revealed significant differences between all comparisons, except between 2Dp:2Ar and 1Dp:3Ar ( Figure 4A ; Table 1 ; Supplementary Table S6 ).

During the experimental period, the average concentration of NO_x-N showed an opposite trend comparatively with NH4-N, i.e., it increased in the aquaria with higher densities of Diopatra. The maximum concentration of NO_x-N was 4.44 mg L⁻¹, measured at day 16, in the monospecific treatment 4Dp:0Ar. In most treatments, NO_x-N concentrations started to increase after the 5^th day. At the end of the experiment, the lowest concentration of NO_x-N was 0.61 ± 0.15 mg L⁻¹ (mean ± s.d.), found in the 0Dp:4Ar treatment, while the highest, 4.26 ± 0.22 mg L⁻¹ (mean ± s.d.), was found in the 4Dp:0Ar treatment. For day 21, the concentration of NO_x-N in the Control was 1.64 ± 0.15 mg L⁻¹ (mean ± s.d., Figure 4B ). Contrarily to what was observed for NH₄-N, significant differences were observed in the mean concentration of NO_x-N at day 21, between the Control and the aquaria with higher ratios of Diopatra (4Dp:0Ar, 3Dp:1Ar) and also between the aquaria with 0Dp:4Ar. No differences were found the concentration of NO_x-N between the Control, 2Dp:2Ar and 1Dp:3Ar treatments ( Figure 4B ; Table 1 ; Supplementary Table S7 ).

For PO₄-P, the Control showed a reduction from day 0 to day 21, while the monospecific 4Dp:0Ar treatment had only a marginal decrease. The remaining treatments showed an increase in the concentration of PO₄-P from day 0 to day 21 ( Figure 4C ). The lowest concentration of PO4-P was 0.02 mg L⁻¹ (day 12, Control), while the highest was 0.24 mg L⁻¹ (day 21, 1Dp:3Ar). The mean concentration of PO₄-P tended to increase with the increasing ratio of Arenicola, achieving a maximum concentration of 0.20 ± 0.02 mg L⁻¹ (mean ± s.d.) for the monospecific Arenicola treatment (0Dp:4Ar) in day 12 ( Figure 4C ). The statistical analysis for the concentration of PO₄-P at day 21 showed no differences between 3Dp:1Ar, 2Dp:2Ar, and 1Dp:3Ar ratios (respectively 0.14 ± 0.03, 0.13 ± 0.10. and 0.19 ± 0.04 mg L^-1, mean ± s.d.), and between 1Dp:3Ar and 0Dp:4Ar (0.19 ± 0.01 mg L^-1, mean ± s.d.) ratios, ( Figure 4C ; Table 1 ; Supplementary Table S8 ).

The principal components ordination clearly shows the evolution of the mesocosms water concentration along time, clustering the samples from days 0, 2 and 5 along Axis 1, which explains almost 60% of the total variance ( Figure 5 ). The nutrients concentration vectors are spread also along Axis 2 and showed that the Control, the monospecific scenario 4Dp:0Ar and the ratio 3Dp:1Ar were characterized by higher concentrations of NO_x-N, while the other treatments (2Dp:2Ar, 1Dp:3Ar, 0D:4Ar) were defined by higher concentrations of NH₄-N and PO₄-P, particularly for the dates after day 8 ( Figure 5 ).