Taxonomic and functional diversity of benthic macrofauna associated with rhodolith beds in SE Brazil

Study area and sampling design

The study area is located within the limits of the Costa das Algas Marine Protected Area (MPA) on the Eastern Marine Ecoregion of Brazil (Fig. 1; Table 1; Spalding et al., 2007). This is a tropical region characterized by rainy summers, with predominantly NE and E winds, and dry winters (Bernardino et al., 2015). The continental shelf on Eastern Brazil is influenced by the Tropical Water (TW) of the Brazil Current, with temperatures above 22 °C and salinities above 36 (Mazzuco et al., 2019), and eventual seasonal summer upwelling of the South Atlantic Central Water (SACW) into the shelf, with temperatures between 6 °C and below 20 °C and salinities between 34.6 and 36 (Quintana et al., 2015; Mazzuco et al., 2019). The continental shelf on the Espírito Santo basin includes a mixed system of terrestrial and carbonate sediments with rhodolith beds extending down to the shelf break at depths over 80 m (Figueiredo et al., 2015; Holz et al., 2020).

Sampling was carried out by SCUBA diving in January 2019 and the density of the rhodolith beds was a determining factor of the sampling design. Based on preliminary images of the area, sampling stations were classified into two categories (beds): high-density (H1, H2, and H3) and low-density (L1, L2, and L3; Fig. 1). High density beds had 100% of the seafloor covered by rhodoliths, and low-density beds had scattered patches of sediment between rhodolith nodules. The differences in bed structure were further confirmed by measurements of rhodolith morphology, internal volume, and ramification (branch density, see below). Abiotic metadata (temperature, salinity, depth, and water visibility) were obtained at the time of sampling using a CTD and Secchi disk. In each station, scuba divers sampled manually all rhodoliths on the surface (50–100 rhodoliths in high-density beds; 4–40 rhodoliths in low-density beds) within a 0.25 m² quadrat from three replicates randomly distributed along a 20 m random transect (Fig. 2). Occasional megafaunal organisms were observed (macroalgae, ophiuroids) but were not sampled and were thus excluded from our analysis. Triplicate samples of the underlying unconsolidated sediment within each quadrat were collected using PVC corers of 10 cm in diameter (10 cm depth) with sealing lids. The rhodolith nodules were packed in cloth bags (<0.5 mm mesh) and sealed to prevent loss of macrofauna during recovery on board, where they were immediately fixed with formaldehyde (10%) and borax to avoid carbonate degradation. Macrofaunal invertebrates in the sediments were also fixed with formaldehyde (10%) and borax. In addition to the biological data, sediment underlying the rhodolith beds was also sampled for grain size, organic matter, and biopolymers analysis using triplicate corers (10 cm diameter), preserved in ice on board, and frozen until processing. Sampling of benthic organisms was authorized by the permit N 24700-1 (MMA-ICMBio) and voucher specimens are deposited at the Unicamp Zoology Museum (http://www.splink.org.br).

Laboratory analysis

The classification of rhodolith morphology was determined by measuring the largest, intermediate, and minor axis of each nodule, which resulted in a mean nodule diameter and sphericity for each station (Bosence & Pedley, 1982). The morphological rhodolith dataset was plotted on a TRIPLOT spreadsheet developed by Graham & Midgley (2000), and drawn on the pebble shape diagram of Sneed & Folk (1958) that discriminates rhodoliths in spheroidal, discoidal, or ellipsoidal shape. The rhodolith bed density was estimated from the number of nodules sampled within each quadrat (nodules·m⁻²). The mean rhodoliths diameter within each station was measured from all nodules in each replicate and averaged per station from the three replicated samples.

The ramification of rhodolith nodules were determined semi-qualitatively from each site based on their relative branching density (Bosence, 1983). Nodules were classified into four groups: I = a single branch; II = few branches; III = frequent branching; IV = dense and solid branching. The average internal volume of the rhodoliths in each site was determined by water volumetric displacement (Steller et al., 2003). Rhodoliths were covered with a plastic film and then submerged in a graduated container of a known volume. Later, rhodoliths were again submerged but without the plastic film. The difference between the initial volume (rhodoliths with film) and the displacement of the liquid with rhodoliths without the film was considered the rhodolith’s internal volume.

In the laboratory, rhodolith nodules were broken and epifaunal organisms were sieved (500 μm) and preserved in 70% ethanol until sorting. All organisms were identified to family or the lowest possible taxonomic level under a stereomicroscope. Macrofaunal trophic group analysis followed the main feeding types including deposit feeders, detritivores, suspension feeders, filter feeders, and omnivores, carnivores or others (OCO) according to Arruda, Domaneschi & Amaral, 2003, Jumars, Dorgan & Lindsay (2015), Queirós et al. (2013) and Macdonald et al. (2010).

For sediment grain size, carbonate content, and organic matter analysis, the samples were thawed and placed in an oven at 60 °C for 48 h. The dry sediment was macerated and taken to a stirrer, where the grain size was determined by sieving it between −1.5 phi (Φ) sieves and 4 Φ, with 1 Φ intervals. Subsequently, the values of Φ were added to the SysGran 3.0 software (Camargo, 2006) to analyze the granulometric properties (i.e. average grain size and the total percentage of gravel, sand, silt, and carbonate). The carbonate contents of the sediment were determined by combustion in a muffle (550 °C for 4 h) with an additional hour at 800 °C. The sediment organic matter content was determined by combustion in a muffle (550 °C for 4 h). Due to sample loss, sedimentary analysis was not done in station L3. All sedimentary organic biopolymers (carbohydrates, lipids, and proteins) analysis were made in triplicates, following the methods in Danovaro (2010). Total protein analysis (PRT) was carried after its extraction with NaOH (0.5 M, 4 h) and was determined according to Hartree (1972), modified by Rice (1982), to compensate for phenol interference. Total carbohydrate (CHO) was analyzed according to Gerchacov & Hatcher (1972). Total lipids (LPD) were analyzed according to the protocol described in Marsh & Weinstein (1966), being extracted from 1 g of homogenized sediment lyophilized by 11 ultrasonication (20 min) in 10 ml of chloroform:methanol (2:0 1 v/v). Blanks for each analysis were taken with pre-combusted sediments at 450 °C and 480 °C for 4 h. The concentrations of total protein, carbohydrate, and lipid were displayed as bovine serum albumin (BSA), glucose, and tripalmitin equivalents, respectively. The concentrations of total protein, carbohydrate, and lipid were converted into carbon equivalents assuming a conversion factor of 0.49, 0.40, and 0.75, respectively (Fabiano & Danovaro, 1994). Also, protein to carbohydrate (PRT:CHO) and carbohydrate to lipid (CHO:LPD) ratios were applied to assess the state of biochemical degradation processes (Galois et al., 2000). The sum of biopolymer concentrations was added to the analysis as a measure of biopolymeric carbon (BPC; Fabiano, Danovaro & Fraschetti, 1995; Hadlich et al., 2018).

Statistical analysis

Benthic assemblages were compared across rhodolith beds with high and low densities (bed), and between rhodoliths and sediments (habitats). Within this sampling design we compared the nodule’s epifauna with the infaunal macrofauna of the underlying soft sediments. Macrofaunal density and diversity (alpha, gamma, and beta) were compared across bed (high and low density) and habitat (nodule and sediments) levels. We used the approach of additive partitioning to estimate alpha (α) and gamma (γ) species richness based on the sum of richness from each sample (α) within each habitat (Josefson, 2009). In addition, the macrofaunal trophic diversity was determined based on the classification of species according to their feeding modes (see Laboratory analysis description) and their density. This matrix was used to calculate the assemblage functional richness (FRic), functional dispersion (FDis), functional evenness (FEve), and entropy (FRaoQ; Mason et al., 2005); which were also tested over the spatial scales above. FRic indicates the amount of niche space (feeding mode category) filled by species in the community, FEve describes the evenness of functional distribution in a the trait space, and FDis and RaoQ are indices quantifying how functionally similar is the community spatially (Botta-Dukát, 2005; Mason et al., 2005).

Spatial differences in rhodolith bed structure (nodule density, internal volume, morphology and ramification) as well as in sediment parameters (total organic matter, carbonate, biopolymers, and granulometry) were also compared across bed and habitat levels. Spatial analyses for bed structure, sediments, and macrofaunal assemblages were made either through analysis of variance (ANOVA; Underwood, 1997) for univariate parameters or by a permutational multivariate analysis of variance (PERMANOVA; Anderson, 2017) for multivariate parameters. These tests were hierarchically designed with one fixed factor (beds, two levels: high and low); or two fixed factors, adding station (nested in bed, three levels: stations 1, 2, and 3) or habitat (fixed, two levels: rhodoliths and sediment). PERMANOVAs were based on a Bray–Curtis resemblance matrix under a reduced residuals model and data was square-root transformed to give more weight to rare taxa in the analyses (Clarke & Gorley, 2006). Post-hoc pairwise tests (Tukey for ANOVAs or PERMANOVA pairwise) were performed to identify significant differences within factor levels (Underwood, 1997; Anderson, 2008).

Macrofaunal assemblage composition was also tested based on the total dissimilarity (i.e. beta diversity) between beds and habitats. Dissimilarity levels were calculated from both species’ turnover (βSIM; total replacement of species) and nestedness (subsets of species among sites, βSN; Baselga, 2010), based on macrofaunal presence-absence data (Sørensen index). A cluster dendrogram was also applied using the average abundance of all taxa from the Bray-Curtis similarity matrix. Additionally to that, a non-metric Multidimensional Scaling Analysis (nMDS) was applied to visualize the similarities of macrofauna assemblies between densities and habitats, using the square-root abundance of all taxa from the Bray-Curtis similarity matrix. A canonical analysis of principal coordinates (CAP; Anderson & Willis, 2003) was performed to determine the association between environmental variables and benthic assemblages between beds and habitats. Graphic design and analysis were performed using R Project (R Core Team, 2014) with packages: ‘ggplot2’ (Wickham, 2016), ‘oce’ (Kelley & Richards, 2017), ‘stats’, ‘vegan’ (Oksanen et al., 2018), MASS (Ripley et al., 2019), mgcv (Wood, 2012), MuMIn (Barton & Barton, 2013), FD (Laliberté, Legendre & Shipley, 2014), and ‘ggdendro’ (Andrie & Ripley, 2020). Macrofaunal raw data is openly available through the Ocean Biogeographic Information System (OBIS) portal (Stelzer et al., 2021).

Rhodolith bed characteristics

The maximum depth of sampling stations varied between 39 and 55 m and did not differ significantly between rhodolith beds of high and low density (Table 1). During sampling, surface and bottom temperature ranged from 26–28 °C and 19–23 °C, respectively, while salinity ranged from 37.7–38.3 (Table 1; Fig. S1). Stations H1 and H2 had water column profiles with a marked halocline in the first 10 m, whereas at stations H3, L2, and L3 the halocline occurred at 15 to 35m depth. Temperature showed a similar bathymetric profile between stations. Secchi’s depth varied between 20–35 m deep and the incidence of light in the water column reached greater depths in station L2 (Table 1).

The sediments under the high and low-density beds had a similar grain size composition with predominance of coarse and medium sand (F = 3.51, p = 0.07; Table 2). The percentage of sediment total organic matter varied between 2% and 5%, with no significant differences between high and low-density beds (F = 0.038, p = 0.849). Sediment carbonate content ranged from 3% and 8%, with higher carbonate content under the high-density beds (F = 5.74, p = 0.0323). Two sediment organic biopolymers (proteins and lipids) had higher concentrations in high-density beds (F = 27.3, p = 0.0002; Table 2, Fig. 3), whereas carbohydrate concentrations were similar. Protein (PRT) and lipid concentrations in sediments of high-density beds ranged from 0.43 to 0.85 mg·g⁻¹ and 0.02 to 0.30 mg·g⁻¹, respectively (Fig. 3). The total biopolymeric carbon and biopolymer ratios (protein: carbohydrate, carbohydrate:lipid) did not differ at the scale of beds, and showed a wide range among samples (0.23 to 2.81 protein:carbohydrate, 2.4 to 111.3 carbohydrate:lipid). Biopolymeric carbon concentration did not vary in sediments under rhodolith beds.

Rhodolith nodule density in high-density beds (204 ± 18.7 m⁻²) was over 3 times higher when compared to the low-density beds (61 ± 27.1 m⁻²; F = 47.9, p < 0.0001; Table 3). High-density beds were dominated by rounded (36%) and elongated (15%) nodules, whereas in low-density beds the nodules were predominantly discoidal (22%) and spherical (14%; F = 3.12, p = 0.05; Table 3, Fig. 4). Overall, nodule mean diameter and the ramification pattern were similar between high-density and low-density beds, with small nodules (3.7 to 5.3 cm mean diameter) with dense and solid branching (88% of ramification type IV; p < 0.01; Tables 1 and 3). An exception of this patter was observed at H1 where nodule diameter was significantly lower than at the other stations (F = 7.69, p = 0.0071, F = 5.17, p = 0.05; Tables 1 and 3). The diameter of rhodolith nodules at low-density beds were more heterogeneous than at the high-density beds (F = 7.69, p = 0.0071; Table 3). The average rhodolith nodule internal volume in high-density beds (312.2 cm³·m⁻²) was similar with the low-density beds (206.3 cm³·m⁻²), but with local differences between stations H3 and L3 (F = 5.68, p = 0.0183; Table 3).

Macrofaunal assemblages

We sampled a total of 11,421 macrofaunal organisms associated with the rhodolith beds (epifauna) and underlying sediments (infauna), and registered significant differences between high and low-density beds and habitats (rhodoliths and sediment). Rhodolith nodules in high-density beds had a similar macrofaunal density (2,736 ± 161.7 ind·m⁻²) when compared to low-density beds (2,339 ± 554.9 ind·m⁻²; F = 0.27, p = 0.601; Table S1, Fig. 5); but in both cases, macrofaunal density was significantly higher on nodules (2,538 ± 288.7 ind·m⁻²) when compared to the sediments below nodules (166 ± 38.8 ind·m⁻²; F = 66.18, p < 0.0001). We identified 151 different macrofaunal taxa within nodules and sediments, with pronounced contrasts between both habitats. There was a higher alpha and gamma diversity in rhodoliths when compared to the sediment infauna (F = 38.6, p = 0.0002).

Macrofaunal composition changed markedly at the scale of beds (F= 3.41, p = 0.02) and between habitats (F = 20.0, p = 0.01; Table S1). Rhodolith nodules were mostly dominated by Annelida (47%) and Crustacea (44%) with a marked difference between high and low-density beds. Annelida (Syllidae, Nereididae, and Lysidice sp) dominated high-density beds (51%), whereas Crustacea (Gammaridae, Melitidae, and Elasmopus sp) dominated (64%) low-density beds (Fig. S2). In contrast, sediment macrofaunal assemblages were relatively similar under high and low-density beds being dominated by Crustacea (63%; Ostracoda, Melitidae, and Leptochelia sp) and Mollusca (22%; Meioceras sp and Cardiidae; Fig. S2).

Macrofaunal trophic richness exhibited a distinct dominance between habitats (Fig. S2). The dominant trophic groups on the rhodolith nodules were omnivores, carnivores, and others feeders (OCO, 74.6%), while sediments were dominated by filter feeders (38.5%) and OCO (28.5%). Macrofaunal trophic richness (FRic) was higher on the rhodolith nodules (7 ± 0.2), when compared to the sediment underneath (5 ± 0.4; F = 27, p < 0.0001; Fig. 6). Sediments presented a more homogeneous trophic evenness (FEve = 0.32 ± 0.03) when compared to the rhodolith nodules (FEve = 0.09 ± 0.02; F = 75.39, p < 0.0001; Table S2). Macrofaunal trophic dispersion (FDis) and entropy (FRaoQ) differed between habitats only at low-density beds, which were both higher in sediments (Table S2; Fig. 6).

There was a high degree of patchiness with high dissimilarity among macrofaunal assemblages between rhodolith nodules and sediments (Sørensen Index = 0.92; Fig. S3). Macrofaunal composition in rhodoliths and underlying sediments was marked by low nestedness and high taxa turnover (Table 4). The CAP analysis supported the spatial dissimilarity in macrofaunal assemblages inside nodules, suggesting a strong effect of bed density on nodule’s macrofaunal density (F = 7.97, p = 0.002; Fig. 7, Table S3). On the other hand, the sediment macroinfauna under rhodolith beds were significantly related to sediment carbonate content (F = 1.71, p < 0.01; Table S4, Fig. 7). The CAP ordination showed that bed density influenced the rhodolith macrofaunal composition with greater contribution of Syllidae and Sipuncula in high-density beds, whereas low-density beds were associated with gammarid amphipods (Fig. 7). The underlying sediments at high-density beds were predominantly dominated by ostracods, but crustaceans Melitidae and Amphiuridae dominated at low-density beds (axis 1; >29%).