Biological Abundance and Diversity in Organic-Rich Sediments From a Florida Barrier Island Lagoon

Organic-rich sediments in estuaries and the coastal ocean are often a product of land clearing, runoff of excess nutrients and other human activities. They can harbor pollutants, oxygen-consuming microbes and toxic hydrogen sulfide (H2S), thereby creating a hostile environment for infauna. In one barrier island lagoon, the Indian River Lagoon (IRL), Florida, layers of organic-rich sediments have increased substantially in thickness and areal extent over the past 60 years. Geochemical properties of these muddy sediments have been described; however, less is known about their habitability. We analyzed infauna and geochemical properties of 102 samples taken during wet and dry seasons at 17 locations spanning 60 km of the lagoon. We quantified infaunal abundance and diversity (Shannon-Wiener, H′) and determined Pearson’s correlation coefficients for effective number of species (ENS = eH′) vs. sediment porosity (ϕ = 0.69–0.95), organic carbon (1–8%), nitrogen (0.1–0.7%), silt + clay (16–99%), porewater H2S (5–3,600 μM), and other environmental variables. Small bivalves accounted for 70% of the organisms collected, followed by gastropods, polychaetes and other biota. The bivalves were predominantly Macoma spp., Mulinia lateralis and Parastarte triquetra with average abundances of 3,896, 2,049, and 926 individuals per m2, respectively. High abundance of some species, such as Macoma, showed that these opportunists had adapted to poor quality sediments. More than two-thirds of the 35 species collected were present at <100 individuals per m2 of sediment. Cluster analysis identified four groups of stations with significantly different geochemical properties. Permutation analyses of variance indicated that the four groups also represented statistically different infaunal communities. Diversity decreased with increasing sediment concentrations of organic carbon, nitrogen and silt + clay; however, community richness at our most prolific station along the perimeter of muddy deposits was ∼7 times lower than found previously in sandy sediments from the IRL. The results identified areas where infaunal communities have experienced the greatest stress due to accumulation of organic-rich sediments. Results from this study help support management plans for remediation of organic-rich mud and improvement of sediment and water quality, especially in areas identified with low ENS.


INTRODUCTION
Estuarine eutrophication has long prompted a global call for action by environmental scientists and managers (e.g., Rhyther and Dunstan, 1971;Paerl, 2006). Impacts from eutrophication in estuaries and the coastal ocean include nutrient loading, disruption of marine ecosystems and increases in sediment organic matter (OM) (Zimmerman and Canuel, 2000;Andersen et al., 2006;Lu et al., 2020). Organic-rich sediments harbor contaminants, store nutrients, and create anoxic benthic environments (National Research Council, 1989;Southwell et al., 2010). Increased deposition rates for OM in sediments also alter microbiomes and induce significant shifts in sediment chemistry, including releases of dissolved nutrients to the overlying water column (Cloern, 2001). Management strategies that improve water quality and ecosystem health require a better understanding of connections between environmental variables and biological abundance and diversity.
Human impacts within watersheds, including agriculture, land clearing and urban development, can increase sedimentation rates in estuaries and lead to progressive increases in nutrient loading and eutrophication (Cronin and Vann, 2003). Such impacts create poor water quality, loss of seagrasses, decreased fisheries and increased potential for harmful algal blooms (HABs) and fish kills (Unsworth et al., 2018). Estuaries, especially barrier island lagoons, are subject to these effects because they often have limited exchange with the coastal ocean. Therefore, future management decisions regarding remediation of organic-rich sediments should be based on a detailed understanding of how and why benthic biota respond to these deposits.
Bacterial decomposition of OM in organic-rich sediments creates anoxia and toxic concentrations of dissolved ammonium (NH 4 + ) and hydrogen sulfide (H 2 S) (Gu et al., 1987;Hitchcock et al., 2010). Exposure to elevated concentrations of NH 4 + and H 2 S greatly shorten times for the onset of adverse health effects and mortality in benthic communities (Kemp et al., 2005;Long and Seitz, 2009). Anoxia or a low concentration of dissolved oxygen (DO) has been linked to reduced growth rates and biomass for many marine organisms; low DO also is lethal to macrofauna when concentrations remain <1.4 mg/L for 1-2 weeks (Josefson and Widbom, 1988). In Chesapeake Bay, increased periods of low DO were strongly and negatively correlated with the integrity of benthic communities (Dauer et al., 2000).
Multiple investigations have identified natural and anthropogenic impacts, including increased OM and sedimentation rates, that alter benthic communities by limiting diversity and allowing more aggressive, opportunistic organisms to excel (Pearson and Rosenberg, 1978;Alongi, 1989;Norkko et al., 2006). As sandy sediment becomes buried and replaced by higher porosity, organic-rich sediment, benthic community abundance and diversity decrease because biota are unable to cope with OM degradation and resulting hypoxia/anoxia (Chou et al., 2004;Hale et al., 2016). Various forms of the Pearson-Rosenberg (P-R) model have shown that decreasing numbers of taxa, abundances of organisms and biomass are functions of increased OM and contaminants (Pearson and Rosenberg, 1978;Thompson and Lowe, 2004;Hyland et al., 2005).
Inputs of OM from rivers have, in some cases, increased macro-faunal abundance; however, when sediment total organic carbon (TOC) exceeded 1.0-1.5%, the number of taxa was reduced (Burd et al., 2008). Numbers of taxa and abundances also decreased when concentrations of sediment total nitrogen (TN) exceeded 0.4% (Burd et al., 2008). One study conducted in a subtropical estuary near Hong Kong found that the diversity of the microbenthic communities correlated negatively with TOC, total Kjeldahl nitrogen, NH 4 + and total phosphorus (TP) in sediments (Shin et al., 2008). Sediment type and grain size also have been shown to be major factors affecting community diversity and abundance (Lin et al., 2018). For example, Chou et al. (2004) concluded that family richness and abundance decreased in benthic habitats with 10-65% increases in silt + clay.
Our study focused on organic-rich sediments in the Indian River Lagoon (IRL), a barrier island lagoon in Florida. Organicrich sediments in the IRL typically contain >60% silt + clay, >75% water by weight, >10% OM, and >100 µM dissolved H 2 S (Trefry and Trocine, 2011;Fox and Trefry, 2018). These sediments release dissolved nutrients that fuel algal blooms (Fox and Trefry, 2018). Along with increased algal blooms, the IRL has experienced reductions in seagrass, decreased water quality and fish kills (Phlips et al., 2015). The goal of this study was to improve our understanding of the link between environmental variables and the abundance and diversity of benthic biota in organic-rich sediments in the IRL. Such knowledge supports management decisions related to successful remediation of finegrained, organic-rich sediments in the lagoon and elsewhere.

Study Area
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 2 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.

Field Sampling
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 4 + and H 2 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 3 , TN, TOC, and porosity. All samples were kept cold (∼4 • C) until analyzed or freeze-dried.

Physical and Chemical Analysis of Sediments
Four mini-cores from each site were placed into N 2 -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 2 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 2 gas before being analyzed for NH 4 + to limit analytical interference from H 2 S. Concentrations of NH 4 + were determined colorometrically using the standard phenate method #4500-NH 3 (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 2 -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 2 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). Freezedried 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 3 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).

Analysis of Sediment Infauna
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 2 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 1 , 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 2 ≥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.

Statistical Analysis
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: Range standardized value = Raw value − Minimum value for area Maximum value for area − Minimum value for area (3) 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.

Sediment Physical and Chemical Variables
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 TABLE 1 | Similarities percentage analysis of mean environmental characteristics of the four groups identified from hierarchical cluster analysis with the rank order from highest concentration to lowest for each parameter per group.

Parameter
Group mean ± Standard deviation Rank order A (Ctr) (n = 13) B (Ctr) (n = 2) C (Perim) (n = 13) D (Perim, Ctr)* (n = 4) A B C D Deposit thickness (m) Standard deviation added after SIMPER and for variables without a rank order. Seasonal average for groups A and C recorded for dry (D) and wet (W) when water temperatures near the sediment was 19 ± 3 • C and 29 ± 2 • C, respectively. Groups identified as either center (ctr), perimeter (perim) or both depending on the samples they contained. *3 perimeter, 1 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 2 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 4 + ) and hydrogen sulfide (H 2 S), respectively ( Table 1). Average concentrations of porewater NH 4 + and H 2 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 4 + and H 2 S for all sediment porewater averaged ∼20 and ∼60% higher, respectively, during the wet season than the dry season ( Table 1). Porewater H 2 S values correlated strongly with NH 4 + ( Figure 3A). pH as [H + ] was about 37% lower in the center samples due, in part, to higher concentrations of dissolved H 2 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 2 S and silt + clay with lowest values for perimeter stations (Figures 3B,D). Average CaCO 3 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).

Variation in Benthic Communities
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 2 . 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 FIGURE 2 | Hierarchical cluster analysis created from range-standardized environmental characteristics and Euclidean distance. Four groups identified as shown by samples above red bars (A-D). Samples from center (C) locations are highlighted, but not perimeter (P) samples. Sampling seasons represented by wet and dry. deposits, best represented by group C, had a greater mixture of taxa (Table 4).
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 2 ; spp.1, 28 individuals per m 2 ) ( 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).

Variation in Benthic Communities Related to Environmental Conditions
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 2 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 3 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 2 S increased with increasing TOC (Figure 5B).

Adaptations Support Opportunistic Organisms in Muddy Indian River Lagoon Sediments
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 2 from our IRL study ( Table 2) was much lower than the 74,000 individuals per m 2 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 2 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 2 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 5 | Pearson correlation coefficients (r) for significance between biotic and sediment factors for core/grab data from 34 samples collected at 17 locations in the Indian River Lagoon; abundances of organism groups, richness (S), Shannon-Wiener diversity index (H ), and effective number of species (ENS). 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.

Physical and Chemical Controls on Abundance of Benthic Infauna
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 4 + and H 2 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 4 + values that were less than the acute toxicity level of 1,950 µM (United States Environmental Protection Agency [USEPA], 1989). However, NH 4 + 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 2 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 2 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 2 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 2 S and NH 4 + . 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 2 S and NH 4 + 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 2 S and NH 4 + compared to group B.
Lower pH in center sediments reflected higher concentrations of dissolved CO 2 and H 2 S. Even though concentrations of H 2 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 2 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.

Applications to Lagoon Management
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 4 + 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 2 S. Most sediments we sampled had H 2 S concentrations that exceeded the LC 50 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 2 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 2 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 3 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 3 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 2 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.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
KF and AF: planning, sampling, analysis, data interpretation, and manuscript preparation. JT: planning, sampling, data interpretation, manuscript preparation, and funding. CJ: analysis, data interpretation, and manuscript preparation. All authors contributed to the article and approved the submitted version.