Taghon, G.L.; Ramey, P.A.; Fuller, C.M.; Petrecca, R.F.; Grassle, J.P., and Belton, T.J., 2017. Benthic invertebrate community composition and sediment properties in Barnegat Bay, New Jersey, 1965−2014. In: Buchanan, G.A.; Belton, T.J., and Paudel, B. (eds.), A Comprehensive Assessment of Barnegat Bay–Little Egg Harbor, New Jersey.
Extended time series of estuarine benthic community composition and the chemical and physical properties of sediment are necessary for distinguishing natural variation from possible anthropogenic influences, such as eutrophication. In July 2012, 2013, and 2014, 97 stations, randomly located throughout the Barnegat Bay–Little Egg Harbor estuary, were sampled. Benthic invertebrates were abundant, and the community was, in general, highly diverse. Although there was considerable spatial variability in sediment-particle sizes throughout the estuary, overall the total organic carbon content of the sediments was low (<1%). Comparable historical data from 1965–2010 are spotty in spatial and temporal coverage, limiting comparisons to these recent data. Where comparisons can be made, the abundance and species composition of the benthos and the sediment properties, show few changes in 45 years. Despite high nutrient loading to this coastal bay, its shallow depth and general lack of stratification lead to relatively high dissolved oxygen levels, and it seems likely that heterotrophs in the sediments, both eukaryotes and prokaryotes, are rapidly metabolizing organic matter as it is produced.
The community structure of benthic macroinvertebrates is often used to assess habitat quality and to monitor the effects of human and natural stresses on marine and estuarine ecosystems (Borja, Dauer, and Grémare, 2012; Ellis et al., 2015; Grall and Glémarec, 1997; Hale, Cicchetti, and Deacutis, 2016; Levin et al., 2009). There are several reasons for that. Macroinvertebrates (operationally defined as individuals retained on screens with a given mesh opening (0.3, 0.5, or 1.0 mm) are abundant in most coastal and estuarine sediments, typically individual abundance is on the order of 103–104 m−2. Benthic communities are usually composed of many taxa from different phyla, and quantitative measures of community diversity (e.g., Rosenberg et al., 2004) and the relative abundance of taxa with different feeding behaviors (e.g., Pelletier et al., 2010; Weisberg et al., 1997) can be used to evaluate ecosystem health. Because most benthic invertebrates are sedentary as adults, they function as integrators, over periods of months to years, of the properties of their environment. Thus, changes over time can be used as surrogates of changes in environmental conditions over time, especially when independent measures of habitat quality are unavailable.
The Barnegat Bay ecosystem is potentially under stress from human impacts, which appear to have increased over the past several decades (Kennish et al., 2007). Most of the ecosystem impacts are thought to be related to excess nutrient inputs, leading to eutrophication. Barnegat Bay was classified as moderately eutrophic using the criteria in the U.S. National Estuarine Eutrophication Assessment (Bricker et al., 2007), but more recent assessments (e.g., Fertig et al., 2014) point to ongoing eutrophication. These criteria are water column chlorophyll a concentration, macroalgal blooms, dissolved oxygen concentration, nuisance/toxic algal blooms, and loss of submerged aquatic vegetation. Despite the long history of studies relating benthic macroinvertebrate community structure to habitat quality, such data have not been included in the assessment algorithm for Barnegat Bay. In part, this may be due to a lack of comprehensive data on the benthic communities, especially how they may have changed over time.
In this study, a data set from 97 sites in Barnegat Bay, sampled yearly in 2012, 2013, and 2014, was analyzed. This study was part of a comprehensive, multidisciplinary project to address data gaps on the ecological status of Barnegat Bay. The benthic community structure and sediment properties from these recent surveys were compared with historical data dating back to 1965, a longer period than that considered by Fertig et al. (2014).
Field and Laboratory 2012−14
Potential sampling locations were identified by overlaying a grid, with cell sizes of 400 × 420 m, on a map of Barnegat Bay, which resulted in 1566 potential sampling sites. A random-number generator was used to select 100 of the sites (Figure 1). Each site was sampled in July 2012, 2013, and 2014. Based on differential GPS coordinates, each site was re-located each year with an average precision of 63 m (standard deviation [SD], 49 m).
At each site, on each sampling date, two 0.04-m2 Ted Young Modified Van Veen grabs were taken for benthic infauna. Sediment was immediately sieved through a 0.5-mm mesh, and material remaining on the screen was fixed in a 3.7% formaldehyde solution in seawater, buffered with sodium borate, and containing Rose Bengal to stain organisms. Fixed material was transferred to 70% alcohol and sorted under low-power magnification. Organisms were counted and identified to the lowest practical taxonomic level, usually species. At each site, on each sampling date, one Van Veen grab was taken for sediment properties. The top 2-cm layer of sediment was removed, transferred to a stainless-steel bowl, and homogenized by stirring with a stainless-steel spoon. Subsamples of the homogenized sediment were taken for total organic carbon (TOC), total nitrogen, and total phosphorus, and for grain-size analysis. All sediment samples were placed on ice aboard the boat after collection and during transport to the laboratory, where they were stored at 4°C until analyzed. Sediment for elemental analysis was freeze dried and then ground with an agate mortar and pestle. TOC and total nitrogen were measured with a Carlo Erba NA1500 series 2 elemental analyzer, following Verardo, Froelich, and McIntyre (1990), except that inorganic carbonate was removed before analysis by exposing samples to concentrated HCl fumes overnight. Total phosphorus was measured using the ascorbic acid–molybdate method (Solórzano and Sharp, 1980). Two replicates were analyzed for each element, and the average value used in further data analysis. Reference sediment samples (Reference Material 2702—Inorganics in Marine Sediment, National Institute of Standards and Technology, Gaithersburg, Maryland, U.S.A) were regularly included in analyses, and measured values of all elements were within 90−110% of the reference values. Sediment particle sizes were measured on samples treated with hydrogen peroxide to remove organic matter, followed by wet sieving through a 0.063-mm mesh using distilled water with 10% sodium hexametaphosphate as dispersant. An aliquot of the <0.063-mm fraction was dried and weighed to measure silt and clay content. The >0.063-mm fraction was dried at 60°C, then separated with stacked, graded sieves, into >1.0 mm, 0.5−1.0 mm, 0.25−0.50 mm, 0.125−0.250 mm, and 0.063−0.125 mm fractions, for 6 min on a sieve shaker. Each size fraction was weighed. Grain-size statistics were calculated using the GSSTAT software program (Poppe, Eliason, and Hastings, 2004).
Field and Laboratory 1965−2010
Samples were taken in central Barnegat Bay, in the vicinity of Oyster Creek and Forked River (latitudes, 39°46′48″ N to 39°51′0″ N), over the years 1965−69 as part of a baseline survey of benthic communities before operation of the Oyster Creek nuclear generating station (Phillips, 1967, 1972). Stations were located by visual sighting on prominent landmarks onshore or with navigational aids, such as buoys and lights. Latitudes and longitudes of the Phillips (1967, 1972) stations were estimated by overlaying his maps with Google Earth images. Multiple types of sampling gear were used for benthos in that study, a Peterson grab (0.1 m2 area sampled) and two types of dredges. Sediment collected by all samplers was sieved over a 1.5-mm mesh. Because of the inconsistent sampling procedures and, especially the qualitative nature of the dredge samples, those data were only used for determining the presence or absence of macrofauna. Phillips (1972) took 73 sediment samples (Petersen and Ponar grab samplers, vertical penetration not provided) over the years 1965−69 at 23 locations. Sediments were separated into 1-ϕ increments. The raw data in the appendix of the Phillips (1972) dissertation were used to calculate median grain size and sorting coefficient using the software program GSSTAT (Poppe, Eliason, and Hastings, 2004). Phillips also measured total volatile solids on sediment samples collected in the summer of 1969, calculated as the percentage of mass lost on ignition at 700°C for 1 h. For comparison, loss on ignition was measured (500°C for 5 h) on samples from 18 of this study's stations located near the Phillips stations (mean distance between the Phillips stations and the stations in this study, 840 m; range 300−1600 m).
Sampling continued beyond the Phillips (1967, 1972) studies as part of monitoring studies related to the Oyster Creek nuclear generating station, resulting in a species list of benthic invertebrates found in that area of Barnegat Bay over the period 1965−74 (Jersey Central Power and Light, 1978). Loveland and Vouglitois (1984) continued to sample benthos at three of the original Phillips locations near Oyster Creek and Forked River in 1969−73. They used a 0.05-m2 Ponar grab and took seven grabs per site and sampling time, pooling the samples and sieving through a 1.5-mm mesh. The retained material was maintained at 5°C and sorted live within a week.
Haskin and Ray (1979) sampled the benthos at ∼39°37′48″ N on the western side of Little Egg Harbor in 1973 and 1974. They used a Ponar grab (0.05 m2) and sieved sediments through a 1-mm mesh. Most of their stations were in creeks and lagoons, but three in the vicinity of Cedar Run Creek and Dinner Point Creek were farther out in the bay and were within 700−1300 m of stations sampled in 2012−14 in the current study.
As part of the U.S. Environmental Protection Agency (USEPA) Environmental Monitoring and Assessment Program (EMAP), one location in Barnegat Bay was sampled for benthos and sediment properties each year in 1990, 1991, and 1993. More-comprehensive sampling took place in 2001, under the Regional Environmental Monitoring and Assessment Program (REMAP), when 80 locations were sampled. As part of the National Coastal Assessment (NCA) program, a varying number of locations were sampled in 2000 (3 locations), 2001 (15 locations), 2002 (3 locations), 2003 (5 locations), 2004 (9 locations), 2005 (3 locations), and 2006 (13 locations). Finally, under the National Coastal Condition Assessment program, five stations were sampled in 2010. All these samples were collected using a 0.04-m2 Young-modified Van Veen grab. For benthos, sediment was sieved through a 0.5-mm mesh. The top 2 cm of sediment were analyzed for TOC (USEPA, 1995). For grain sizes, only the percentage of sand (>0.063 mm) and the percentage of silt/clay (<0.063 mm) data are available. Data were accessed https://archive.epa.gov/emap/archive-emap/web/html/index-26.html (1990 data), https://archive.epa.gov/emap/archive-emap/web/html/index-27.html (1991 data), https://archive.epa.gov/emap/archive-emap/web/html/index-29.html (1993 data) https://archive.epa.gov/emap/archive-emap/web/html/index-124.html (2000−2006), and https://www.epa.gov/national-aquatic-resource-surveys/data-national-aquatic-resource-surveys (2010).
Benthic Community Analysis
Validity of each species identified was checked against the continuously updated list in the World Register of Marine Species (2017) or the Integrated Taxonomic Information System (2017). When necessary, species names in older data sets were brought up to date with their current taxonomic status, for consistency with currently acceptable standards (International Code of Zoological Nomenclature). Organisms from groups of typically smaller size (such as nematodes, ostracods, and harpacticoid copepods, etc.) were not enumerated. Organisms with a primarily epibenthic life position (such as caprellid amphipods, spirorbid and serpulid polychaetes, and mussels), and highly mobile taxa that are not reliably sampled with a Van Veen grab (such as decapod and mysid crustaceans) were not enumerated. Any taxon for which <10 individuals were collected in all samples combined from all data sets was omitted from further analysis. Stations at three locations farther upstream in Toms River sampled in 2012−14 were omitted from analyses because those locations have low salinities (average, 12.8) and are not representative of the bay.
Rarefaction curves were calculated based on the total number of individuals and number of taxa collected for a given sample or combination of samples. Hurlbert's E(Sn), the expected number of species for a given number of individuals (Hurlbert, 1971), was calculated using Primer version 6 software (Clarke and Gorley, 2006). For samples taken in 1990−2014, the number of each taxon per grab sample was available, allowing a straightforward calculation of total numbers collected. The older data are tabulated as abundances per square meter, but because the sample sizes and total number of samples collected were also given, the total number of individuals of each taxon collected each year at each site was recalculated to use in calculating rarefaction curves. Furthermore, because of differences in sieve-mesh sizes used to process samples collected in 1969−74 (1.0 or 1.5 mm) and 1990−2014 (0.5 mm), any taxon from the newer data sets that was not recorded in an older data set was omitted when making comparisons among years. This conservative step was necessary because of the possibility that a taxon was “missed” in the earlier studies because of a larger mesh size.
Multivariate analyses of community structure were performed using nonmetric multidimensional scaling (NMDS) on fourth-root–transformed Bray-Curtis similarity matrices of average abundances. Because only average abundance per square meter data were available for the older data sets, abundances in the newer data sets were expressed on a square meter basis when necessary; this conversion was unnecessary for among-year comparisons of the newer data sets. Similarity profile permutation tests (SIMPROF [Clarke, Somerfield, and Gorley, 2008]) were used to compare abundances for taxa combined at the Class taxonomic level. Sample groups were considered to belong to the same cluster if they were not different at a 5% probability level. For dissimilar groups, a similarity percentage test (Clarke, 1993) was used to identify the dissimilarity due to each Class.
In 2012−14, median grain size was highly variable with location and, at some locations, variable from year to year (Supplemental Figure S1). Generally, coarser sediments were present on the E side of the bay and near Barnegat Inlet and Little Egg Inlet, whereas finer sediments were present on the W side of the bay. There was no N–S trend. Two adjacent stations on the eastern side near Seaside Heights, stations 28 and 31, showed a large increase in median grain size in 2013, then a return to smaller sizes in 2014, similar to the sizes in 2012.
For the subset of the stations within 2000 m of the Phillips (1972) stations, median grain sizes were similar (0.138 mm in 1965−69 vs. 0.136 mm in 2012−14), but the range of grain sizes was narrower (0.054–0.429 mm in 1965−69 vs. 0.040–0.295 mm in 2012−14) (Figure 2). The difference between periods was marginally significant (Wilcoxon rank sum test, p = 0.1). The silt/clay content in 1965−69 and nearby stations in 2000−06 and 2012−14 were similar (Figure 2; Kruskal-Wallis one-way nonparametric analysis of variance [AOV], p = 0.34). The larger data set of bay-wide samples from 2000−06 and 2012−14 showed a slight increase in silt/clay content in the most recent samples (Figure 2), but the difference was not considered significant (Wilcoxon rank sum test, p = 0.28).
The sediment-sorting coefficient showed temporal and spatial trends in 2012−14 (Supplementary Figure S2). At the extremes, very poorly sorted sediments (sorting coefficient, >2) have a broad range of particle sizes and are often associated with nonselective transport processes, such as landslides; very well sorted sediments (sorting coefficient, <0.35) have a narrow range of particle sizes usually associated with selective transport processes, such as strong currents or wave action. At all locations, sediment became better sorted from 2012 through 2014. There was a latitudinal gradient, with sediment becoming better sorted from N to S. There was also more spatial variability in sorting from stations in N Barnegat Bay, whereas stations in the S were more uniform. The median value of the sorting coefficient of sediment was 0.88 in 1965−69 and 0.74 at nearby stations in 2012−14 (Figure 3). Sediments in 2012−14 were significantly better sorted (Wilcoxon rank sum test, p = 0.023).
TOC content of sediment in 2012−14 varied from 0.015% to 5.71% (Figure 4). Highest values occurred at locations in the north, near major freshwater sources (Metedeconk and Toms Rivers) and along the W margin, whereas lowest values occurred at locations near Barnegat Inlet and Little Egg Inlet and along the E margin of the Bay. There was little temporal variation at most locations. C:N ratios were mostly in the 8−15 range, with little spatial variability (Supplementary Figure S3). C:N increased at most locations in 2012−14. C:P ratios varied widely among stations in the northern part of the bay, from 30 to >200 (Supplementary Figure S4). C:P ratios were less variable, 10−140, at locations south of Toms River. C:P ratios throughout the bay generally did not vary among years.
Total volatile solids in sediments collected in 1969 had a median value of 3.38% dry mass (range, 1.61−6.85) (Supplementary Figure S5). Sediments from nearby locations in 2012−14 had a substantially lower median value of 1.44% (range, 0.23−5.78; p = 0.0006, Wilcoxon rank-sum test). The median value for TOC of sediments collected in 1990−2006 was 0.74% (range, 0.018−7.20) (Figure 5). Total organic content of sediments from 2012−14 was similar (median 0.57%, range 0.015−5.71) and was not significantly different (p = 0.2, Wilcoxon rank-sum test).
Benthic Community Structure
Total abundances of benthic macroinvertebrates in 2012−14 varied among stations by nearly 40-fold (Figure 6). The general trend was an increase from N to S, although abundances were low at several of the southern-most stations near Little Egg Inlet. There was no significant variation among years (Kruskal-Wallis one-way nonparametric AOV, p = 0.13). Taxon richness varied from 11 to 79 taxa sample−1 (Figure 7). Stations at and N of Toms River tended to have lower richness, with uniformly higher richness S of Toms River. There was no difference among years (p = 0.38). In all years, Polychaeta represented the greatest numbers of individuals and taxa, followed by Malacostraca (Figure 8). Although Polychaeta were still the numerical dominants in 2014, they were less abundant than in previous years, and the abundance of Malacostraca increased relative to previous years. Bivalvia and Gastropoda were less abundant in all years and had fewer species. The remaining taxa (Others) consisted of Clitellata, Pycnogonida, Ascidiacea, Anthozoa, Holothuroidea, Enteropneusta, and Nemertea.
Many more individuals were collected in 1969−73 by Loveland and Vouglitois (1984) than at nearby locations in 2001−14 because of the greater number of samples they took, but the number of taxa was similar in both periods (Supplementary Table S1). For similar numbers of individuals, more-recent samples from nearby stations had similar numbers of species as the older samples (Supplementary Figure S6). The NMDS analysis and SIMPROF identified six main groups of stations: group d, all stations sampled in 1969; group f, Stouts Creek stations in 1970−72; group b, Oyster Creek stations in 1970−73; group c, Forked River stations in 1970−73; group h, REMAP station 048 in 2001 and this study's station 23 in 2012−14; group i, this study's stations 47 and 48 in 2012−14 (Supplementary Figure S7). Three groups had only a single station each. A dissimilarity matrix based on similarity percentage (SIMPER) analysis of species abundances in the six main groups showed that the most dissimilar groups were b and h, b and i, c and h, and d and i (Table 1A). All these comparisons indicated an effect of time. For example, comparing b and h and b and i, Pectinaria gouldii was more abundant at Oyster Creek in 1970−73 than it was at nearby stations in 2001 and 2012−14, whereas Sabaco elongatus, Glycera americana, and Clymenella torquata were abundant in 2001 and 2012−14 but were not present in 1970−73 (Table 1B). Comparing d and i, Mulinia lateralis and Pectinaria gouldii were more abundant in 1969 than they were at nearby stations in 2012−14, whereas Elasmopus levis and Microdeutopus gryllotalpa were abundant in 2012−14 but were not present in 1969. There were also temporal differences within the 1969−73 timeframe, although the dissimilarity was less. Most notably, all three stations in 1969 were similar, whereas, in subsequent years, each station diverged from both the 1969 group and from each other. Mulinia lateralis and Pectinaria gouldii were extremely abundant in 1969 but then decreased drastically in 1970−73. Abundance of Ampelisca spp. was low in 1969 and increased in following years at all three stations but was especially high at the Stouts Creek station, indicating a distance effect; Stouts Creek was the northernmost station, approximately 2.4 km from Forked River and 4 km from Oyster Creek.
More individuals were collected in July and August of 1973 and 1974 by Haskin and Ray (1979) because of the greater number of samples taken then than at nearby locations in 2012−14, but the number of taxa was similar in both periods (Supplementary Table S2). For similar numbers of individuals, more-recent samples from nearby stations had similar numbers of taxa as the older samples (Supplementary Figure S8). NMDS and SIMPROF analysis of species abundances identified four groups of stations (Supplementary Figure S9). Group d consisted of all stations sampled in 1974 and a single 1973 station, groups a and b were each single stations sampled in 1973, and group c was two nearby stations sampled 2012−14. Groups c and d, representing the main effect of time, were 56% dissimilar (Table 2A). Of the 11 taxa contributing most to the dissimilarity, nine were less abundant or absent in 1973−74, relative to 2012–13 (Table 2B). Ampelisca spp. were the only taxa substantially more abundant in 1973−74.
Samples from 1990−2014 were collected from stations throughout Barnegat Bay and were processed similarly. For clarity in data presentation and analysis, the results are aggregated into five time “bins”: 1990−93 (3 stations total), 2000−02 (102 stations), 2003−06 (30 stations), 2010 (5 stations), and 2012−14 (97 stations each year). Given the differences in the numbers of stations sampled in each of those intervals, there was much variation in the number of taxa and number of individuals collected (Table 3). Rarefaction curves for species of Polychaeta and Malacostraca, the two classes that numerically dominated the benthos, show that sampling efforts were sufficient to capture taxa likely to be present in the Bay in 2000−14 (Figure 9). Polychaete species diversity and malacostracan species diversity were greater in 2012−14 than they were in any other period. The rarefaction curve for Bivalvia indicated that species diversity was greatest in 2003−06, slightly lower in 2012−14, and still lower in 1990−93. Gastropod species diversity was also lowest in 1990−93. The species richness curve for 2000−02 samples rose quickly, then leveled, whereas the curve for 2012−14 did not plateau.
Between 1990 and 2014, 243 locations were sampled using consistent methods, some more than once. Most of those stations, 199 (82% of all), were sampled in two periods 2000–02 and 2012–14. SIMPROF analysis of the species abundance for the four most-abundant classes at those stations identified 12 groups. Four of these groups, which comprised 131 of the 199 locations, were selected for SIMPER analysis to identify temporal differences and spatial differences, especially associated with latitude and the strong salinity gradient in Barnegat Bay.
Group E consisted mostly of stations from 2000−02 that were mainly south of latitude 39°57′0″ N, the region of the bay where salinity increases from typical values of 20 further N to 28 S of latitude 39°47′59.9994″ N (Figure 10). Only two stations from 2012−14 were in group E. Group G also consisted of stations south of latitude 39°57′0″ N, mostly from 2012−14 (52 stations) and only six from 2000−02. Most of the dissimilarity in these two groups arose from the much greater abundance of Malacostraca at these transitional to high salinity locations in 2000−02 and, secondarily, from the greater abundances in 2012−14 of Bivalvia, Gastropoda, and Polychaeta (Table 4).
Stations in the northern part of the bay, where average salinity is 20, also showed a difference over time. Group L consisted only of stations from 2000−02, most at or north of latitude 39°57′0″ (Figure 10). Group I consisted mostly of stations from 2012−14 in the north; there were only two stations from 2000−02 in that group. Most of the dissimilarity in these two groups was due to the low abundances of Gastropoda and Bivalvia in 2000−02, compared with their moderate abundances in 2012−14 (33 and 18 sample−1, respectively) (Table 4).
There was a strong latitude effect in 2000−02. Stations mainly in the middle and south sections of the bay (Figure 10) had much greater abundances of Malacostraca than stations had in the northern part of the Bay (∼16 times greater), and moderate abundances of Gastropoda and Bivalvia, whereas these two Classes were in low abundance in the north (Figure 10, Table 4). There was also a latitude effect in 2012–14, but it was not as strong as in the earlier samples. Malacostraca were again more abundant in the mid and S Bay than it was in the N, but only about sixfold.
Change in benthic community structure over time is commonly used to evaluate both natural and anthropogenic changes in coastal environments. There are difficulties in comparing results among studies that occurred decades apart. Coastal benthic community structure is strongly affected by seasons, especially in temperate zones. In this study, data used were, for the most part, from samples collected in July–August, with the exception of the Loveland and Vouglitois (1984) data, which were annual averages. Variations in sampling gear and sample processing hindered comparisons. Different types of grab samplers recover different amounts of sediment, both in terms of surface area and depth. A standard practice when reporting results is to “normalize” species abundances to an area of 1 m2, although that is never the actual surface area sampled at any single time and place. Such extrapolations to a larger area are probably never justified; distributions of benthic invertebrates are notoriously patchy on spatial scales <1 m2 (e.g., Eckman, 1979; Hogue, 1982). Here, because the older data from 1969−74 were as abundances in square meters, the data from 1990 onward had to be converted to the same units when comparing the older and newer data. Once sediment samples were collected, there can also be inconsistencies in processing, usually in which mesh size is used to separate organisms from sediment. This is a well-known problem in benthic ecological studies (Reish, 1959), but there is no obvious way to correct for it. Here, the older data (1969−74) were based on mesh sizes of 1.0 mm or 1.5 mm, whereas the newer data are based on a 0.5-mm mesh. The conservative approach used here, when comparing data sets that were based on different mesh sizes, was to drop from analyses any species that were not collected with the larger mesh size. Finally, consistent identification of species is critical. Expertise of those who are identifying species varies widely. Without voucher specimens, often not available even for more-recent surveys, it is impossible to assess whether species were consistently identified (Gallagher and Grassle, 1997).
Species accumulation (rarefaction) curves for stations sampled in 1973 and 1974 by Haskin and Ray (1979) (Supplementary Figure S8) illustrate the problem of comparing samples that were processed using different sieve-mesh sizes. When only taxa found by Haskin and Ray are considered, species accumulation curves for most of the older and the two recent, nearby stations rose at similar rates initially, but with increased sample sizes, there were ultimately more species present in 1973−74 than in 2012−14. Summing all the stations and years for the two periods shows there were 46 taxa present in 1973−74 vs. 29 in 2012−14 (Figure S8). Based only on the criterion that stressful environmental conditions lead to lower species richness, that result could lead to the conclusion that conditions at this location in Barnegat Bay had declined. A different picture emerges if all individuals collected at that location in 2012−14 are considered. The species accumulation curve rises more rapidly, is always above the 1973−74 curve, and continues to rise, suggesting that more species would have been found if more samples were taken. It is not possible to conclude whether this result is due to the smaller mesh size used in 2012−14 or to an actual increase in species richness. The abundance data for the individual species is subject to the same uncertainty for many, but not all, species (Table 2). Ampelisca abdita was more abundant in 1973−74, which could not have been an artifact of different mesh sizes; population dynamics of Ampelisca spp. are discussed below. Species more abundant in 2012−14 could be due to the smaller mesh size used, but this is unlikely for all cases. The maldanid polychaete Sabaco elongatus was rare in 1973−74 but was common in 2012−14. That large species would likely have been retained on the 1-mm mesh used by Haskin and Ray. The same argument applies to the gastropods Turbonilla interrupta and Acteocina canaliculata. This study concludes, conservatively, that species diversity in 2012−14 was similar to that found in 1973−74 (and probably greater). Beyond the species list, the community structure was clearly different between the two periods (Supplementary Figure S9). The community in 1973−74 was dominated by the amphipod Ampelisca abdita, whereas species abundances in 2012−14 were more uniform.
There are also difficulties comparing the data in Loveland and Vouglitois (1984) with more-recent data. They sampled three stations over the 5-year period from 1969−73. Unfortunately, in addition to the difference in mesh size, the data for species abundances in their article include only the 29 taxa that were in the top-10 most abundant for at least one of the 15 data sets (three stations × 5 y). There were many more samples and taxa collected near their stations during this monitoring program in the central part of Barnegat Bay. Phillips (1972) lists 98 taxa as occurring there, but there are no quantitative abundance data for most of them. Based on analyzing only the 29 taxa for which quantitative data are available, the species accumulation curves for stations in 2012−14 are similar to those from 1969−73 at two of the three stations (Supplementary Figure S6). The Forked River location was more diverse for three of the years. Summing all the stations and years for the two periods shows that there were 29 species present in 1969−73, but only 23 of them were collected in 2001−14 (Supplementary Figure S6). Again, however, if all taxa collected are considered, there were 143 present in 2001−14, closer to the total number of 98 collected by Phillips (1972).
The Loveland and Vouglitois (1984) results exemplify that changes in benthic communities can occur from one year to the next, in the absence of any obvious environmental change. The community structure at their stations in 1969 was nearly as dissimilar to the same stations in 1970 as to nearby stations sampled 32 to 45 years later (Supplementary Figure S7, Table 1). Much of this difference was due to the high abundances of the bivalve Mulinia lateralis and the polychaete Pectinaria gouldii in 1969 and its decline in subsequent years. In contrast, Ampelisca spp. were least abundant in 1969.
Samples from 1990 onward were all collected and processed similarly, making comparisons more straightforward. Diversity in samples collected in 2012−14 was greater than for earlier collections, especially for species in the classes Polychaeta and Malacostraca (Figure 9). Although the lower diversity in the previous periods 1990−93 and 2003−06 could be due to fewer locations sampled at those times, that was not the case for the 2000−02 period when sample density and spatial coverage were similar to 2012−14. None of the “new” species collected in 2012−14 are invasive species; in fact, 21 of them were previously recorded in Barnegat Bay in 1965−69 by Phillips (1972). The diversity of benthic invertebrates throughout Barnegat Bay in 2012−14 was greater than it has been since 1990 and, most likely, since 1969.
Evaluating changes in community structure over time is complicated by the strong salinity gradient in Barnegat Bay. North of latitude 39°53′59.9994″ N, salinity is low (15−20), because of the two major sources of freshwater to the Bay, Toms River and the Metedeconk River. South of latitude 39°53′59.9994″ N, salinity increased to 25–30 because of the oceanic influence through Barnegat Inlet and Little Egg Inlet. The significant change in community structure between 2000−02 and 2012−14 south of latitude 39°53′59.9994″ N (Figure 10) was due mainly to the nearly 10-fold decline in the abundance of malacostracans, which were the dominant class in 2000−02. There were 58 species of Malacostraca in Barnegat Bay, but the numerical dominants are three species of Ampelisca, which comprise 51% of all individuals in this Class. Polychaeta, in contrast, increased in abundance south of latitude 39°53′59.9994″ N and took over as the numerical dominants in 2012−14. There were 104 species of polychaetes, and that Class was dominated by one species, Mediomastus ambiseta, which accounted for 38% of all individuals.
The relative abundances of Mediomastus ambiseta and Ampelisca spp. have fluctuated in Barnegat Bay since 1990 (Figure 11). Until 2006, Ampelisca spp. were more abundant at most locations, and then by 2012 and 2013, there was a shift to dominance by M. ambiseta. In 2014, there was a shift back to similar abundances. Those two taxa are important components of benthic communities throughout the Virginian Province. Ampelisca spp. are well known for substantial swings in abundance between years (e.g., abundance in central Barnegat Bay declined precipitously between 1969 and 1970). Mediomastus ambiseta may be a relative newcomer to the Virginian Province, but that is debatable. The species was not recorded in Barnegat Bay in 1965−74. It does not appear to have been present in Narragansett Bay, Rhode Island, before 1975, but since then, it has become a dominant species in the middle regions of that bay (Frithsen, 1989). Frithsen (1989) discounts the possibility that M. ambiseta was present earlier but not recorded because of differences in sampling methods, sieve sizes (large specimens are up to 20 mm long), or taxonomic expertise, instead suggesting its rapid establishment was a result of greater organic enrichment. On the other hand, the species was not described in the New England area until 1971 (Hobson, 1971). Gallagher and Grassle (1997) suggest that M. ambiseta may have been present previously but was misidentified as Heteromastus filiformis, another capitellid polychaete. Heteromastus filiformis was collected by Haskin and Ray (1979); unfortunately, there are no voucher samples to check their identification. Kinner and Maurer (1977) collected both M. ambiseta and H. filiformis in Delaware Bay in 1974. It is probable that M. ambiseta was present in Barnegat Bay before 1990. Mediomastus ambiseta is now the most abundant macrobenthic species in the entire Virginian Province (Gallagher and Grassle, 1997).
The presence and abundance of various benthic taxa have been used as indicators of coastal environment quality, in particular, the effects of eutrophication (e.g., Pearson and Rosenberg, 1978). Benthic taxa are placed into groups reflecting the degree to which they are adversely affected by, or can tolerate, the effects of eutrophication, usually expressed as organic enrichment of sediment and related effects, such as low dissolved oxygen concentration. A widely used framework is to place taxa into one of five ecological groups (e.g., Grall and Glémarec, 1997). For the eastern United States, of the three species of Ampelisca, A. verrilli is classified as group 2 (“indifferent to enrichment, always present in low densities with nonsignificant variations in time”) and A. abdita and A. vadorum are classified as group 3 (“tolerant of excess organic matter… may occur in normal conditions but their populations are stimulated by organic enrichment”) (Gillett et al., 2015). Populations of Ampelisca spp. rapidly expanded in Boston Harbor after reductions of organic loading when wastewater outfalls were moved offshore. In Narragansett Bay, populations of Ampelisca spp. expanded between 1988 and 2008, attributed, in part, to improved sewage treatment and reduced organic loading (Shumchenia, Guarinello, and King, 2016). Classification of Mediomastus ambiseta varies among studies. Weisberg et al. (1997) and Llanso et al. (2002) classified that species as pollution sensitive, whereas Dauer, Rodi, and Ranasinghe (1992) classified it as an opportunist, and Gillett et al. (2015) placed it into group 4 (“second-order opportunistic species, small species with a short life cycle; adapted to life in reduced sediment”). Thus, variations in the abundances of those two common and abundant taxa would be expected to be related to sediment organic carbon concentration. There is little evidence for that, however, in Barnegat Bay. Both taxa occur over nearly a 1000 ranges of sediment organic carbon concentration (Figure 12). Peak abundances were at <1% sediment organic carbon. Furthermore, the cumulative abundance of four species of maldanid polychaetes showed a similarly broad distribution with a peak at 1% sediment organic carbon (Figure 12). These maldanids are all classified in group 1 (“species very sensitive to organic enrichment”) and increased greatly in abundance between 1990−2010 and 2012−14.
A threshold of 1−2% sediment organic carbon concentration is typically used to indicate “good” conditions; at greater concentrations, benthic communities are negatively affected (Hale, Paul, and Heltshe, 2004; Hyland et al., 2005; Pelletier et al., 2010). By that criterion, most of Barnegat Bay was in good condition in 2012−14 and was better than in prior years. Sediment TOC was low and generally only exceeded 2% at stations adjacent to rivers and creeks, where fine-grained sediments are transported into the bay. Sediment TOC was lower, on average, in 2012−14 than in 1990−2006 (Figure 5), and sediment concentrations of total volatile solids was substantially lower in 2012−14 than it was in 1969 (Supplementary Figure S5). One problem with relating abundance or diversity to sediment TOC is that TOC is a “relict” variable; it represents the organic matter not metabolized by all heterotrophs in the sediment, eukaryotes as well as prokaryotes. Whether a value measured at a given time represents a steady state between supply and metabolic breakdown or will decrease over time as organic matter is respired and converted to biomass or will increase if supply of newly deposited organic matter exceeds use is unknown. In Barnegat Bay, the high abundances of infauna and the low sediment carbon content point to efficient cycling of organic matter under aerobic conditions. Average bottom-water dissolved oxygen concentrations were 5.8 mg L−1 in 2012, 6.2 mg L−1 in 2013, and 6.7 mg L−1 in 2014 (unpublished data), consistently >5 mg L−1, the threshold concentration typically used to indicate “good” conditions (Hale, Paul, and Heltshe, 2004; Pelletier et al., 2010). The shallow water depth and lack of density stratification in Barnegat Bay lead to high dissolved oxygen levels for the benthos. Sediment C:N ratios were consistently >6−8, indicating that high-quality organic matter was not accumulating in the sediments but was aerobically degraded. High aerobic activity was also indicated by the low sediment C:P ratios, usually well below the Redfield 106:1 ratio. Under aerobic conditions, phosphorus released from organic matter decomposition is bound to iron and manganese oxides and stored in sediments.
There are numerous difficulties in evaluating whether, and how, benthic community structure responds to changes in the environment. Lack of long-term measures of environmental properties; inadequate understanding of the effects of ecological processes, such as competition and predation on community structure; and the present inability to incorporate the role of climate regime shifts on local community structure, all contribute to these difficulties (Nichols, 2003). All are in urgent need of further study. A meta-analysis by Gamfeldt et al. (2015) concluded that there is a generally positive relationship between community diversity and ecosystem functioning. In Barnegat Bay, macrofaunal species diversity in 2012−14 was as high as, or higher than, it was in the past. Benthic community structure should be included in assessing human effects on estuaries.
Compared with previous data dating back to 1969, benthic macroinfauna in Barnegat Bay over the period 2012−14 was generally at least as abundant throughout the bay. Species diversity in the most recently collected samples was as high, or higher, than it was in earlier samples. Polychaete worms and malacostracan crustaceans continued to dominate the benthos, with lesser contributions from gastropods and bivalves. Sediments in 2012−14 were low in organic matter, in some cases, lower than reported in 1969. Sediment particle size varied considerably with location in the bay, but there were no apparent long-term trends. Taken together, these data indicate little change in the composition and functioning of the benthic community in Barnegat Bay during the past 49 years. Barnegat Bay harbors an abundant, diverse benthic community that appears to efficiently metabolize organic matter and prevent its accumulation in bottom sediments. In addition to measures of nutrient concentrations and planktonic processes, the structure and function of benthic communities should be included in assessments of human effects on estuaries.