The IRL is a barrier island lagoon that extends ∼250 km along the central east coast of Florida with a width of 1–8 km, mostly wind-driven currents and limited oceanic exchange via just five inlets to the Atlantic Ocean. Salinities in the IRL typically range from 15 to 40 (Morris and Virnstein, 2004). The average depth of the IRL is ∼1.7 m (Smith, 2007).

Samples were collected from nine areas along ∼60 km of the IRL from Cocoa to Saint Sebastian River (Figure 1). These samples included mouths of tributaries plus areas known from previous studies to have anoxic, organic-rich, muddy sediment. The average water depth for areas with muddy sediments was 2.5 ± 1.0 m with a range of 0.9–4.9 m. These muddy sediments cover an area of <10% of the lagoon; however, they were estimated to be the source of >25% of total dissolved inorganic N and P fluxes to the lagoon (Trefry and Trocine, 2011; Fox and Trefry, 2018; Tetra Tech Inc and Closewaters LLC, 2021).

Two stations were sampled in each of eight areas of the lagoon, plus one station in South Valkaria, for a total of 17 stations (Figure 1). Station locations identified with a P (e.g., SS-P) were on the perimeter of an organic-rich, muddy deposit where mixing with sand was more likely. Station locations identified with a C (e.g., SS-C) were in the center of muddy deposits that were especially rich in OM and H₂S. Therefore, two nearby stations with significantly different sediment textures and compositions were sampled in each area. These two stations had a high likelihood of being colonized by the same species; however, sediment composition for the two stations was typically very different.

This study was conducted from November 2015 to July 2016. All 17 stations were sampled twice, once when water temperatures averaged 19 ± 3°C (dry season) and once at 29 ± 2°C (wet season). Each pair of stations (e.g., CO-P and CO-C) was sampled on the same day during each sampling period. Continuous water column profiles for temperature, DO, % DO saturation (DO_sat), salinity and pH were obtained using a Yellow Springs Instrument 6600 V2.0 Sonde. The sonde was calibrated prior to each day-long sampling trip following manufacturer’s specifications.

Water depths and thicknesses of muddy sediments overlaying sand layers were determined at each station using a marked and capped 4.2-cm diameter polyvinyl chloride pole. The pole was inserted vertically into the water until it touched the surface of the sediment; that depth was recorded to the nearest 0.1 m as the water column depth. Then, the pole was pushed into the muddy sediment until it reached hard bottom (usually sand), with the resulting total depth representing the depth of the water column plus the thickness of muddy sediment. Thicknesses of muddy sediments were determined by subtracting water column depth from total depth.

Triplicate sediment samples were collected at each station using an Ekman grab (15 cm × 15 cm × 19 cm). Sediment from each grab was sieved through a 500-μm, stainless steel mesh and organisms were placed in separate 700-mL containers with water from the sample site. A fourth grab was collected at each station for chemical analysis of the sediment. A clear acrylic tube (7 cm diameter, 24.5 cm long), with holes every 2 cm along its length for subsequent electrode insertion (plugged during deployment), was placed vertically into the center of the fourth grab, then removed and sealed for determination of DO, Eh, and pH. Four 10-mL syringes with the forward ends cut off (1.3 cm diameter, 8 cm long) were used to collect mini-cores from the chemistry grab (Fox and Trefry, 2018). These syringe samples were sealed immediately with parafilm, kept on ice and returned to the laboratory for determination of porewater concentrations of dissolved NH₄⁺ and H₂S. Sediment also was collected from the upper 5 cm of the chemistry grab and placed in double Ziploc bags. This sediment was used to determine water content, grain size, loss on ignition at 550°C (LOI), CaCO₃, TN, TOC, and porosity. All samples were kept cold (∼4°C) until analyzed or freeze-dried.

Four mini-cores from each site were placed into N₂-purged tubes and centrifuged for 7 min to obtain porewater. The supernatant (porewater) was then filtered through Whatman, 0.45-μm polypropylene, membrane filters. The filtered water was immediately analyzed for dissolved H₂S by Hach Method #8131 with a precision of 10% relative standard deviation (RSD). Filtered water was stored at 4°C and purged with N₂ gas before being analyzed for NH₄⁺ to limit analytical interference from H₂S. Concentrations of NH₄⁺ were determined colorometrically using the standard phenate method #4500-NH₃ (Clesceri et al., 1989) with an analytical precision of 3% (RSD) and 3% accuracy.

Concentrations of DO in sediments were determined using an OM-4 oxygen microelectrode. Two-point calibration was carried out using air-purged and N₂-purged water samples prior to each use of the electrode. Dissolved oxygen concentrations were determined for surface water, the water-sediment interface, and then at 1-mm intervals into the sediment to a depth of 5 mm. Oxygen concentrations were determined prior to Eh and pH measurements to avoid disturbing or aerating the sediment. A depth of 5 mm was selected because DO concentrations were below the detection limit of 0.1% O₂ at deeper than 1 mm in all sediments. Eh and pH were then determined every 2 cm in the 15-cm long core using a calibrated Orion 250A meter with Eh (platinum, ±5 mV) and pH (Ag-AgCl, ±0.05 pH units) electrodes. The probes were inserted into pre-cut holes in the tube after removing the tape attached prior to sampling.

Sediment samples from Ziploc bags were freeze-dried in 74-mL vials and water loss (by mass) was determined to calculate porosity (ϕ, the fractional volume of water). Freeze-dried sediments were homogenized and then analyzed using the LOI methods of Heiri et al. (2001) to determine OM after heating to 550°C and CaCO₃ content after subsequent heating to 950°C. TOC and TN were determined using homogenized, freeze-dried sediments treated with 10% (v/v) hydrochloric acid to remove carbonate. Organic N was not lost during this process based on our laboratory testing. Analysis for TOC and TN was carried out using a Leco TruMac CNS/NS/Carbon/Nitrogen/Sulfur Analyzer following methods provided by the manufacturer (Version 1.3x, part number 200–753, August 2014). All values obtained for the certified reference material (Leco CRM 502-309) were within the 95% confidence interval (11.98 ± 0.44% TOC and 0.93 ± 0.04% TN). Laboratory precision was 6% for TOC and 2% for TN. Grain size was determined by the method of Folk (1974).

Samples of benthic organisms were stored at 4°C and sorted within 72 h of collection using an Olympus SZ Binocular Stereo Zoom Microscope with 7×–40× zoom magnification. Movement differentiated living vs. dead organisms within the 72-h window. Identifications were reported to the family level when genus and species could not be determined; all taxa were used to determine diversity. Abundances were reported as individuals per m² to facilitate comparisons with other studies. High abundances of tolerant organisms make sediments seem enriched with life, especially in stressed environments; however, this is not necessarily representative of the overall health of a benthic community (Warwick and Clarke, 1995). Therefore, community diversity was determined using the Shannon-Wiener diversity index (H′, Eq. 1). Effective Number of Species (ENS, Eq. 2) was used to express the number of equally occurring individual species needed in a community to produce the observed diversity index; values were reported as number of species instead of an index (Dauby and Hardy, 2012). ENS has been studied extensively for use in ecology; it has advantages over other diversity indices including versatility and direct relationship to species composition (Jost, 2006; Chao et al., 2010). Other diversity indices can be easily transformed into ENS values as shown for H′ (Eq. 2).

Community diversity for each station was calculated by first determining the Shannon-Wiener (H′) diversity index value:

where S refers to the species richness and p_i is the fraction of S of a single species (i) to the S of all species recorded. The ENS, also referred to as the Hill number (N₁, Hill, 1973), was calculated as:

Number of species per grab was recorded as richness (S). When biological variables such as diversity were subjected to regression analysis with our environmental data, the exponential form, ENS, had the best fit (most R² ≥0.5). Therefore, ENS values were used in this study as the main biological parameter for identifying relationships between environmental properties and benthic community composition.

Values for geochemical properties of sediments (e.g., TOC) were range standardized within sampling areas by converting the values to a 0–1 scale using:

A resemblance matrix was created from the standardized values using Euclidean distance. Hierarchical cluster analysis was used to explore groups of samples, with similarity profile analysis (SIMPROF) indicating data that could be permuted without any sample moving into a different group (PRIMER 6 and PERMANOVA). Similarities percentage analysis (SIMPER) was used to determine which geochemical parameters contributed to differences among groups.

Uncommon taxa (i.e., fauna present in <5 samples), were culled from the dataset. Faunal abundances were range standardized within sampling areas using Eq. 3 to balance the influence of less and more abundant groups. A permutation analysis of variance (PERMANOVA) that treated the groups identified from geochemical characteristics as a fixed factor was performed on a resemblance matrix created with a Bray-Curtis similarity measure and an arbitrary value of 1 added to eliminate undefined distances. Significantly different groups (p < 0.05) were analyzed with SIMPER to determine the taxa that contributed to statistically significant differences. To further identify relationships between biological and environmental variables, Pearson’s correlation coefficients (r) were used to compare values (e.g., TOC vs. ENS) and tested for significance using a two-sample, two-tailed t-test. Correlations were identified as mild for 0.2 ≤ r < 0.4, moderate for 0.4 ≤ r < 0.6, moderately strong for 0.6 ≤ r < 0.8, strong for 0.8 ≤ r < 0.9, and very strong for ≥0.9.

Concentrations of physical and chemical variables for the full suite of organic-rich sediments varied widely (Table 1). For example, silt + clay and sand content ranged from 6–99 to 1–84%, respectively. TOC and TN ranged from 1–8 to 0.1–0.7%, respectively. Therefore, environmental characteristics were compared using a hierarchical cluster analysis with SIMPROF. Four groups were identified as statistically different (Figure 2). Groups A and B were strictly center muddy deposits, whereas group C contained only perimeter muddy sand. Group D was the only mixed group with three perimeter deposits and one center deposit (Table 1); the anomalous center deposit contained 63.5% sand. SIMPER indicated that all environmental characteristics contributed to differences among groups and mean values for each parameter provided insight about differences in environmental conditions (Table 1).

Group A with the lowest sand content and highest concentrations of contaminants (i.e., TOC, TN, and H₂S) represented sandy mud, center deposits. Group C showed the opposite trend and represented muddy sand, perimeter deposits (Table 1). Samples from center deposits (group A) contained an average of 2.3-fold higher silt + clay contents than perimeter sediments (group C, Table 1). Groups B and D had similar values for porosity, sand and silt + clay that were between values for groups A and C (Table 1).

All sediments were anoxic, even within the top 1 mm. Redox potential (Eh) in the top 8 cm ranged from −43 to −176 mV in more sandy perimeter sediments relative to −98 to −190 mV in muddier center sediments. Negative Eh values showed the range of reducing conditions and supported high porewater concentrations of reduced nitrogen and sulfur as ammonium (NH₄⁺) and hydrogen sulfide (H₂S), respectively (Table 1). Average concentrations of porewater NH₄⁺ and H₂S were an overall ∼5- to 12-fold higher in the muddy central area sediments than in sandier perimeter deposits (Table 1). Seasonality also had an impact on these concentrations due to increased microbial activity during warmer, wetter summer months. Concentrations of NH₄⁺ and H₂S for all sediment porewater averaged ∼20 and ∼60% higher, respectively, during the wet season than the dry season (Table 1). Porewater H₂S values correlated strongly with NH₄⁺ (Figure 3A). pH as [H⁺] was about 37% lower in the center samples due, in part, to higher concentrations of dissolved H₂S and HS^– (Table 1; Peltzer et al., 2016).

Average values for TOC and TN in the top 5 cm of sediment from the center of muddy deposits were 1.9- and 2.2-fold greater, respectively, than for sediments from the sandier perimeter (Table 1). TOC values correlated very strongly with TN for all samples; the molar C/N ratio for all samples averaged 13.0 ± 2.6 and ratios were similar in the perimeter and center sediments (Figure 3C and Table 1). TN correlated positively with H₂S and silt + clay with lowest values for perimeter stations (Figures 3B,D). Average CaCO₃ content was not statistically different (p = 0.77) in perimeter vs. center locations because sizeable shell hash, when found, was present at both stations in the same area (Table 1).

The presence and abundance of key species varied considerably throughout our study area as well as between perimeter and center deposits. A total of 23,522 individuals and 35 species were collected from all organic-rich sediments (Table 2). The number of individuals ranged from 1 to 1,009 per grab with an average of 231 ± 285 per grab or 9,900 ± 12,300 per m². The organisms were divided into bivalvia, gastropoda and polychaeta, plus a combined group (others) (Table 2). Bivalvia was the most abundant class followed by gastropoda, others and polychaeta. Ranges in abundances for all four groups were large (Table 3).

A PERMANOVA showed that faunal abundances varied significantly among the four different environmental groups (pseudo-F = 7.29, PERMANOVA P-value = 0.001), and follow-up PERMANOVAs indicated that all four groups were significantly different from each other (P-values < 0.05). SIMPER identified the main taxa comprising these four groups (Table 4). Group A, which is most representative of center deposits, yielded ∼70% Macoma spp. and ∼14% Apolochus brunneus; whereas perimeter deposits, best represented by group C, had a greater mixture of taxa (Table 4).

Just three species, Macoma spp., Mulinia lateralis and Parastarte triquetra, comprised >99% of the bivalves and 70% of all individuals collected (Table 2). These three species showed somewhat different distribution patterns for perimeter vs. center deposits (Table 4). Macoma spp. was ubiquitous (39% of all individuals) and collected from all but one station (Table 2). Abundances of Macoma were 1.5-fold greater in perimeter, sandier sediments (group C) than center, muddy deposits (Table 4, organism B1). The bivalve M. lateralis comprised 20% of all individuals collected and was found at 12 of the 17 stations (Table 2). However, 25-fold more M. lateralis individuals were collected from perimeter deposits and the largest mean concentration was found in group D sediments (Table 4, organism B2). The third most abundant bivalve, P. triquetra (9% of all individuals; Table 2), also was predominantly found in more sandy perimeter deposits from group C (Table 4, organism B3). The other four bivalve species made up <0.4% of all bivalves (Table 2). Collectively, bivalve abundances were ∼3 times greater in the perimeter than the center of mud deposits (Table 3). Bivalvia abundance for all groups also was about 3.5-fold higher than the average for gastropoda (Table 3).

Gastropoda (sea snails and slugs) made up 20% of total abundance and included Acteocina spp., Crepidula fornicata and a gastropod in the family Naticidae (Table 2). Gastropod abundances were ∼4 times higher in the sandier perimeter than muddy center deposits, with the greatest number of gastropods in group D (Table 3). Acteocina spp. was the most abundant gastropod in this study; ∼6-fold more of these organisms were collected from sandier perimeter locations of group C than center deposits of group A (Tables 2, 4, organism G1). A species of Naticidae was almost as abundant as Acteocina spp. and found in 12 of 17 IRL areas (Table 2, organism G2). The abundance of snails from the family Naticidae was 2.4 times higher in perimeters than centers of muddy deposits (Table 4).

Although polychaetes were the least abundant (2%) class of organisms, they were ∼3 times more abundant in perimeter than center deposits with the greatest abundance found in group C (Tables 2, 3). The most abundant polychaetes were Diopatra cuprea, Schistomeringos pectinate and two species of Lumbrineris (spp. 2, 32 individuals per m²; spp.1, 28 individuals per m²) (Table 2). The only polychaete that was more abundant in the central area of the muddy deposits was Lumbrineris spp. 1 (organism P4) with ∼60% of the individuals collected in muddy center deposits of group A (Tables 2, 4).

The others category was made up of four species of malacostracans and a brittle star (class Ophiuroidea). This mixed group made up 8% of all the organisms collected. A. brunneus, an amphipod and the most abundant organism in the others group, was collected at 8 of 17 stations (Table 2). A. brunneus also was the fourth most abundant organism in this study (Table 2, organism O1). About 52% of the A. brunneus were collected from perimeter deposits (group C) relative to 38% from center deposits (group A), with the remainder from groups B and D (Table 4).

Pearson’s correlation coefficients were moderate to moderately strong between metrics for biota and environmental data (Table 5). Negative correlations for ENS vs. TOC, TN, LOI, H₂S, porosity and silt + clay were determined (Table 5). Thus, higher values for each of these six environmental covariates corresponded with less habitable sediments and lower ENS values. In contrast, positive correlations were obtained for ENS vs. sand content (Table 5). Gastropoda and polychaeta had the most significant correlations with environmental parameters, including a significant positive correlation with % sand. The bivalvia and others groups correlated significantly and positively only with CaCO₃ and bottom water temperature, respectively (Table 5).

We also investigated the significance of individual environmental variables in controlling ENS using correlation plots (Figure 4). The greatest difference for species abundance and diversity was for perimeter sediments from station TC vs. the center deposit at station EH (Figure 4). Station TC-P, with an ENS value of 7–9, had <0.2% TN, <1.5% TOC, and >80% sand. In contrast, station EH-C with an ENS of 1, had >0.7% TN, >6% TOC, and <10% sand (Figures 4A,B,E). Overall, higher ENS values were found where there was lower porosity and more sand (Figures 4C,E). Sediments from the IRL with the highest ENS had a porosity of 0.65–0.75 and 65–89% sand (Figures 4C,E). Sediment samples with >6% TOC, a porosity >0.9, and <30% sand, had ENS values <4 (Figures 4B,C,E).

As in the P-R model (Pearson and Rosenberg, 1978), data for abundance were plotted as a function of TOC (Figure 5). Lower abundances of polychaeta and gastropoda were found with increasing TOC (Figure 5A). This trend was not observed for the bivalves that had highest abundance at a TOC of 4–7% (Figure 5A). When ENS was plotted vs. TOC ranges, a clear trend of increasing ENS with decreasing TOC was found (Figure 5B). H₂S increased with increasing TOC (Figure 5B).

Increased loading of sediment OM can stress benthic infaunal communities (Cloern, 2001). This stress leads to shifts in community structure with opportunists becoming the predominant species amidst the loss of other organisms (Norkko et al., 2006). We can use our data plus global examples, to help improve our understanding of these shifts and identify opportunistic organisms that are adapted to fill the niche of organic-rich sediments.

Similarities percentage analysis showed that Macoma best characterized benthic communities in the areas sampled (Table 4). Macoma spp. was collected from all stations except EH-C, one of the least favorable habitats due to high TOC and TN (Table 2 and Figure 3C). Macoma are well suited for muddy and even anoxic environments due to their long siphons that can extend above the sediment surface (Seitz et al., 2001). Macoma have been reported to live in organic-rich effluent from a wastewater treatment plant (Hummel et al., 2000). Furthermore, Korpanty and Kelley (2014) noted that Macoma were better known as opportunists and capable of living among death assemblages, much like we found in shell hash layers. Therefore, their survival in organic-rich, muddy sediment in the IRL, was consistent with other studies.

Like Macoma spp., M. lateralis also is well adapted and commonly found in high porosity, anoxic, and sulfidic sediments because they can function (feed, digest, and grow) under such conditions with the specialized adaptation of anaerobic glycolysis to produce energy (Shumway et al., 1983). Their average abundance of 2,049 individuals per m² from our IRL study (Table 2) was much lower than the 74,000 individuals per m² reported for Hillsborough Bay, Florida, in the mid-1970s (Santos and Simon, 1980). Conditions in Hillsborough Bay, only 10% silt + clay and higher DO in bottom water, created more favorable environments than muddy deposits in the IRL (Santos and Simon, 1980). Unlike Macoma and M. lateralis that survived under poor conditions, P. triquetra have been found in sand or around seagrass beds in high densities (Virnstein et al., 1983). P. triquetra recolonized at a density of 6,600 individuals per m² in sediments with 90% sand in the adjacent Banana River Lagoon; however, their abundance decreased with decreasing DO and they were not considered opportunists (Grizzle, 1984).

Virnstein et al. (1983) found gastropods at greater densities in seagrass beds in the IRL than on bare, sandy substrate. The more robust gastropod communities in the TC area potentially survived due to proximity to a healthier benthic habitat where seagrass beds were present. TC is the only area that had no single dominant species contributing to the benthic community, but instead a more evenly distributed group of organisms indicating a greater diversity (Supplementary Table 1). Although gastropods made up only 20% of total abundance, they accounted for 46% of all species collected in this study (Table 2). Two of the most abundant gastropods, Acteocina and a snail from the family Naticidae, were much more abundant in sandier perimeter deposits in the IRL, as previously noted by Virnstein and Howard (1987; Table 4). Acteocina spp. have been reported in other lagoons and on mudflats such as the Bolsa Chica Lagoon (California), where ∼4,000 individuals per m² inhabited sediment with 3–83% sand and 3–12% OM, a habitat more like perimeter deposits in the IRL (Levin et al., 1998; Table 1 where LOI = OM). Much like Acteocina spp., the naticid also has been found in sandy, vegetated areas; however, because of its large foot, it is well suited to traverse soft sediments (Cernohorsky, 1971).

Polychaetes have been used as indicator species for environmental stress in benthic communities because they have a shorter life span and lower tolerance for contaminants (Dean, 2008). Low abundances for polychaetes in our study may highlight the stressful character of anoxic, muddy sediments in the IRL. However, D. cuprea (organism P1), our most abundant polychaete, has adapted by building tubes and completely replacing them within days, making them less likely to become smothered in soft sediments and die (Thomsen and McGlathery, 2005). Mangum et al. (1968) concluded that D. cuprea have an exceptionally low oxygen consumption and their abundance was not found to correlate with grain size. Fast tube building, the ability to survive in fine and course sediments and the capacity to cope with low oxygen concentrations make D. cuprea well suited for anoxic mud.

Unlike D. cuprea, Pectinaria gouldii create a cone-shaped tube from sand grains over a life cycle, not days, and their abundance is dependent on sediment type (Busch and Loveland, 1975). Therefore, they may become more easily buried or have lower preference for settlement in the muddy deposits found in the IRL than D. cuprea (P1). Less sand and higher water content in muddy sediments decreases sediment bulk density and organisms more easily sink into the organic-rich sediments and die. Hinchey et al. (2006) reported that when organisms sink and become buried by sediment beyond a manageable depth, respiration via the oxic water layer above the sediment was inhibited. The only polychaete in our study to be found more commonly in the sandy mud of center deposits was a Lumbrineris spp. Some species of Lumbrineris have been collected in the eastern Mediterranean Sea where they were found surviving in sandy mud habitats like our center deposits (Carrera-Parra et al., 2011).

More amphipods were collected during the wet season, most likely due to seasonal peaks of benthic amphipods as they have an inverse relationship to macrophyte biomass, which declines in hotter, wetter months (Stoner, 1980). A. brunneus, our most abundant amphipod, is typically found in sandy sediments and is exclusively a substrate deposit feeder (Conradi and López-González, 1999; de-la-Ossa-Carretero et al., 2012). However, A. brunneus (organism O1) were roughly equally found in perimeter and center deposits (Table 4). This lesser degree of selectivity may result from an overall greater mobility for amphipods than most organisms in the study; likewise, they could escape anoxic sediments after feeding.

The P-R model was formulated to describe benthic community responses to environmental stressors, such as increased OM or contaminants (Pearson and Rosenberg, 1978). This approach can monitor the overall health of an ecosystem because the benthic community is consistently present and can be studied temporally and spatially. Areas experiencing eutrophication also experience hypoxia and increased NH₄⁺ and H₂S; therefore, these chemical components have more recently been added to the P-R model (Thompson and Lowe, 2004; Hyland et al., 2005).

All 34 IRL samples had NH₄⁺ values that were less than the acute toxicity level of 1,950 μM (United States Environmental Protection Agency [USEPA], 1989). However, NH₄⁺ concentrations were greater than the USEPA chronic level of 250 μM in 26 of the 34 samples (Figure 3A). The 8 samples with less than chronic values were collected from sandier perimeter sediments. Seventy percent of the 34 sediment samples had H₂S values greater than the 4-day tolerance of 176 μM for Macoma (Figure 3A; Caldwell, 1975). Only bivalves and gastropods were found at stations with high concentrations of porewater H₂S and low ENS.

Pearson and Rosenberg (1978) also distinguished infauna as opportunists when only a few species contributed to overall abundances in a polluted environment. Only two taxa had a mean abundance >10 in the center deposits in group A with increased contamination (Table 4). In contrast, six taxa from group C had abundances >10. Group A also had the lowest mean sand content at ∼8% and highest TOC (∼6%) and TN (0.6%) relative to ∼55% sand, ∼3% TOC and 0.27% TN in group C, the main group of perimeter samples (Table 1 and Figure 4). One extensive study of TOC vs. species abundance and diversity across 7 coastal regions and over 900 stations around the world typically found higher species richness when sediment TOC values were ∼1% (Hyland et al., 2005). Conversely, lower species richness was found more frequently when sediment TOC values were >3.5% (Hyland et al., 2005). This observation by Hyland et al. (2005) is consistent with our observations that sediments with TOC ≥5.0% had ENS ≤5; whereas sediments with <2% TOC had ENS values >7, with richness as high as 24 species (Figure 4B). Burd et al. (2008) also found species richness declined as sediment TOC and TN increased. Elevated concentrations of TOC and TN yield an anoxic, degraded habitat that is more reducing with elevated concentrations of H₂S (Figure 3B). Biota from low ENS sediments were either 100% bivalves or gastropods, but not both. Even at our most prolific site (TC-P-dry), ENS values were ∼7 times lower than found by Mikkelsen et al. (1995) from a previous study of sandy sediments in the IRL.

Hierarchical clustering, along with the SIMPROF, and SIMPER analyses, show how the environmental parameters vary across the groups and in concert with each other (Figure 2 and Table 1). Correlation plots for ENS vs. environmental parameters help validate groupings from the omnibus statistical tests (Figure 4). When samples were grouped by environmental characteristics, the results supported distinct separation of sediments from center and perimeter locations (Figure 2 and Table 1). Groups B and D were similar in geomorphic data; however, there was a large difference in mean concentrations of H₂S and NH₄⁺. This difference was observed because group B selected exclusively wet season samples, when temperatures of bottom water were ∼10°C warmer and increased microbial activity increased concentrations of H₂S and NH₄⁺ in porewater from organic-rich sediments (Fox and Trefry, 2018). Group D selected samples collected during the cooler, dry season, and had a slightly greater sand concentration, consistent with lower porewater concentrations of H₂S and NH₄⁺ compared to group B.

Lower pH in center sediments reflected higher concentrations of dissolved CO₂ and H₂S. Even though concentrations of H₂S varied significantly between warm (wet) and cool (dry) seasons due to differences in temperature of bottom water, no significant differences were found for abundance, richness, diversity and ENS between the two sampling seasons (Table 5). Therefore, a 10°C difference in temperature had no discernible effect on abundance, richness, diversity or ENS at any given station. Overall, ENS was a valuable parameter for showing organismic responses to environmental variables and determining a hierarchy of habitable organic-rich sediments.

Groups C and D, mostly with samples from perimeter deposits, were more representative of transitional organic-rich sediments where more taxa contribute to total abundance (Table 4; Borja et al., 2000; Norkko et al., 2006). When abundance was compared with TOC, the bivalves follow a trend of a tolerant taxa as defined by the enhanced P-R model of Thompson and Lowe (2004); bivalve abundance was high while abundances of other taxa decreased (Figure 5A). In this study, the tolerant organism was Macoma spp. due to its high average abundance of ∼86 individuals per grab in group A with an average TOC of ∼6% (Tables 1, 4). This result was predicted by the enhanced P-R model of Hyland et al. (2005) which showed that increased OM had a negative impact on the overall health of benthic faunal communities. Sediments with more OM are well known to experience increased bacterial decomposition, decreased bottom water DO, and increased concentrations of H₂S, making these environments much less habitable for benthic life and more favorable for opportunistic, adapted species (Hansen et al., 2000; Schulz et al., 2000; Schultz and Urban, 2008). These data show that even at low OM content, proximity to muddy sediments can lead to decreased diversity due to potential for mud to mix with sand.

Fine-grained, organic-rich sediments have been accumulating in the IRL since the early 1950s, a consequence of increased human activity (Trefry et al., 1990; Trefry and Trocine, 2011; Fox and Trefry, 2018). These sediments help deplete oxygen in surrounding water, create less habitable benthic environments, store pollutants, and release N and P to the overlying water (Trefry et al., 1990). Although the total area of the lagoon covered with muddy deposits is estimated to be <10%, benthic fluxes of NH₄⁺ and phosphate to the overlying water column from these muddy sediments yield ≥25% of total fluxes of N and P to the IRL (Fox and Trefry, 2018; Tetra Tech Inc and Closewaters LLC, 2021). Dredging has been the primary form of remediation for muddy sediments in the IRL since the late 1990s. Such activity has focused on sediment-filled creeks and harbors to limit sediment transfer to the open lagoon. We believe that results from our study of species abundance and diversity in high-porosity, organic-rich sediments fill an information gap and thereby promote enhanced discussion of alternate forms of sediment remediation.

Our collection of organic-rich sediments from perimeters and centers of muddy deposits along 60 km of IRL provided a representative selection available infauna from sediments with a large range in concentrations of TOC, TN and porewater H₂S. Most sediments we sampled had H₂S concentrations that exceeded the LC₅₀ values for many infauna. Total abundance was ∼3-fold greater for more sandy perimeter deposits than for more muddy center areas. Richness from sandy sediments adjacent to mud were at least ∼7-fold lower than previously reported for sandy areas in the open IRL (Mikkelsen et al., 1995). Observed higher abundance, yet low diversity from this study, was due to greater abundances of a few bivalves, including the clam Macoma. ENS values from our IRL study ranged from 1 to 9 (i.e., H′ of 0–2.2). Species surviving in these deposits were adapted to survive extreme conditions, most likely by using a siphon as reported for Macoma species and other clams found in the muddy deposits (Hummel et al., 2000; Seitz et al., 2001). Negative correlations were obtained for ENS vs. TOC, TN, porosity, and porewater H₂S, whereas a positive correlation was obtained for ENS vs. % sand. Therefore, ENS could be predicted from correlation plots of ENS vs. various environmental variables. Center areas of muddy deposits had the lowest ENS values (<3); ENS increased to ∼9 with decreasing concentrations of TOC, TN, and H₂S, plus increasing sand content. These data for environmental and biological variables in muddy IRL sediments can be used to help predict the overall health of the benthic biota and determine where sediment treatment/remediation is needed to improve the overall health of the ecosystem. Post-treatment data also can help assess the success of remediation.

Dredging and capping contaminated sediments have been primary forms of sediment remediation globally for many years (National Research Council, 1989). Leaving contaminated sediments in place is generally more economically viable and often preferred due to concerns for contamination during removal or issues at containment and disposal sites on land. Sediment contaminants in IRL mud are primarily N and P that support large fluxes of dissolved ammonium and phosphate to the overlying water (Fox and Trefry, 2018). From 1998 to 2009, ∼2 million m³ of fine-grained, organic-rich sediment, along with some underlying sand, were dredged from Crane Creek, Turkey Creek and the St. Sebastian River. From 2016 to 2020, ∼600,000 m³ of sediment were dredged from Turkey Creek, the Eau Gallie River and some smaller areas (Fox and Trefry, 2018). All these locations were adjacent to or upstream of the IRL. Removal costs and recovery of dredged spoils have prompted discussion of alternative, in situ treatment of muddy sediment in the IRL. Capping of anoxic, fine-grained, organic-rich sediment with sand in Lake Worth Lagoon (LWL), Florida, has shown positive results (United States Army Corps of Engineers, 2018). Most sites where capping was carried out in LWL were close to shore or in dead-end passages. However, sand capping has been, thus far, deemed successful in LWL.

Sand is clearly a key variable that helps increase ENS in IRL sediments because it promotes decreased H₂S in porewater through lower TOC and TN (Figure 4). We now know that mixing of sand with anoxic, fine-grained, organic-rich sediments leads to a significant increase in ENS (Figure 4E). We believe the process of testing sand capping in open water areas should proceed carefully in barrier island lagoons with the best engineering designs and high quality, remote monitoring to track performance over time. Integration of sand capping, coupled with monitoring the boundary between sand and organic-rich sediments, would reduce nutrient loading to the lagoon system and improve water quality and benthic habitat.