Benthic macrofaunal structure and secondary production in tropical estuaries on the Eastern Marine Ecoregion of Brazil

Study area

The study was carried out in three tropical estuaries with a microtidal, semidiurnal tidal pattern within the Eastern Brazil Marine Ecoregion (Spalding et al., 2007; Fig. 1). The northernmost estuary, Piraquê-Açu-Mirim estuary (PAE; 19°57′S 40°09′W) is within a municipal conservation unit and covered by extensive and well-developed mangrove forests with an area of over 19 km² (Bernardino et al., 2018; Servino, Gomes & Bernardino, 2018). The Vitória Bay estuarine system (VIB; 20°16′S 40°20′W) is the largest estuary in the region with approximately 18 km² of mangrove forests and surrounded by a densely populated metropolitan area with high levels of sewage input and industrial activities (Jesus et al., 2004). The southernmost estuary, Benevente estuary (BEN, 20°48′S 40°39′W), has well preserved mangrove forests that cover an area of approximately 4.6 km² with minor urban settlement (Pereira et al., 2009; Petri et al., 2011). Mangrove forests of the three estuaries are composed by Rhizophora mangle, Laguncularia racemosa and Avicennia schaueriana species.

Sampling and sample processing

Benthic macrofaunal assemblages were sampled on a nested spatial design on vegetated (V—mangrove forests) and unvegetated (NV—tidal flats) habitats on the mesohaline sectors (salinity range between 18 and 5; Venice system, 1958) of the three estuaries (VIB, PAE and BEN, ICMBIO permit N 24700-1). Salinity sectors in the sampled estuaries were measured with either a multiparameter or with semi-continuous (5–20 days) deployment of data-loggers. Sampling occurred in one sampling event in each estuary between August and September 2014, during low tides and on the dry season. Each estuary was divided in two sites distanced in the scale of hundreds of meters containing adjacent vegetated and unvegetated habitats (Fig. A1). Three sampling plots were randomly established in each habitat and site, parallel to the waterline and separated by tens of meters. Three replicate faunal samples were sampled within each plot, distanced by approximately 1 m from each other using a PVC corer with 15 cm diameter (0.0177 m² area) and to a sediment depth of 10 cm. Additionally, one composite sediment sample was collected at each plot for sediment analysis (grain size, total organic matter and chlorophyll-a), by mixing three samples of 7 cm diameter and 5 cm depth. Superficial water temperature and salinity were measured in each sampling area.

Faunal samples were preserved in 4% formalin and posteriorly washed through a 1 mm sieve and the retained material was stored in 70% ethanol. In the laboratory, samples were sieved through a stacked series of sieves (1, 1.4, 2, 2.8 and 4 mm), using the methods described by Edgar (1990a). Macrofauna was sorted in each sieve size and identified at family level, considering that this level of identification is satisfactory to identify differences in macrofaunal assemblages for the aims of this study (Warwick, 1988; Chapman, 1998; Olsgard, Somerfield & Carr, 1998; De Biasi, Bianchi & Morri, 2003). During sorting of samples, the plant material was separated for plant biomass (plant detritus and living roots) determination (dry weight) after drying at 60 °C during 48 h.

Sediment subsamples were treated with hydrogen peroxide (H₂O₂), to eliminate organic matter, and mud content was determined by wet sieving samples through a 0.063 mm mesh size. After drying, the sediment >0.063 mm was sieved through a series of sieves and grain size was classified following the Wentworth scale (Suguio, 1973). Sediment total organic matter (TOM) content was estimated by weight loss after combustion at 500 °C during 4 h. Chlorophyll-a (Chl-a) and phaeopigments were extracted from the sediment with acetone and analyzed using a spectrophotometer before and after acidification with HCl (Lorenzen, 1967; Quintana et al., 2015).

Faunal biomass and secondary production

Macrofauna was wet weighed within each taxonomic group, generally family, by each sieve size (1, 1.4, 2, 2.8 and 4 mm) after identification. Macrofaunal biomass (mg wet weight m⁻²) was converted into ash-free dry weight (mg AFDW m⁻²) using the conversion factors compiled in Brey (2001) and Brey et al. (2010). Shells of mollusks were excluded from biomass estimation. Conversion factors from Brey (instead of estimate by methodology used by Edgar (1990a)) were chosen to avoid overestimation of AFDW and consequently of production, mainly in the larger sieve size, since some individuals with elongated bodies and low weights were retained in the sieves.

The secondary production of benthic macrofauna was estimated using the general equation (P) = 0.0049∗B^0.80∗T^0.89 of Edgar (1990a), which relates daily macrobenthic production P (µg day⁻¹) to ash-free dry weight B (µg) and water temperature T (^∘C). The water temperature measured in all three estuaries during faunal sampling had a small range so the temperature (T) used in the estimation of benthic secondary production was standardized at 23.5°C. The use of a standardized temperature does not show seasonal, daily or spatial variations, thus limiting secondary production estimates (Edgar, 1990a; Tumbiolo & Downing, 1994). Although with a limited temporal applicability, we used a standardized temperature to indicate the relative secondary production between different estuaries (Edgar, 1993; Edgar & Barrett, 2002). Production was calculated for each taxon (Polychaeta, Oligochaeta, Kalliapseudidae, Other Crustacea, Mollusca and Others) within each sieve size and total production per sample was calculated as the sum of these values. The annual production to biomass ratio (P/B) for each habitat in each estuary was calculated from mean production divided by the mean macrofaunal biomass. This parameter can be considered a measurement of biomass turnover rate (Dolbeth et al., 2012).

Data analysis

The spatial variability in benthic macrofaunal density, biomass and secondary production were evaluated at multiple spatial scales in different habitats using a nested and orthogonal analysis of variance (ANOVA). Habitat was defined as a fixed factor and orthogonal to spatial factors (estuary, site, plot). Spatial factors were treated as random and included three estuaries, sites (N = 2) nested in estuary, plots (N = 3) nested in site and samples (N = 3) collected at plots. Spatial differences on sediment properties and plant biomass were assessed by ANOVA across scales of estuary, and site (nested in estuary), due to the lack of sample replication at plots. This ANOVA also included habitat factors orthogonal to spatial factors since both vegetated and unvegetated habitats were sampled. A Cochran’s test was performed previously to each ANOVA to assess homogeneity of variances and when necessary data was transformed.

Differences among macrofaunal assemblages were assessed by a Permutational Multivariate Analysis of Variance (PERMANOVA) that was performed using Bray–Curtis similarity coefficients on square-root transformed abundance data (9999 permutations, Anderson, Gorley & Clarke, 2008). A non-metric multidimensional scaling (nMDS) performed using Bray–Curtis similarity matrix and square-root transformed data was used to visualize variation in macrofaunal assemblages. A Similarity Percentage analysis (SIMPER) was used to identify the taxa that most contributed to dissimilarities among habitats. The relationships between sediment properties (TOM, Chl-a, Mud, plant biomass) and macrofaunal density were investigated using a Canonical Correspondence Analysis (CCA). In this analysis, the density of the top five dominant taxa (comprising over 90% of total density) was used instead of the complete data. In the CCA, the sum of macrofaunal density data from all replicates was used so that the number of samples from density and sediment properties was the same. Part of statistical analyses were performed in the software R (R Core Team, 2015). The GAD package was used to perform ANOVA analysis (Underwood, 1997; Sandrini-Neto & Camargo, 2014) and the Vegan package was used to perform nMDS and CCA analysis (Oksanen et al., 2017). PERMANOVA was carried out using the software PRIMER 6 with the PERMANOVA + add-on package (Clarke & Gorley, 2006; Anderson, Gorley & Clarke, 2008).

Sediment properties and plant material

The sediment was predominantly composed of mud in all estuaries and habitats (Fig.  2). When comparing the three estuaries, the sediment mud content, mean grain size and total organic matter differed significantly among sites and in the interaction between habitat and site (Table 1). These results represent spatial variation at local scales. Chl-a and phaeopigments differed significantly between habitats and estuary, respectively, with higher sediment Chl-a at unvegetated habitats and lower phaeopigments in the BEN estuary (Table 1, Fig. 2). Plant biomass differed significantly among estuaries and in the interaction between habitat and site (Table 1). VIB presented over two times higher total plant biomass when compared to similar sectors in the BEN and PAE estuaries (Fig. 2F).

Macrofaunal density, biomass and secondary production

A total of 18,036 individuals belonging to 37 taxa were sampled at the three estuaries, of which 7,989 individuals were Kalliapseudidae (Tanaidacea). BEN estuary had a total of 11,481 individuals, distributed in 23 taxa. In PAE estuary 1,728 individuals were collected and distributed in 27 taxa. VIB had a total of 4,827 individuals, distributed in 28 taxa. Within the mesohaline sector of the three estuaries, total macrofaunal density was significantly different at the plot and estuary spatial scales, and in their interactions with habitat (Table  2). The BEN estuary presented higher macrofaunal density in unvegetated habitats, but this pattern was opposite to the VIB and PAE estuaries that had no density differences between habitats (Fig. 3A). Macrofaunal densities varied over 30-fold between unvegetated habitats at BEN and PAE estuaries (33,022 ± 14,709 ind m⁻² and 1,033 ± 1,558 ind m⁻², respectively; Fig. 3A). Kalliapseudidae (Tanaidacea) was dominant in unvegetated tidal flats at BEN estuary, whereas Polychaeta and Oligochaeta were more abundant in similar habitats at PAE and VIB estuaries. Vegetated habitats in the three estuaries had higher densities of Oligochaeta and Polychaeta (Fig. 4A).

Significant differences in macrofaunal biomass and estimated secondary production were observed only in the interaction between habitat and estuary when comparing the three estuaries (Table 2). Biomass and production followed patterns of macrofaunal density and were not clearly distinct between unvegetated or vegetated habitats among the three estuaries studied (Figs. 3B and 3C). The lowest macrofaunal biomass and production were observed at unvegetated tidal flats in the PAE and VIB estuaries (Figs. 3B and 3C).

The contribution from each macrofaunal group to total assemblage biomass and secondary production varied greatly between estuaries and habitats (Fig. 4). Large individuals including bivalve molluscs and brachyuran crabs contributed greatly to benthic biomass and production in vegetated habitats at the three estuaries despite their low density (Figs. 4 and 5). At vegetated habitats in VIB estuary, Mollusca (mainly Mytilidae and Solecurtidae) contributed to most of the biomass (1,832.5 ± 2,780.5 mg AFDW m⁻²) and production (28.3 ± 37.8 mg m⁻² day⁻¹), with Oligochaeta and Polychaeta representing second and third groups respectively. At vegetated habitats of the PAE estuary, Mollusca (mainly Mytilidae; 2,864.6 ± 8,115.6 mg AFDW m⁻², 35.1 ± 82.9 mg m⁻² day⁻¹) and Crustacea (mainly Panopeidae; 1,199.4 ± 4,331.9 mg AFDW m⁻², 15.6 ± 49.3 mg m⁻² day⁻¹) were the most representative groups in biomass and production. Crustaceans (mainly Ocypodidae; 1,897.8 ± 3,682.9 mg AFDW m⁻², 28.3 ± 46.5 mg m⁻² day⁻¹) contributed to over 70% of the macrofaunal biomass and production in vegetated habitats at the BEN estuary with Polychaeta as the second group.

In general, macrofaunal biomass and production of estuarine habitats were mainly derived from large size classes (Figs. 5B and 5C). Vegetated habitats had over 70% of its production from large size classes (>4 mm), whereas unvegetated habitats had more variable contribution (43–87%) of the other size classes from 1 to <4 mm (Figs. 5B and 5C). At unvegetated habitats in the mesohaline sector of VIB (329.4 ± 759.3 mg AFDW m⁻², 6.6 ± 12.6 mg m⁻² day⁻¹ of Mollusca) and PAE (51.3 ± 193.2 mg AFDW m⁻², 1.3 ± 4.2 mg m⁻² day⁻¹ of Mollusca) estuaries, Mollusca (mainly Solecurtidae) and Polychaeta (mainly Capitellidae) contributed significantly to total macrofaunal biomass and production (Figs. 4B and 4C). Kalliapseudidae was the dominant taxa at unvegetated habitats in BEN estuary (7,315.7 ± 5,343.6 mg AFDW m⁻², 126.8 ± 86.8 mg m⁻² day⁻¹) and contributed greatly to biomass and production (over 90%; Figs. 4B and 4C).

The mean estimated community annual production to biomass ratio (P/B) varied among estuaries and habitats. The highest P/B ratio was observed at unvegetated habitats at PAE estuary (12.6 y⁻¹), whereas vegetated habitats in this estuary had the lowest P/B ratio (5.3 y⁻¹). P/B ratios did not vary significantly between habitats at BEN (6.4 and 6.6 y⁻¹ for V and NV respectively) and VIB estuaries (7.5 y⁻¹ and 9.3 y⁻¹ for V and NV habitats respectively).

Assemblage composition

Macrofaunal assemblages differed markedly between vegetated and unvegetated habitats and between estuaries (Table 3). The numerically dominant taxa in vegetated habitats in the three estuaries were Oligochaeta and Capitellidae (>90%). In the unvegetated habitats the numerically dominant taxa were more variable among the estuaries. At BEN estuary Kalliapseudidae and Oligochaeta (>98%) were dominant. In unvegetated habitats at VIB Spionidae and Capitellidae (>80%) were more abundant, whereas at PAE estuary Capitellidae and Oligochaeta (75%) dominated in unvegetated habitats. Although differences among the dominant taxa between unvegetated habitats at BEN, VIB and PAE, all three estuaries had most taxa shared between them.

The macrofaunal assemblage composition was significantly different in several spatial scales within the mesohaline sector of the three estuaries (PERMANOVA; Table 4). These significant differences occurred in the interaction among habitat and all the spatial scales analyzed (estuary, site and plot) and the spatial scales within estuaries (site and plot). Faunal distribution patterns in nMDS ordination evidenced differences between unvegetated and vegetated habitats, with more heterogeneous assemblages in the former if compared to tightly grouped vegetated samples (Fig. 6). This pattern of higher spatial variability was also observed among the three estuaries, where macrofaunal assemblages at unvegetated habitats had lower similarity if compared to vegetated habitats.

Dissimilarities were high (>60%) between habitats inside each estuary and among estuaries in the unvegetated habitat (SIMPER). Kalliapseudidae, Oligochaeta, Capitellidae and Ampharetidae were the taxa that most contributed to the observed differences among habitats in the mesohaline sector of BEN (SIMPER; Table A1). At VIB and PAE, Oligochaeta, Spionidae, Capitellidae, Nereididae and Pilargidae were the taxa that most contributed to the observed differences among habitats in the mesohaline sector (SIMPER; Table A1). The dissimilarity between unvegetated habitats among estuaries within the mesohaline sectors occurred mainly due to differences in abundance of Kalliapseudidae (BEN), Spionidae (VIB) and Oligochaeta (PAE; SIMPER; Table A2).

Relationships between sediment properties and macrofaunal assemblages

Macrofaunal assemblages were related to sediment mud content, TOM, plant biomass and Chl-a, with the first and second canonical axes explaining 26% and 17.2% of the variation in the data, respectively (CCA; Fig. 7). These relationships also explained the differences in assemblage composition between vegetated and unvegetated habitats. The CCA evidenced differences between habitats and estuaries and three groups of samples were formed in the CCA. The first group was corresponding to unvegetated habitat in VIB, the second group to unvegetated habitats in BEN, and the third group was formed by samples from both habitats in PAE, vegetated habitat in VIB and in BEN. Vegetated habitats of the three estuaries were related to higher TOM content, higher plant biomass and to higher densities of Oligochaeta and Capitellidae. Nereididae was also a family with high densities at vegetated habitats in PAE. Unvegetated habitats were more heterogeneous between estuaries, with VIB exhibiting higher Chl-a and dominated by Spionidae, whereas at PAE Capitellidae was dominant. At BEN, Kalliapseudidae was abundant at unvegetated sediments with high mud content and relative low plant biomass and TOM content.