Sediment Diatoms as Environmental Indicators in New Jersey Coastal Lagoons

Shallow coastal lagoons or lagoonal estuaries are dominant landscape forms of the mid-Atlantic coastline. These are exceptionally productive ecosystems highly valued for providing habitat to fish and shellfish, storm protection, and recreational opportunities (Nixon, 1982). Eutrophication and climate change are major threats to their ecological integrity, affecting both ecosystem structure and function (Anthony et al., 2009). The impact of increasing nutrient loading is especially strong in poorly flushed lagoons with highly developed watersheds, such as Barnegat Bay–Little Egg Harbor (BB-LEH) estuary (Baker et al., 2014; Defne and Ganju, 2014). In the last two decades, Barnegat Bay has experienced many of the negative consequences of eutrophication, such as harmful algal blooms, the loss of seagrass beds, extirpation of native and proliferation of exotic and often harmful biological species (Bricker et al., 2008; Kennish et al., 2007; Kennish, Fertig, and Lathrop, 2012; Olsen and Mahoney, 2001).

In order to monitor ecosystem changes, it is important to develop reliable biological indicators specific for coastal lagoons (Birk et al., 2012; Borja et al., 2008; Borja and Dauer, 2008; Borja et al., 2013). Diatoms, siliceous microscopic algae, are among the major primary producers in coastal areas and serve as valuable environmental indicators, commonly used to monitor water quality and infer past environmental conditions (Smol and Stoermer, 2010). While the diversity and ecological properties of the freshwater diatoms are relatively well known, coastal diatom flora is less studied. Most applied studies of coastal diatoms exploit their sensitivity to water salinity (Cholnoky, 1968; Kolbe, 1927), which allows using the composition of fossil diatom assemblages to reconstruct marine transgressions and regressions (Horton and Sawai, 2010; Horton et al., 2007; Snoeijs and Wekström 2010). Such studies do not necessarily require fine-resolution species-level diatom taxonomy, as many diatom genera are specific to either fresh or saline waters. Responses of diatom assemblages to human-induced pollution, on the other hand, require species-level characterization, which is often challenging because of the exceptionally high diversity and poorly understood species boundaries of many coastal diatoms.

A number of studies explored the potential use of sediment diatoms as indicators of human impacts in estuaries, other shallow coastal areas, and salt marshes (see reviews by Admiraal, 1984; Trobajo and Sullivan, 2010; Trobajo et al., 2011). Several field and laboratory experiments demonstrated shifts in diatom species composition, preferential growth of particular species, or decrease in community diversity in response to nutrient or organic enrichment (Admiraal and Peletier, 1979; Peletier, 1996). However, effects of nutrients on biomass, diversity, and species composition of diatoms were not always consistent (e.g., Sullivan, 1978; Sullivan and Daiber, 1975). Differential responses of various diatoms to sediment toxicity were also confirmed in laboratory settings (Moreno-Garrido et al., 2003, 2007). To develop diatom indicators of various types of pollution suitable for particular geographic areas, it is important, however, to conduct observational studies across human-impact gradients in natural ecosystems located in these areas. Datasets compiled as a result of such observational studies can be used to model species distributions in response to various environmental variables and to construct models for inferring concentrations of pollutants from sediment diatom assemblage composition. Such inference models have been constructed for inferring water-column nutrients in the Gulf of Finland, Baltic Sea (Weckström and Juggins, 2005) and Biscayne Bay, Florida (Wachnicka, Gaiser, and Boyer, 2011).

However, the problem in determining diatom assemblage responses to nutrient or other types of pollution from observational datasets lies in the existence of strong natural environmental gradients driving patterns of species distributions in coastal environments. In addition to salinity, which is the most obvious factor explaining patterns of species composition in coastal diatom assemblages (e.g., Juggins, 1992; Saunders et al., 2007; Sherrod, 1999), there are other habitat characteristics that influence species distributions including depth or marsh elevation, sediment texture, type of substrate, exposure to waves, tidal range, irradiance, and temperature (Aleem, 1950; Amspoker and McIntire, 1978; Nelson and Kashima, 1993; Ribeiro et al., 2013; Rovira et al., 2012; Sherrod, 1999; Ulanova, Busse, and Snoeijs, 2009; Ulanova and Snoeijs, 2006; Weilhoefer and Nelson, 2015; Whiting and McIntire, 1985; Zong and Horton, 1998). Natural and human-impact gradients often coincide in coastal environments. For example, pollution sources are often located in the upstream portions of the estuaries, so that it is difficult to disentangle effects of pollution and salinity on diatoms (e.g., Underwood, Phillips, and Saunders, 1998). A careful examination of diatom responses to various factors is thus required to develop diatom-based environmental indicators for specific regions and types of ecosystems.

Little has been known so far on sediment diatoms of the New Jersey coast. Hein and Koppen (1979) investigated diatoms on artificial substrates installed at the Oyster Creek Nuclear Station, located in Barnegat Bay, and reported the effect of thermal pollution on the abundance of common species. Potapova, Desianti, and Enache (2015) studied the relationships of subtidal surface sediment diatoms and contaminants in New York and New Jersey coastal lagoons and found an increase of small planktonic diatoms on the most polluted sediments, although the contribution of other factors, such as eutrophication, to this shift in diatom assemblages could not be ruled out. Phytoplankton of Barnegat Bay was studied by Mountford (1971) and Olsen and Mahoney (2001), who noted an increase of small-celled species as a result of nutrient enrichment. Such an increase of small planktonic diatoms in sediment diatom assemblages caused by eutrophication has been reported from many coastal areas around the world, including such locations as Chesapeake Bay (Cooper, 1995) and Long Island Sound (Cooper, Gaiser, and Wachnicka, 2010; Varecamp et al., 2010).

This study belongs to the series of projects initiated by the New Jersey Department of Environmental Protection, with a common goal to develop indicators for measuring the ecosystem health of New Jersey's shallow coastal lagoons. The objectives of this research was to determine (1) which factors were the main drivers of the composition of sediment diatom assemblages in New Jersey coastal lagoons and (2) which properties of the diatom assemblages could be used as indicators of environmental conditions, including those induced by human impacts in these water bodies.

Our study was conducted in two coastal lagoons: Barnegat Bay–Little Egg Harbor estuary and Great Bay. BB-LEH is composed of three segments (Baker et al., 2014): north, central, and south (Figure 1). It is a narrow water body extending approximately 70 km along the central New Jersey coastline and is approximately 2 to 6 km wide and up to 7 m deep. The area of BB-LEH is 279 km², while the watershed covers approximately 1,700 km². The Bay is a back barrier island lagoon system with three connections to the ocean (Manasquan, Barnegat, and Little Egg inlets). The land-use types vary throughout the Bay watershed. The predominant land-use type in the northeastern mainland area is residential urban land. which includes major population centers such as Toms River and Lakewood. The southeastern mainland area is less heavily developed than northeastern watershed area and contains several protected wildlife refuges and management areas. However, the amount of developed area has increased over the last 70 years in both southeastern and northeastern portions of the watershed. The complex of barrier islands on the eastern shore of the estuary is highly developed, with the exception of Island Beach State Park. Much of the western portion of the watershed lies in the Pinelands National Reserve, a protected area under the Pinelands Comprehensive Management Plan. It is characterized by forested land and low-density development (Baker et al., 2014; Hunchak-Kariouk and Nicholson, 2001).

The Great Bay represents the estuary of the Mullica River. It is located south from the BB-LEH and is connected to the ocean via the Little Egg Inlet (Figure 1). Average water depth in the Great Bay is 1.5 m. Extensive areas of the bay bottom are covered by macroalgae and seagrasses. The adjacent land comprises a variety of habitats such as intertidal marshes, mudflats, and sandflats. In comparison to the BB-LEH, the Great Bay watershed is considerably less developed so that the Bay is a relatively pristine, unaltered ecosystem. This productive estuary supports a high diversity of aquatic and terrestrial habitats and species, especially marine and estuarine fisheries populations, bird colonies on the salt marsh islands, migrating and wintering waterfowl, rare brackish and freshwater tidal wetland communities, plants, and invertebrates (Dowhan et al., 1997).

One hundred sampling sites were organized in 33 north-to-south transects positioned approximately perpendicular across the Barnegat Bay–Little Egg Harbor and Great Bays coastline, with three sampling sites located along each transect (Figure 1). Samples from the Great Bay were collected to represent reference conditions, in contrast to samples from BB-LEH that has considerably more developed watershed. The distance between adjacent transects was approximately 2 km. Within each transect one sampling site was in the intertidal zone and two sites were in subtidal zone: one near shore and another farther offshore. In the intertidal zone mudflats, tidal channels and salt pans were sampled. One sampling site was located approximately in the center of the Great Bay.

Field sampling was conducted in June 2012 in BB-LEH and in mid-July 2012 in the Great Bay. Sediment samples were collected using a variety of different devices depending on water depth and nature of sediment. Short cores were collected using a ∼8 cm diameter acrylic core barrel at intertidal sites. The barrel was slowly pushed into the sediment to minimize compaction. In subtidal locations, a Glew-modified gravity corer (Glew, 1989) or an Ekman Grab was used depending upon sediment consistency (organic or unconsolidated, respectively). The upper 1 cm layer of the core/Ekman Grab was extruded and placed into a precleaned bottle and stored according to parameter of interest (diatoms or chemistry). Samples for sediment chemistry were kept in the dark and on wet ice before return to the laboratory and were then stored frozen until preparation and analysis. Samples for diatom analysis were stored in the dark in a refrigerator and not frozen. Water samples were collected into precleaned HDPE bottles and stored on wet ice in a cooler until returning to the shore-based laboratory or facility for filtration. In subtidal locations water samples were collected directly above the sediment sampling site, while in intertidal sites they were collected from the adjacent waterway.

Water temperature, salinity, conductivity, and pH were measured in situ with the YSI 556 hand-held meter. Surface water samples were filtered through prerinsed and preweighed Whatman GF/F filters (47 mm diameter, 0.7 μm nominal pore size) for total suspended solids and dissolved nutrients, and through precombusted 25 mm GF/F filters for particulate organic carbon and total nitrogen measurements. Turbidity was determined by nephelometric method using a HACH 2100P turbidimeter (U.S. Environmental Protection Agency Staff, 1997). Total suspended solids (TSS) were determined gravimetrically after drying the residue retained on a glass fiber filter at 103–105°C. Chlorophyll a concentration in the water column was determined on a Turner Design fluorometer after filtration (25 mm GF/F filters) then extraction with acetone:water (90:10). Dissolved ammonia + ammonium–nitrogen (NH₃) was determined by an Alpkem Autoanalyzer (RFA 300), using the colorimetric phenate method (U.S. Environmental Protection Agency Staff, 1997). Total Kjeldahl Nitrogen (TKN) was determined by Alpkem Autoanalyzer, using semiautomated block digester and colorimetric phenate method (U.S. Environmental Protection Agency Staff, 1997). Dissolved nitrate + nitrite–nitrogen (NO₃) was determined by an Alpkem Autoanalyzer, using cadmium reduction of nitrate to nitrite, followed by diazotization (U.S. Environmental Protection Agency Staff, 1997). Total phosphorus (TP) and total dissolved phosphorus (TDP) were determined by persulfate digestion. The resulting orthophosphate concentration was measured on the Alpkem Autoanalyzer by the ascorbic acid colorimetric method (U.S. Environmental Protection Agency Staff, 1997).

Total sediment organic carbon (Cs) and total sediment nitrogen (Ns) were measured using a CE Flash Elemental Analyzer following the guidelines in EPA 440.0, manufacturer instructions and ANSP-PC SOP. Samples were pretreated with acid to remove inorganic carbon. Total sediment phosphorus (Ps) was determined using a dry oxidation method modified from Aspila, Agemian, and Chau (1976) and Ruttenberg (1992). Solubilized inorganic phosphorus was measured with standard phosphate procedures using an Alpkem Rapid Flow Analyzer. Standard reference material (spinach leaves) and procedural blanks were analyzed periodically during this study. All concentrations were reported on a dry weight basis. Sediment contaminants were measured at 50 subtidal sites (Figure 1). Samples for organic contaminant analyses were mixed with sulfuric acid and extracted with dichloromethane using a Soxhlet extractor for 18 hours. Solid–liquid chromatograph using florisil (hexane as the eluent) was performed as an additional clean-up step to remove lipids and other compounds. Polycyclic aromatic hydrocarbons (PAH) were quantified using a capillary gas chromatograph coupled with a mass spectrometer in the electron impact mode after a clean-up procedure employing liquid–solid chromatography with alumina as the stationary phase (Ashley and Baker, 1999). After PAH determination, samples were further cleaned using liquid–solid chromatography with florisil as the stationary phase. Internal standards were added to all the samples and calibration standards prior to instrumental analysis: 2,3,6-trichlorobiphenyl (CB30) and 2,2′,3,4,4′,5,6,6′-octachlorobiphenyl (CB204) for polychlorinated biphenyls (PCB). One hundred and eleven PCB congeners, either singly or coeluting, were analyzed using a Hewlett Packard 6890 gas chromatograph equipped with a ⁶³Ni electron capture detector and a 5% phenylmethyl silicon capillary column. The carrier and make-up gases were hydrogen and argon/methane, respectively. The temperature program began with 100°C for 2 min, 100–170°C at 4°C min⁻¹, 170–280°C at 3°C min⁻¹, and 5 min at 280°C. The injector and detector temperatures were 225 and 285°C, respectively. Two-microliter samples were injected with an auto sampler (HP 7673) in the splitless injection mode. The identification and quantification of PCB congeners followed a previously published method (Swackhamer, 1987) in which the identities and concentrations of each congener in a mixed Aroclor standard (25:18:18 mixture of Aroclors 1232, 1248, and 1262) were determined by calibration with individual PCB congener standards. Reported concentrations were for the sum of 110 congeners (some coeluting expect for congener 1 and 3), 40 aromatic hydrocarbons (not including naphthalene), and 22 pesticides (i.e. DDXs, chlordanes, BHC, and others). Trace metals were digested and analyzed at the Geochemical and Environmental Research Group (GERG), Texas A&M University. Wet samples were homogenized in their container, and an aliquot was freeze-dried and homogenized to a fine powder. Approximately 1.00 gram of powdered sediment was mixed with acids and heated to 95°C. After acid digestion, the digestant was diluted to specific volume with reagent water, mixed, and either filtered or allowed to settle overnight before analysis if required. The solutions were then analyzed for trace metals using appropriate instrumental method (ICP-OES or ICP-MS). Grain size distribution was analyzed by the sieve method and proportions of size fractions of >63 μm and <63 μm were reported as “Sand + Gravel” and “Silt,” respectively.

Spatial analysis of land use for each transect was conducted using existing 30 m resolution maps of land cover (National Land Cover Database [NLCD] 2006), 30 m resolution maps of geology, and recently updated 2 m resolution maps of land cover. Land cover was analyzed at several intervals of distance from each raster grid cell. Watersheds were delineated using 10-m resolution digital maps of elevations (DEM). Medium-resolution maps of streams and other human made channels (National Hydrography Dataset) were burned or embedded into the DEM. Burning streams into the elevation model corrected artificial “water dams,” such as bridges, and enhanced the modeled flow of water to the drainage point. All minor depressions in the elevation model were then filled using the ArcGIS sink procedure. Once the DEM had been preconditioned, each sampling location was inspected and moved to the watershed drainage point. Each watershed was delineated separately using the ArcGIS watershed function and inspected to verify that proper watershed was captured. The tasks of checking and then adjusting the delineation took several iterations.

Two hundred training points were used to classify 2010, high-resolution National Agriculture Imagery Program (NAIP) imagery (http://www.fsa.usda.gov/FSA/apfoapp?area=homeandsubject=progandtopic=nai). The NAIP imagery was 1-m resolution and was collected during the growing season. Our training points came from onscreen digitization. Training data captured 14 types of land use. The land-use classes were aggregations of the standard land-use classes used in the government's NLCD land-use maps. They were developed–open space, developed–low intensity, developed–medium intensity, developed–high intensity (aggregated into developed land cover); barren land, deciduous, evergreen and mixed forest (aggregated into forest land cover); shrub/scrub, grassland/herbaceous, pasture/hay (aggregated into grassland land cover); cultivated crops (agriculture land cover); woody wetlands, and emergent herbaceous wetlands (aggregated into wetlands land cover) (Fry et al., 2011). Land use was then summarized for each watershed.

About 1 g of sediment (wet weight) was boiled with 70% nitric acid to oxidize the organic component and was then repeatedly rinsed with deionized water until it reached a neutral pH. The cleaned material was then dripped onto a microscope cover slip and dried. Cover slips were then mounted onto slides using a high refractive index mounting medium (Naphrax™). Diatoms were counted and identified using a Nikon Eclipse 80i microscope equipped with DIC optics. Five hundred valves were counted for each slide at 1000× magnification. All identifications were made to species/variety level when possible, mainly using diatom floras of Hustedt (1955), Krammer and Lange-Bertalot (1986, 1988, 1991a,b), Simonsen (1987), and Witkowski, Lange-Bertalot, and Metzeltin (2000), but also many other taxonomic and floristic works, including Giffen (1963), Hohn and Hellerman (1966), and Sullivan and Reimer (1975). The genus Chaetoceros was represented mostly by spores and small vegetative cells that could not be identified to species level; therefore, its representatives were lumped into “Chaetoceros spp.” category. Species of the genus Skeletonema were likewise not distinguishable under light microscopy and were lumped into Skeletonema spp.

To characterize environmental variability in the study area, a principal component analysis (PCA) was carried out using the dataset of 24 water-quality and land-use variables measured for 100 sites from Barnegat and Great Bays. Most variables were log or square-root transformed prior the analysis.

The relationships between diatom assemblage composition and the environmental variables were explored using nonmetric multidimensional scaling (NMDS) and canonical correspondence analyses (CCA). NMDS was used to determine the major patterns of variation of the samples based on differences in diatom species composition; environmental parameters were superposed as passive variables on an NMDS diagram to visualize their relationships with diatom assemblages. A series of CCAs were carried out to determine the strength of the relationships between diatom assemblage composition and specific environmental variables, one at a time. Species data were square-root transformed for apparent large values, and Wisconsin double standardization was performed on data with values larger than common abundance class scales. Only diatom taxa with high occurrence (in more than 50% of the total samples) were displayed on the diagram. NMDS was carried out in R using vegan package (Oksanen, 2015) PCA and CCA were carried out with the CANOCO software (ter Braak and Šmilauer, 2002).

Diatom transfer functions were constructed using the weighted-averaging–partial least squares (WA-PLS) method as implemented in the C2 software (Juggins, 2003). Bootstrapping was used to obtain objective error estimates for the models. Model performance was estimated by coefficients of determination (R² and R² boot) and the root-square mean errors of prediction (RMSEP). Species optima were calculated as weighted averages and tolerances as weighted-average standard deviations.

The dataset of 100 sites covered a wide range of land use, sediment texture, salinity, chlorophyll a, and nutrient concentrations in water and sediments (Table 1). The first two axes of the PCA of the environmental characteristics of the sites accounted for 44.6% of the total variance, with axes 1 and 2 explaining 24.8% and 19.8% of the total variance, respectively. Two environmental gradients were apparent in the dataset (Figure 2A). One was related to habitat differences and separated subtidal sites with coarser-grained sediments (left side of the PCA diagram in Figure 2A) from intertidal sites with relatively fine sediments and higher concentrations of organic matter in both sediments and the water column (right side of the PCA diagram). The second gradient largely corresponding to the second PCA axis encompassed north–south variation in salinity, TDP, and ammonia concentrations (Figure 2A). The largest variation in the land use was also related to the same north–south gradient, since the northern part of BB-LEH watershed has the highest proportion of the urban (“developed”) land use.

The first two axes of the PCA of the environmental characteristics of the subtidal sites accounted for 52.2% of the total variance, with axes 1 and 2 explaining 38.4% and 13.8% of the total variance, respectively (Figure 2B). Sediment contaminants, such as PAH, PCB, As, Cu, Ag, and Ni, were positively correlated with fine sediments and higher concentrations of organic matter in both sediment and the water column (right side of the PCA diagram). Other contaminants (not included in the PCA) followed the same trend and were the highest in the northern part of BB-LEH, especially in fine-grained sediments rich in organic matter. The contaminants measured in this study were compared to the published sediment quality guidelines, such as “effects range low” (ERL) and “effects range median” (ERM) (Long et al., 1995). For trace metals, no location had concentrations above the ERM. There were concentrations that exceeded the ERL for arsenic, cadmium, copper, lead, and zinc. For total PCBs, only 4 out of 54 stations had concentrations above the ERL, and none above the ERM, similarly only three stations were above the total PAH-ERL (4.0 μg/g dw). Overall, concentrations were low compared to the published guidelines of Long et al. (1995).

A total of 603 diatom taxa belonging to 111 genera were found. Of these taxa, 271 could not be identified to species level (Figure 3). Taxa richness ranged from 43 to 122, with an average of 78 taxa per sample. Shannon-Wiener diversity index based on the natural logarithm ranged from 2 to 4.3. Only 110 taxa reached at least 1% in five or more samples. The most common (occurrence >75%) taxa were Cyclotella atomus var. gracilis Genkal & Kiss, Navicula salinicola Hustedt, Cyclotella choctawatcheeana Prasad, Thalassiosira proschkinae Makarova, Planothidium delicatulum (Kützing) Round & Bukhtiyarova, Nitzschia frustulum (Kützing) Grunow, and Opephora sp. 2 COAST. Chaetoceros spp. occurred in 54% of samples, mostly in the northern part of Barnegat Bay. The most diverse genera were Navicula (128 taxa), Nitzschia (68 taxa), Fallacia (24 taxa), and Cocconeis (22 taxa).

NMDS analysis of the dataset of 100 subtidal and intertidal sites showed that diatom assemblage composition was strongly related to the north–south gradient of salinity and land use (Figure 4A,B). The habitat was also strongly related to diatom assemblages as indicated by the long arrow for depth and good separation of the subtidal and intertidal sites.

The dominant diatoms in the northern part of Barnegat Bay characterized by low salinity and the high degree of urbanization were Chaetoceros spp., Cyclotella choctawatcheeana, Fallacia cryptolyra (Kützing) Stickle & Mann, Fragilaria amicorum Witkowski & Lange-Bertalot, Nitzschia dissipata (Kützing) Grunow, and Nitzschia palea (Kützing) Smith. The southern part of Barnegat Bay and Great Bay had large populations of Cyclotella atomus var. gracilis, C. striata (Kützing) Grunow, Cocconeis cf. scutellum, C. stauroneiformis (Smith) Okuno, Paralia sulcata (Ehrenberg) Cleve, and Thalassiosira proschkinae. Epipsammic diatoms, including Adlafia sp. 3, Amphicocconeis disculoides (Hustedt) Stefano & Marino, Amphora staurophora Juhlin-Dannfelt, Cocconeis neothumensis var. marina Stefano et al., Opephora sp. 2, and Opephora sp. 8, were most abundant on subtidal coarser-grained sediments, while in the intertidal marshes and mudflats Denticula subtilis Grunow and several species of Halamphora and Navicula were prevalent.

NMDS of the 50 subtidal sites for which sediment contaminant data were available (Figure 5) shows that the main gradient in diatom data was associated with north–south salinity and land-use gradient. Contaminant concentrations in sediments were the highest in the northern part of Barnegat Bay, which also had higher chlorophyll a and particulate phosphorus concentrations in the water column. The high proportion of small planktonic diatoms, such as Chaetoceros spp. and Cyclotella choctawatcheeana, in sediment samples from northern sites coupled with high chlorophyll a and particulate phosphorus in the water column indicate their abundance in plankton.

Strong association of diatom assemblages with north–south salinity gradient was also evident from the NMDS of 34 intertidal sites (Figure 6); although, in this dataset, the north–south gradient also comprised an increase of sediment carbon and nitrogen in the northern portion of BB-LEH. Planktonic diatoms, although clearly less prominent than in the subtidal sites, were nevertheless relatively abundant in intertidal sediments.

CCA testing of the effect of specific environmental variables confirmed that salinity had the strongest influence on the composition of sediment diatom assemblages (Table 2). Other variables that had a high degree of association with diatom assemblage variation were TDP, carbon and nitrogen content of sediment, chlorophyll a, and depth. Diatom species distributions were also associated with sediment contaminants, especially PAH, Cu, and Cd (Table 3). Strong correlations of salinity with chlorophyll a and TDP coupled with known sensitivity of diatoms to salinity indicate, however, that the response was primarily to this latter factor. The effects of carbon and nitrogen in sediments as well as sediment contaminants were impossible to separate, since these variables were highly correlated (Figure 2A,B). F ratios calculated for intertidal and subtidal subsets (Table 2) were smaller than for the 100 sites dataset because of the lower number of observations, but the general order of importance of environmental variables was approximately the same.

Inference models were developed for salinity and sediment nitrogen content (Table 3; Figure 7). The models with best performance for both variables were WA-PLS 2nd component. Predictive power of the models was quite high, especially for the salinity model, but both models tended to overestimate low values and underestimate higher values of the inferred variables (Figure 7). Optima and tolerances of dominant diatom species by weighted averaging (Table 4) show that species that are the best indicators of nigh nitrogen content in the sediment are Fragilaria amicorum, Chaetoceros spp. and Cyclotella choctawatcheeana, in low-salinity areas and Navicula salinicola, Nitzschia frustulum, and Planothidium delicatulum in higher salinity situations.

This study revealed highly diverse diatom assemblages in the subtidal and intertidal habitats of New Jersey coastal lagoons. Species richness was higher than reported in other studies of the comparable habitats (e.g., Benito, Trobajo, and Ibáñez, 2015; Sullivan, 1975; Wachnicka et al., 2010). One reason for such a high diversity may be that the studied lagoons provided relatively well-sheltered habitat protected from wave action and strong currents. Another reason is that some previous studies of marsh diatoms (Sullivan 1975) used different sample collection methods that excluded nonmotile diatoms from the analysis.

The analyses showed that diatom assemblages in BB-LEH and Great Bay were mostly affected by salinity, the types of habitat (subtidal vs. intertidal), sediment texture, and organic matter content expressed either in terms of sediment carbon or nitrogen. Our findings are in agreement with observations made in other coastal ecosystems where considerable variation in these environmental characteristics exists (e.g., Amspoker and McIntire, 1978; Benito, Trobajo, and Ibáñez, 2015; Ribeiro et al., 2013). The information on species composition of diatom assemblages found in specific habitats and at different salinity levels is valuable for paleoreconstructions of the environmental conditions in the study area.

The confounding effects of salinity and habitat characteristics, however, create difficulties for the development of indicators of human-induced pollution. For instance, in this dataset the complex north-to-south ecological gradient combined gradients of salinity, TDP, ammonia, chlorophyll a, sediment contaminants, and land use. TDP increased toward the south, while higher nutrient loadings are characteristic for the northern part of BB-LEH (Baker et al., 2014). Elevated TDP may not be related to cultural eutrophication but may instead be caused by the release of phosphates from terrigenous sediments with increased salinity, a phenomenon commonly observed in estuaries (Jordan et al., 2008). The interpretation of species responses is even more obscured by the collinearity between salinity and variables clearly indicating human impacts, such as percentage of developed land use, contaminants, and chlorophyll a. Relatively high water-column chlorophyll a concentrations in the northern part of BB-LEH indicate increased primary productivity due to nutrient enrichment. The elevated proportions of small planktonic diatoms, such as Chaetoceros spp. and Cyclotella choctawatcheeana in the sediment assemblages of this area also support this hypothesis. Blooms of small planktonic diatoms and other small-celled phytoplankton have been commonly associated with increased nutrient availability (Jaanus et al., 2009). Many studies of sediment cores report a considerable increase of small centric diatom species due to eutrophication of coastal areas (Andrén, 1999; Andrén, Andrén, and Kunzendorf, 2000; Cooper, 1995). Chaetoceros spp. is often associated with nutrient enrichment due to upwellings (Nave, Freitas, and Abrantes, 2001), while C. choctawatcheeana has been recognized as a signature species of coastal “anthropogenic assemblages” (Leśniewska and Witak, 2011). Presence of allochthonous diatoms, such as planktonic species in sediment diatom assemblages, is considered an impediment for the investigations of the relationships between benthic diatoms and environmental factors and for reconstruction of past environmental conditions from sediment diatom assemblages (Vos and de Wolf, 1993), and it is often recommended to exclude these species from the analyses (Simonsen, 1969). However, the proportions of the small-celled planktonic diatoms in these sediment samples were so high that it made no sense to discount them in the hope of elucidating responses of benthic species. Benthic diatoms may decline in areas experiencing phytoplankton blooms because of the light limitation, so that the high relative abundance of planktonic species is by itself a good indicator of a regime shift caused by eutrophication.

Although the data did not reveal direct links between water nutrient concentrations and benthic diatom assemblage composition, there was a relationship of diatoms to sediment organic content. Such a relationship has been also found by Amspoker and McIntire (1978) in the Yaquina Estuary, Oregon. These findings are not surprising considering that diatoms are sensitive indicators of the dissolved organic matter, and this sensitivity is the basis of the “saprobic” system widely used to monitor water quality in freshwaters (Kolkwitz and Marsson, 1902; Sladecek, 1986). The increase of sediment organic content, however, may be attributed to a number of factors, not necessarily linked to nutrient pollution. It is definitely relatively high in vegetated marshes in comparison to mudflats and sandflats, and especially subtidal sediments, and may depend on the degree of exposure to wave action and winds and the intensity of decomposition. At the same time, there is evidence that nutrient pollution leads to increased carbon and nitrogen in coastal sediments (Frascari, Matteucci, and Giordano, 2002; Savage, Leavitt, and Elmgren, 2004). Since nitrogen and carbon content of sediments in the northern, most impacted part of BB-LEH was elevated in comparison to other areas in the study, the conclusion is that diatom responses to the organic matter content may be exploited for the purposes of bio-indication. Organic matter content is usually higher in fine-grained sediments, while in this study area, intertidal sediments in the northern part of BB-LEH were relatively coarse. This further suggests that relatively high organic matter content may have resulted from eutrophication. Since nitrogen and carbon in sediments were highly correlated, their effects could not be separated and interpreted at this point. The authors hope that future studies will disentangle these relationships and refine indicative properties of coastal diatoms.

Inference models developed in this study have relatively high predictive power, especially the salinity model, which is often the case in the coastal settings. While diatom response to salinity is unquestionable, inference models developed for water-column nutrients are less straightforward (Juggins, 2003; Wachnicka et al., 2010). Responses to nutrients are often confounded by the presence of other strong gradients in the observational datasets. A strong response to nutrients may be observed when salinity variation in the dataset is limited, as in the study conducted by Weckström and Juggins (2005) in the Gulf of Finland. In this study, the diatom response to water-column nutrients was weak in comparison to other environmental variables. This is probably the consequence of several factors. As Figure 2 indicates, there was no pronounced increase in the water-column nutrients along the gradient of human impact expressed by developed land use. Dissolved nutrients may be immediately taken up by algae and thus not present in the available forms in the water column. The relatively high concentration of chlorophyll a in the most impacted northern part of BB-LEH confirms this possibility. The model developed for sediment nitrogen content may be nevertheless useful for both past environmental reconstructions and for monitoring eutrophication. Further refinements of both inference models are possible. For example, additional sampling from various habitats may lead to the development of separate, more accurate models for subtidal and intertidal habitats.

This is the first detailed study of sediment diatoms in New Jersey coastal lagoons. A rich diatom flora comprised of 603 taxa was found in studied samples, with 271 taxa not yet reported in the literature. A strong north-to-south salinity gradient in the study area was the main factor explaining variation in the composition of diatom assemblages. Subtidal sediments, especially in highly productive areas in the northern part of the Barnegat Bay, contained a high proportion of planktonic diatoms, especially small-celled Cyclotella and Chaetoceros species. Habitat type and physical properties, such as particle size and depth, were other important factors structuring diatom assemblages. Water-column nutrients and sediment contaminants did not show much effect on sediment diatoms either because of the overriding effect of salinity or low variability in nutrient/contaminant concentrations. The organic matter content of sediments, however, was significantly related to diatom species composition. Both sediment carbon and nitrogen were considerably higher in the northern part of the Barnegat Bay and other areas experiencing strong human impacts. This study developed diatom inference models for salinity and sediment nitrogen content and concluded that diatom species with relatively high optima for sediment nitrogen may be used as indicators of nutrient enrichment in studied lagoons.

This research was funded by the New Jersey Department of Environmental Protection, Division of Science, Research and Environmental Health through an agreement with the New Jersey Sea Grant Consortium (2011-30). We would like to thank several colleagues who carried out field sampling and laboratory analyses. William Wallon helped collecting sediment and water samples in the field. Paul Kiry, Linda Zaoudeh, and Paula Zelanko carried out sediment and water chemistry analyses. Sylvan Klein conducted sediment grain size analysis. We are grateful to Nick Procopio of NJDEP and three anonymous reviewers for the constructive comments that helped to improve this paper.