Spatial and Temporal Distributions of Live Salt-Marsh Foraminifera in Southern New Jersey: Implications for Sea-Level Studies

Salt-marsh foraminifera have played a crucial role in reconstructing past relative sea-level (RSL) change, which began with the recognition of the tightly constrained vertical zonation of modern marsh foraminiferal assemblages (Scott, 1976; Scott & Medioli, 1978, 1980a). Investigations of modern foraminifera from salt marshes along North American coastlines validated the utility of foraminifera as reliable sea-level indicators due to their vertical distributions’ robust relationship to tidal elevation (e.g., Thomas & Varekamp, 1991; Varekamp et al., 1992; Gehrels, 1994, 1999; Nydick et al., 1995; Kemp et al., 2009). Salt-marsh foraminifera are ideal proxies for quantitative analysis to reconstruct RSL due to their low diversity, high abundance assemblages that are preserved in coastal sediment (e.g., Gehrels, 2007). The development of foraminiferal-based transfer functions quantified distributions of modern assemblages with tidal elevation to apply to fossil counterparts enumerated from sediment cores (e.g., Gehrels, 1999; Horton et al., 1999). The transfer function approach has been used to produce estimates of paleomarsh elevation and, when combined with a core chronology, continuous records of RSL at decimeter vertical resolution (e.g., Gehrels, 2000; Horton & Edwards, 2006; Kemp & Telford, 2015; Cahill et al., 2016).

In order to produce reliable estimates of RSL derived from foraminiferal-based transfer functions, a detailed understanding of salt-marsh foraminifera distributions in the modern environment over different spatial and temporal scales is required (e.g., Horton, 1999). While foraminifera have often only been analyzed through one sampling period without replicates, modern foraminiferal assemblages have been found to vary both spatially (sub-meter scales) (e.g., Buzas, 1968; Swallow, 2000; Morvan et al., 2006; Kemp et al., 2011) and temporally (seasonally and interannually; e.g., Buzas et al., 2002; Hippensteel et al., 2002; Martin et al., 2002; Horton & Edwards, 2003; Horton & Murray, 2006; Berkeley et al., 2008). In particular, since early studies of foraminifera, it was suggested that detailed observations of living (at time of collection) assemblages over months to years are necessary to fully understand the foraminiferal population (Buzas, 1968). These observations can then be used to determine how the live assemblage is incorporated into the dead assemblage and affects the total (live + dead) modern foraminifera assemblage (e.g., Horton, 1999), and to determine which assemblage (live, dead, or total) is best suited in a foraminiferal-based transfer function to reconstruct RSL (e.g., Scott & Medioli, 1980b; Murray, 2000; Murray & Bowser, 2000; Culver & Horton, 2005). Furthermore, the composition and elevational range of foraminiferal assemblages can vary among sites and regions due to factors such as climate and salinity (e.g., de Rijk, 1995; Hayward et al., 1996; Wright et al., 2011), so it is necessary to have a full understanding of the environmental variables influencing foraminiferal distributions in order to use them as a sea-level indicator.

We conducted a three-year monitoring study of live salt-marsh foraminifera from four high marsh stations along a salinity gradient in southern New Jersey. Previous investigations at the sites established the relationship of dead assemblages (Walker et al., 2020), but the variability of living assemblages and their impact on the total assemblage is unknown. We collected samples at each station every three months to assess seasonal and interannual changes and to evaluate small-scale spatial variability. In addition, we measured porewater and sedimentary variables and stable carbon isotopes from each monitoring station. Stable carbon isotopes were specifically analyzed to evaluate seasonal and interannual variability because of their application as a secondary proxy in foraminiferal-based transfer functions from salt-marsh environments (Middelburg et al., 1997; Johnson et al., 2007; Kemp et al., 2012b). This study is one of the longest seasonal/interannual monitoring assessments of salt-marsh foraminifera compared to other studies ranging from several months to several years (e.g., Scott & Medioli, 1980b; Hippensteel et al., 2002; Horton & Edwards, 2003; Horton & Murray, 2006). We found that elevation was the dominant controlling factor influencing live, dead, and total foraminifera distributions, with secondary influences from salinity and substrate. There was no clear correlation between variability in live assemblages and measured environmental variables over time. However, we identified four distinct clusters of live foraminifera assemblages that correspond to the stations from which they were sampled and also correspond to the dead and total assemblages from the same stations (Walker et al., 2020). Our study continues to highlight the applicability of salt-marsh foraminifera to reconstruct RSL nearly 50 years after the inferences of Scott (1976) and Scott & Medioli (1978; 1980a).

The salt-marsh field sites are located in southern New Jersey on the U.S. mid-Atlantic coast (Fig. 1). Barrier island systems line the coast where inlets allow water to exchange between the Atlantic Ocean and inland bays (Ferland, 1990). Along the coast, extensive modern salt marshes form gently sloping platforms (Ferland, 1990). Our field sites are in the Mullica River-Great Bay estuary near the Rutgers University Marine Field Station in Tuckerton, New Jersey. The 1474-km² watershed provides an excellent monitoring location because the coastlines encompass extensive pristine Spartina-dominated marshes relatively unaltered by human activities (Good & Good, 1984; Lathrop et al., 2000; Kennish, 2004). The estuary is part of the Jacques Cousteau National Estuarine Research Reserve, has low human population density with a lack of development, and a natural inlet (Little Egg Inlet) from Great Bay to the ocean (Fig. 1B; Chant et al., 2000; Kennish, 2004). The Mullica River-Great Bay estuary has semidiurnal tides with microtidal (<2 m) ranges that vary from 0.7 m (in Little Egg Harbor) to 1.1 m (near the mouth of Great Bay) estimated by VDatum (Yang et al., 2008).

We obtained climate data from the Jacques Cousteau National Estuarine Research Reserve meteorological station at Nacote Creek, about 12 km from the monitoring stations (Fig. 1). These data are available from the National Oceanographic and Atmospheric Administration (NOAA) National Estuarine Research Reserve System's Centralized Data Management Office. Over the three-year sampling timeframe from August 2014 to June 2017, average monthly air temperatures ranged from –5 to 25ºC (Fig. 2 A). Air temperature minima and maxima occurred each year in January or February, and July or August, respectively. Additionally, six separate months during the three-year study period were recorded as the statewide warmest in New Jersey from 1895–2018 (Office of the New Jersey State Climatologist). These record warm months were May, November, and December 2015; August 2016; and February and April 2017. The meteorological station at Nacote Creek recorded total monthly precipitation ranging from <20 mm to >200 mm (Fig. 2A). Our monitoring included two months with extreme precipitation for New Jersey's 1895–2018 period of record. In spring 2015, the third driest May was immediately followed by the fourth wettest June on record (Office of the New Jersey State Climatologist). The Atlantic City tide gauge, located ∼19 km from Station 1 and ∼27 km from Stations 2, 3, and 4, had monthly mean sea levels with annual lows in February and March and annual highs in September and October (Permanent Service for Mean Sea Level; Holgate et al., 2013; Fig. 2B). The intra-annual variability in mean sea level heights is driven by natural seasonal cycles, as well as by fluctuations in salinities, winds, and currents (e.g., Barbosa et al., 2008). Further, a significant winter storm flooding event occurred in January 2016, which was in the top ten highest historic crests for the Atlantic City tide gauge, according to NOAA.

The sampling methodology follows that described in Walker et al. (2020). We established a 1 m × 1 m plot at four salt-marsh monitoring stations from high marsh/high marsh-upland transition along a salinity gradient (Fig. 1), which we sampled from September 2014 to June 2017. Station 4 was established in March 2015. We chose to study high marsh sites (predominantly used for sea-level studies) because they provide more precise elevation and chronological estimates (e.g., Gehrels, 2000; Kemp et al., 2011). We sampled every three months (September = Fall; December = Winter; March = Spring; June = Summer) over three years to examine temporal variability of salt-marsh foraminifera. We also collected four replicate surface sediment samples during each sampling period to assess small-scale spatial variability. Samples were of a standardized volume of 10 cm³ (10 cm² by 1 cm thick) to allow comparison with similar studies (e.g., Scott & Medioli, 1980b; Horton & Edwards, 2006; Kemp et al., 2012a). Additionally, we sampled a different quadrant of each plot every three months, following Horton et al. (2017), to allow the marsh surface time to recover before it would be sampled again the following year.

We surveyed elevations of each station to NOAA tidal benchmarks using a total station and referenced elevations to the North American Vertical Datum (NAVD88). Because Station 1 is located <500 m from the USGS tide gauge at Little Egg Inlet (station number 01409335), we converted the elevation at Station 1 to tidal datum levels using this tide gauge, as well as VDatum (Yang et al., 2008). For Stations 2, 3, and 4, we deployed two automatic water-level loggers (Solinst Levelogger Edge) in tidal channels within 100 m of each of the stations and leveled them to NOAA tidal benchmarks. We then used the NOAA tide gauge at Tuckerton Creek (station number 8534080) located approximately 1 km north of Stations 2, 3, and 4 to correlate water levels with the logger data. Finally, we converted the tidal elevations from each station into a standardized water level index (SWLI), following Horton et al. (1999) because of differences in tidal range between Station 1 and Stations 2, 3, and 4.

At each monitoring station, we counted live and dead foraminifera from four replicate samples from each sampling period over all three years. The four replicates were separate 10-cm³ surface samples taken within a ∼50-cm² area within a quadrant of each plot. The samples were stained with rose Bengal immediately after collection in order to identify live foraminiferal tests (Walton, 1952) and then stored in a buffered ethanol solution and refrigerated (Scott et al., 2001). We wet sieved each sample to obtain the 63–500 µm size fraction and then split each sample using a wet splitter (Scott & Hermelin, 1993). Identifications of foraminifera (Appendix A) were made under a binocular microscope with samples immersed in distilled water. We did not identify specimens of the genera Haplophragmoides and Ammobaculites to the species level due to difficulties in identifications of these species (Kemp et al., 2009). However, all identifications of foraminifera were confirmed with type specimens at the National Museum of Natural History, Smithsonian Institute, Washington, D.C. While this study focuses on the live assemblages, further results concerning the dead assemblages are available in Walker et al. (2020). All live foraminifera data and coefficients of variation (CV; the ratio of the standard deviation to the mean) to measure the variability among replicate samples are available in Appendix B.

To analyze the composition of the live, dead, and total foraminifera assemblages from the four monitoring stations, we used partitioning around medoids (PAM; Kaufman & Rousseeuw, 1990). We used the entire dataset of relative abundances of live, dead, or total foraminifera for three separate analyses, using four groups to correspond to the four monitoring stations. We completed the analysis using the ‘cluster’ package in R and graphically represented the data with a silhouette plot (Rousseeuw, 1987). Silhouette widths between –1 and 1 provide an estimate of how well a sample fits within its assigned group. Sample values that are close to 1 indicate that the sample was assigned to an appropriate group where dissimilarity within the group was less than the dissimilarity among the four groups. Conversely, sample values that are close to –1 indicate that the sample was not appropriately classified.

We analyzed the stable carbon isotope (δ¹³C) values of surface bulk sediment collected from each monitoring station from each sampling period for all three years, as well as from 13 randomly selected replicate samples throughout the sampling period (includes both sub-sampling of one 10-cm³ surface sample and separate samples taken within a ∼50-cm² area). Due to possible contamination of samples from Station 3 in spring and summer of Year 3, we took two additional samples from Station 3 three months and six months after the summer of Year 3. We analyzed bulk sediment δ¹³C values because the dominant carbon input to salt-marsh sediment is in situ marsh vegetation. Bulk sediment δ¹³C was measured using cavity ring-down laser spectroscopy (CRDS) following the flash combustion technique described by Balslev-Clausen et al. (2013) at the Departments of Geology and Environmental Studies at Bryn Mawr College. Bulk peat samples were freeze-dried in a Virtis™ benchtop freeze dryer prior to isotopic analysis to remove moisture. The samples were then ground in a Retsch™ ball mill until finely powdered. Approximately 2 mg (±0.5 mg) of the prepared sample was weighed on a Mettler Toledo™ XP56 microbalance with 4 µg precision. The samples were sealed in pressed-wall tin capsules and flash combusted at 980°C in a Costech™ ECS 4010 element analyzer, using N₂ as a carrier gas. The isotopic composition of the CO₂ produced by combustion is analyzed in a Picarro™ G2201-i CRDS instrument. For each sample, the carbon abundance was calculated from the peak ¹²CO₂ concentration measured by the CRDS system and was calibrated using standard reference material (NIST 1547 peach leaf). Reproducibility of carbon mass concentration is ±0.8% (1 s.d., n = 96). The carbon isotopic composition is standardized to Vienna Pee Dee Belemnite (VPDB) using standard reference material USGS40 (glutamic acid). The reproducibility of δ¹³C values is ±0.2‰ based on repeat analyses of NIST 1547 (±0.20‰, 1 s.d., n = 530) and USGS40 (±0.09‰, 1 s.d., n = 128).

We measured several environmental variables at each of the four monitoring stations. Using a handheld YSI meter, we measured porewater salinity, temperature, and pH at each station plot during each sampling period. We analyzed organic matter by loss on ignition (LOI) for samples from each station from each sampling period. We dried the samples in an oven and ignited the samples in a muffle furnace following the methods of Plater et al. (2015). Additionally, we analyzed the grain size distribution of samples from each station. We used 30% H₂O₂ to digest the organic fraction of the samples to prepare them for grain-size analysis (Donato et al., 2009). Grain size distributions were measured using a Malvern Mastersizer 3000 laser particle-size analyzer and were described after Folk & Ward (1957).

To analyze the relationship between foraminiferal assemblages [live, dead, and total (live + dead)] and all measured environmental variables (elevation, porewater salinity/temperature/pH, LOI, grain size) and to determine the type of response displayed by the species distribution to environmental gradients (unimodal or linear), we used detrended canonical correspondence analysis (DCCA) ordination methods using the program Canoco, version 5 (ter Braak & Smilauer, 2012). For each assemblage, we used principal component analysis (PCA) because the gradient length measured in standard deviation units was <2, which suggests a linear species-environment relationship (Birks, 1995). Each arrow on the PCA biplot represents a measured environmental variable and indicates the strength and direction of influence of the variable on the distribution of species. The longer the arrow, the more strongly correlated that variable is to species composition. We then used redundancy analysis (RDA) to quantify the total variance in each assemblage that can be explained by the measured environmental variables, and we determined the contribution of each individual variable using interactive-forward-selection with variables with a p-value <0.01 considered as major controls on the assemblage.

In the 184 samples taken over the three-year study, which included 12 separate sampling periods with spatial replicate samples, we counted 10,260 live foraminiferal tests, compared to 72,804 dead foraminiferal tests, for a total of 83,064 tests (Fig. 3). The average standing crop of live foraminifera for one sample was 56 ± 60 tests/cm³ (1σ) and ranged from 0 to 282 tests/cm³ compared to an average of 389 ± 221 tests/10 cm³ for the dead assemblages. Of the 184 samples, 18 samples did not contain any live foraminifera. The greatest total number of live (4978 tests) and dead (26,194 tests) foraminifera tests across all four stations was observed in Year 3 (Fig. 3). The standing crop of live, dead, or total (live + dead) foraminifera at each station did not exhibit a clear seasonal pattern over the three years. Across all four stations, on average, 12% of the total assemblage was live during each sampling period, but this ranged from an average of 4% at Station 3 to 25% at Station 4 and reached a maximum of 46% in fall of Year 2 at Station 4.

The live foraminifera consisted of 11 agglutinated species. Three additional species were found in the dead assemblage, but these species only made up 390 tests of all dead tests (<1%) counted. The dominant species in the live assemblage across all four stations (relative abundance of at least 1% across all samples) in order of most to least abundant were Jadammina macrescens (6,220 tests), Trochammina inflata (1,574 tests), Balticammina pseudomacrescens (962 tests), Ammoastuta inepta (588 tests), Tiphotrocha comprimata (548 tests), and Haplophragmoides spp. (222 tests). These were also the six most abundant species in the dead assemblage.

The PAM analysis of the entire live foraminifera dataset, dead foraminifera dataset, and total foraminifera dataset each had an average silhouette width of 0.56, meaning the samples fit well into four groups because the value is close to 1 (Fig. 4). The average silhouette width for the group containing samples from Station 1 for each analysis ranged from 0.61–0.62, compared to 0.50–0.52 for Station 2, 0.50–0.51 for Station 3, and 0.62 for Station 4. The similarity among the average silhouette widths for each group across the three analyses suggests consistency among the live, dead, and total foraminifera assemblages. The greater the silhouette width the more similar the samples are to each other in that group. Therefore, the foraminiferal assemblages from Stations 1 and 4 are more consistent over time and space compared to Stations 2 and 3. Of the 165 samples containing live foraminifera, only 13 samples were not assigned to the group corresponding to the station from which they were sampled. Twelve of the 184 samples of dead foraminifera and 8 of the 184 samples of total assemblages were not assigned to the correct group; these particular eight samples were not assigned correctly in the live or dead assemblage analysis either. All but two of the incorrectly assigned samples were from Stations 2 and 3, which had the groups with lower silhouette widths, further suggesting that the foraminifera assemblages from these two stations were less consistent over time and space compared to Stations 1 and 4.

Station 1 is primarily vegetated by Spartina alterniflora (short form), has an elevation of 212–223 SWLI units (0.68 ± 0.03 m MTL, 1σ), and is the highest salinity site. We identified 10 total foraminifera species from Station 1. Five species, which were also the most abundant, were found in both the live and dead assemblages, of which four had a relative abundance of at least 1% across all live assemblage samples (Fig. 5 A). The dominant species at Station 1 was T. inflata, which made up 63–66% of the live, dead, and total assemblages. The relative abundances of the other most abundant species were similar among the live and dead assemblages, except for T. comprimata, which had an average relative abundance of 10% in the dead assemblage and only 2% in the live assemblage.

The live standing crop ranged from 0 to 214 tests/10 cm³ with an average of 40 ± 47 tests/10 cm³ (1σ; Fig. 5B). Interannually, Year 3 had the greatest number of live foraminifera with 1334 total tests (average 78 ± 59 tests/10 cm³) compared to Year 1 with 178 tests (average 11 ± 16 tests/10 cm³) and Year 2 with 462 tests (average 29 ± 24 tests/10 cm³). Station 1 had the highest average CV among replicates at 70 ± 45%. The three most dominant live foraminifera species were T. inflata, J. macrescens, and B. pseudomacrescens (Fig. 5). Trochammina inflata had an average standing crop of 27 ± 36 tests/10 cm³ and J. macrescens had an average standing crop of 12 ± 15 tests/10 cm³. Balticammina pseudomacrescens was found in very low numbers with only 68 total tests over the three years. While T. inflata and J. macrescens both had standing crops that were greatest in Year 3 (948 and 368 total tests, respectively), B. pseudomacrescens had 0 live tests in Year 3. In fact, B. pseudomacrescens only had live tests present in 13 of the 48 samples, which occurred in Year 1 and the first half of Year 2. None of the species in the live assemblage exhibited a clear seasonal pattern.

Station 2 is primarily vegetated by Spartina patens, has an elevation of 211–224 SWLI units (0.68 ± 0.03 m MTL), and had the largest live standing crop over the three years (3658 total tests). We identified 11 total foraminifera species from Station 2, of which nine species were found in the live assemblages. Four of the nine species in the live assemblages had a relative abundance of at least 1% across all samples (Fig. 6 A). The dominant species at Station 2 was J. macrescens, which made up 72%, 48%, and 51% of the live, dead, and total assemblages, respectively. Balticammina pseudomacrescens and T. comprimata made up the majority of the rest of the live, dead, and total assemblages, accounting for 23% (live) to 46% (dead) of the assemblages.

The live standing crop ranged from 0 to 282 tests/10 cm³ with an average of 76 ± 71 tests/10 cm³ (1σ; Fig. 6). Interannually, Year 3 had the greatest number of live foraminifera with 1680 total tests (average 105 ± 82 tests/10 cm³) compared to Year 1 with 482 tests (average 30 ± 32 tests/10 cm³) and Year 2 with 1496 tests (average 94 ± 67 tests/10 cm³). The average CV among replicates was 55 ± 33%. The three dominant live foraminifera species were J. macrescens, B. pseudomacrescens, and T. comprimata. Jadammina macrescens accounted for 72% of the total live foraminifera tests counted and had an average of 55 ± 59 tests/10 cm³. While J. macrescens and T. comprimata both had standing crops that were greatest in Year 3 (1196 and 246 total tests, respectively), B. pseudomacrescens had standing crops that were greatest in Year 2 (352 total tests). None of the species in the live assemblage exhibited a clear seasonal pattern.

Station 3 is primarily vegetated by Spartina patens and Distichlis spicata but borders a Phragmites australis flora. It has an uneven surface topography and an elevation ranging from 192–227 SWLI units (0.46–0.60 m MTL). We identified 12 total foraminifera species from Station 3, of which nine species were found in the live assemblages. Eight of the nine species in the live assemblages had a relative abundance of at least 1% across all samples (Fig. 7 A), resulting in the most diverse assemblage of the four stations. There was no clear dominant species in the live or dead assemblages. The highest relative abundance in the live assemblage was B. pseudomacrescens at 22%, compared to T. comprimata at 34% in the dead assemblage.

Station 3 had the smallest live standing crop over the three years (1014 total tests). The live standing crop ranged from 0 to 78 tests/10 cm³ with an average of 21 ± 21 tests/10 cm³ (1σ; Fig. 7). Interannually, Year 3 had the greatest number of live foraminifera with 502 total tests (average 31 ± 23 tests/10 cm³) compared to Year 1 with 80 tests (average 5 ± 8 tests/10 cm³) and Year 2 with 432 tests (average 27 ± 20 tests/10 cm³). The average CV among replicates was 64 ± 48%. The three most dominant live foraminifera species were B. pseudomacrescens (average count size of 5 ± 6 tests/10 cm³), T. comprimata (average count size of 4 ± 6 tests/10 cm³), and T. inflata (average count size of 4 ± 7 tests/10 cm³). While T. comprimata and T. inflata both had standing crops that were greatest in Year 3 (102 and 116 total tests, respectively), B. pseudomacrescens had standing crops that were greatest in Year 2 (128 total tests). Interestingly, T. inflata had ≤10 tests in all of the samples, except for one sample in the last sampling period, Summer of Year 3, where there were 48 live tests. While B. pseudomacrescens and T. inflata did not exhibit a clear seasonal pattern, T. comprimata had maximum abundances in Summer of each year.

Station 4 is a high marsh-upland transition site vegetated by Phragmites australis, has an elevation of 248–264 SWLI units (0.61 ± 0.03 m MTL), and is the lowest salinity site. Station 4 was established in March 2015; therefore, Year 1 refers to data only from spring and summer 2015. We identified 10 total foraminifera species from Station 4, of which eight species were found in the live assemblages. Four of the eight species in the live assemblages had a relative abundance of at least 1% across all samples (Fig. 8 A). The dominant species at Station 4 was J. macrescens, which made up 82%, 66%, and 70% of the live, dead, and total assemblages, respectively. The relative abundances of the other most abundant species were similar among the live and dead assemblages.

The live standing crop ranged from 0 to 268 tests/10 cm³ with an average of 90 ± 65 tests/10 cm³ (1σ; Fig. 8). Interannually, Year 2 had the greatest number of live foraminifera with 1858 total tests (average 116 ± 65 tests/10 cm³) compared to Year 1 with 294 tests (average 37 ± 40 tests/10 cm³) and Year 3 with 1462 tests (average 91 ± 62 tests/10 cm³). Station 4 had the lowest average CV among replicates at 35 ± 26%. The three most dominant live foraminifera species were J. macrescens, A. inepta, and B. pseudomacrescens. Jadammina macrescens accounted for 82% of the total live foraminifera tests counted and had an average of 74 ± 60 tests/10 cm³. Jadammina macrescens and B. pseudomacrescens both had standing crops that were greatest in Year 2 (1602 and 60 total tests, respectively), and A. inepta had standing crops that were greatest in Year 3 (240 and 58 total tests, respectively). Balticammina pseudomacrescens had maximum abundances in Summer of each year. For the other species, the later establishment of Station 4 and the curiously low total number of live tests found in the first sampling period of Spring 2015 compared to the rest of the study timeframe impede the identification of any other clear seasonal patterns in the live assemblages.

The δ¹³C values across the stations ranged from –15.7‰ to –28.0‰ (Fig. 9). Stations 1, 2, and 3 all had average δ¹³C values less depleted than –18.9‰, which is associated with a C₄ dominated salt-marsh plant community and corresponds to all three of these sites’ Spartina vegetation. Station 1 bulk sediment samples had an average δ¹³C value of –16.7‰ ± 0.5 over the three years, Station 2 had an average δ¹³C value of –15.7‰ ± 0.5, and Station 3 had an average δ¹³C value of –17.5‰ ± 1.8. Stations 1, 2, and 3 exhibited some replicate variability (<1.5‰) and variability through time (<2.5‰), but δ¹³C values always remained within the range for C₄ plant communities. Station 4 had an average δ¹³C value of –28.0‰ ± 0.2, indicating a C₃ dominated salt-marsh plant community, which corresponds to Station 4’s Phragmites australis vegetation. Station 4 δ¹³C values exhibited the least variability through time (<0.7‰) among the four stations.

Porewater temperature exhibited seasonal variability that was comparable across the four stations, with highs in the summer and fall (average 22.9 ± 2.3°C) and lows in the winter and spring (average 7.5 ± 2.7°C) each year (Fig. 10 A). Station 1 had the highest porewater salinity with an average of 39.6 ± 8.8 psu (1σ) over the three-year sampling period. Stations 2 and 3 had comparable salinities with an average of 13.7 ± 4.9 psu and 13.1 ± 4.9 psu, respectively (Fig. 10B). Station 4 was the lowest salinity site with an average of 2.9 ± 2.0 psu. The salinity at each station exhibited some seasonal variability, with the highest recorded porewater salinities tending to occur in the fall of each year. The average porewater pH was comparable at each station, ranging from 6.4 ± 0.4 at Station 1 to 6.6 ± 0.5 at Station 4, and did not exhibit any clear seasonal trends (Fig. 10C). Station 2 had the highest average LOI over the sampling period at 81.1 ± 3.3%, compared to 64.3 ± 3.8% at Station 3, and 53.4 ± 7.2% at Station 4 (Fig. 10D). The LOI at Station 1 had an increasing trend over the three years, from 17.4% in Fall of Year 1 to 40.0% in Summer of Year 3. The substrate at Stations 1, 2, and 3 is a sandy silt with relative abundances of sand, silt, and clay of 27–33%, 55–57%, and 10–15%, respectively (Fig. 10E). Station 4’s substrate is a fine silt with relative abundances of sand, silt, and clay of 9%, 66%, and 25%, respectively.

Our PCA with the live foraminifera assemblage produced Axis 1 with an eigenvalue of 0.38 and Axis 2 with an eigenvalue of 0.21 that explains 59% of the total variance in the live foraminifera assemblages with 43% explained by the environmental variables. Constrained and interactive-forward-selection analyses show that elevation is the most significant environmental variable for the live assemblages, explaining 18% of the variance in the dataset, followed by percentage clay (12%), salinity (10%), and LOI (3%; Fig. 11). The PCA results for the dead and total (live + dead) assemblages were very similar. This PCA produced Axis 1 with an eigenvalue of 0.35 or 0.36 and Axis 2 with an eigenvalue of 0.30 or 0.31 that explains 66% of the total variance in the dead and total foraminifera assemblages with 71–72% explained by the environmental variables. Constrained and interactive-forward-selection analyses show that elevation is the most significant environmental variable for dead and total assemblages, explaining 31% of the variance in each dataset, followed by salinity (23–24%), percentage clay (15%), and LOI (2%; Fig. 11).

Each of the four monitoring stations has a unique live foraminifera assemblage, as shown with the PAM analysis. The dominant foraminifera species, T. inflata, J. macrescens, T. comprimata, B. pseudomacrescens, Haplophragmoides spp., and A. inepta, were consistent with the dead assemblages from each station (Walker et al., 2020) and also similar to dominant foraminifera from high marsh sites found in previous studies in New Jersey and on the U.S. Atlantic coast (e.g., Culver et al., 1996; Hippensteel et al., 2000; Kemp et al., 2009, 2011, 2012a). The consistency among the live, dead, and total assemblages is evident as the PAM analysis clusters nearly all of the samples over time and space into groups corresponding to the station from which they were sampled.

While the total count of the live foraminifera assemblages and dominant species fluctuated over the three-year study period at each station, the overall assemblages remained relatively consistent temporally and on sub-meter spatial scales among replicate samples. For example, Year 3 had the greatest standing crops of live foraminifera at Stations 1, 2, and 3, whereas Year 2 had the greatest standing crops at Station 4, but the assemblages making up the standing crops remained the same. Although live assemblages commonly show interannual seasonal variability (e.g., Buzas & Hayek, 2000; Swallow, 2000), clear seasonal trends in the live assemblages at these four stations were absent. The only evidence of seasonal variability was that T. comprimata at Station 3 had maximum abundances in summer of each year, as well as B. pseudomacrescens at Station 4, but these species did not exhibit the same pattern at the other stations. For B. pseudomacrescens, there was no indication of a translation of this trend from the live assemblage to the dead assemblage at Station 4 (Walker et al., 2020), which has been noted in previous studies (Horton & Murray, 2006, 2007; Morvan et al., 2006). In the dead assemblage at Station 3, there were maximum abundances of T. comprimata each year in fall and winter, which could be related to the increase in live T. comprimata found in the summer (Walker et al., 2020). Therefore, while the live foraminiferal assemblages do fluctuate on interannual timescales, there is little evidence for consistent seasonal trends. Further, the live foraminifera assemblages did not exhibit any strong correlations with broader weather changes and seasonal sea-level fluctuations, which included record monthly air temperatures, precipitation extremes, and a significant winter flooding event (Appendix C). The lack of correlation with weather extremes is also consistent with the results of the dead assemblages (Walker et al., 2020).

On small spatial scales, the live assemblages also varied, but the overall assemblages remained relatively consistent among replicate samples. For example, at an extreme, in spring of Year 2 at Station 2, one replicate had only 2 live tests, while another had 212 live tests. However, the dominant foraminifera species remained the same among all replicates at each station. The observed small-scale spatial variability is likely due to factors, such as availability of food resources (Alve & Murray, 2001; Fontanier et al., 2003), response to predation (Buzas, 1978, 1982), reproduction (Stouff et al., 1999), or species interactions (Buzas, 1968; Hayward et al., 1996; Scott et al., 2001).

Although all four monitoring stations are located in high marsh environments, each station has a unique foraminiferal assemblage. High marsh assemblages can vary both among and within regions (e.g., Ellison & Nichols, 1976; Wright et al., 2011; Kemp et al., 2012a); although salt-marsh foraminifera distributions are strongly linked with tidal elevation (e.g., Scott & Medioli, 1978; Horton & Edwards, 2006; Kemp et al., 2012a), secondary environmental factors can also influence foraminiferal assemblages (e.g., de Rijk, 1995; Nikitina et al., 2003; Kemp et al., 2009; Wright et al., 2011; Kemp et al., 2012a). For example, across the four stations, we found through constrained and interactive-forward-selection analyses that elevation was the dominant environmental control on the live (18%), dead (31%), and total (31%) foraminiferal assemblages. However, salinity and the substrate (% clay used here) were also controlling environmental factors influencing the modern foraminifera distributions, and LOI also contributed a small percentage to the variance in each dataset. Porewater temperature and pH were not identified as major environmental controls on foraminifera distributions, although other studies have found temperature and/or pH to have an influence on modern salt-marsh foraminiferal distributions (Horton & Edwards, 2006; Horton & Culver, 2008; Avnaim-Katav et al., 2017).

As the four stations were chosen along a salinity gradient, it is not surprising that salinity was a dominant environmental control on the live (10%), dead (23%), and total (24%) foraminiferal assemblages. At Station 1, the highest salinity site, T. inflata was the dominant species, which has previously been associated with higher salinity sites on the U.S. Atlantic coast (Kemp et al., 2009). By contrast, at Station 4, the lowest salinity site, J. macrescens was the dominant species, which has typically been found in maximum abundances in high marsh-upland transition environments like Station 4 with low salinities (e.g., Spencer, 2000; Nikitina et al., 2003; Robinson & McBride, 2006; Horton & Culver, 2008; Kemp et al., 2009). Further, A. inepta, which was also found in high abundances at Station 4, has been found to be a dominant species in similar low salinity environments on the U.S. Atlantic coast (Scott et al., 2001; Culver & Horton, 2005; Kemp et al., 2009, 2013). Stations 2 and 3 had similar low salinities with several species, such as T. comprimata and Haplophragmoides spp., that have been previously associated with lower salinity environments on the U.S. Atlantic coast (de Rijk, 1995; de Rijk & Troelstra, 1997; Kemp et al., 2013). The fact that Stations 2 and 3 had similar salinities but different foraminiferal assemblages and dominant species also highlights that elevation is still the primary environmental control on the foraminiferal distributions.

The δ¹³C values were analyzed because of their application as a secondary proxy in foraminiferal-based transfer functions from salt-marsh environments (Middelburg et al., 1997; Johnson et al., 2007; Kemp et al., 2012b). Monitoring Stations 1, 2, and 3 had average bulk sediment δ¹³C values ranging from –15.7‰ to –17.5‰ over the three-year sampling period, which is associated with a C₄ dominated salt-marsh plant community and corresponds to all three of these sites’ Spartina vegetation. Kemp et al. (2012b) found similar bulk sediment δ¹³C values for Spartina-dominated salt marsh zones in southern New Jersey, ranging from –18.9‰ to –15.4‰. These values are also consistent with other Spartina marshes on the U.S. mid-Atlantic coast (e.g., Ember et al., 1987; Middelburg et al., 1997; Kemp et al., 2010). Station 4 had an average δ¹³C value of –28.0‰ ± 0.2, indicating a C₃ dominated salt-marsh plant community, which corresponds to Station 4’s Phragmites australis vegetation. Kemp et al. (2012b) found slightly less depleted bulk sediment δ¹³C values for brackish transition marsh zones vegetated by Phragmites australis, Iva fructescens, and Typha spp. in southern New Jersey, ranging from –27.0‰ to –22.0‰. In Massachusetts, Middleburg et al. (1997) recorded a bulk sediment δ¹³C value of -24.5‰ in a salt marsh-upland transition zone vegetated by Phragmites australis, Typha spp., and Scirpus spp.

The temporal and small-scale spatial variability observed at all four stations was very minimal and δ¹³C values always remained within the range for C₄ (Stations 1, 2, and 3) or C₃ (Station 4) plant communities. The stability in δ¹³C values suggests inconsequential influence from seasonal, interannual, or small-scale spatial changes. Milker et al. (2015) also found no significant influence of interannual variations in bulk salt-marsh sediment δ¹³C values from three sampling periods over two years.

In order to use stable carbon isotope geochemistry as a sea-level indicator, it is necessary to understand the relationship between bulk sediment δ¹³C values and tidal elevation (Shennan, 1986; van de Plassche, 1986; Khan et al., 2015). For example, the transition between C₃ and C₄ dominated salt-marsh plant communities has been shown to act as the boundary for the mean higher high water (MHHW) tidal datum on the U.S. mid-Atlantic coast (e.g., Middelburg et al., 1997; Johnson et al., 2007; Kemp et al., 2012b). From three sites in southern New Jersey, Kemp et al. (2012b) found that δ¹³C values of salt-marsh bulk sediment >–22.0‰ were found above MHHW, while values <–18.9‰ were found between mean tide level (MTL) and MHHW. Monitoring Station 4 is consistent with these findings since it has δ¹³C values >–22.0‰ and is found above local MHHW. However, Monitoring Stations 1, 2, and 3 in Spartina marshes are also found above their local MHHW but have δ¹³C values <–18.9‰. Site-specific ecological conditions could alter the distribution of dominant vegetation as it relates to a tidal datum (e.g., McKee & Patrick, 1988; Kemp et al., 2012b, 2017). For example, McKee & Patrick (1988) observed that elevation limits of Spartina alterniflora do not always correspond to a consistent elevation relative to a tidal datum. The variation in Spartina alterniflora vertical distribution was attributed primarily to differences in tidal range; however, local differences in salinity, nutrients, or physical disturbance may also have an influence (McKee & Patrick, 1988).

In studies using salt-marsh foraminifera as proxies to reconstruct RSL, it has been debated whether to use the live, dead, or total assemblages (e.g., Scott & Medioli, 1980b; Murray, 2000; Murray & Bowser, 2000; Culver & Horton, 2005). However, it has also been noted that because live assemblages are typically so low in number, they do not significantly contribute to or alter the total assemblage when combined with the dead assemblage (Buzas, 1965). At the four high marsh stations, the standing crop of live foraminifera fluctuated over space and time, but the overall assemblages and dominant foraminiferal species remained consistent both temporally and on small spatial scales. Further, although the live assemblages occurred in much smaller numbers than the dead assemblages and only accounted for <1 to 25% (average of 12%) of the total assemblages across all four stations over the study time period, the live assemblages also correlate with the dead assemblages. All of the species found in the live assemblages occurred in the dead assemblages and only occurred in greater abundances since the dead tests accumulate over time with many generations.

The results from this monitoring study showcase high-marsh environments where the live, dead, and total assemblages are all consistent. Our PAM analysis shows four clusters corresponding to the four monitoring stations regardless of whether the live, dead, or total assemblage is used. The live and dead assemblage analysis only had 13 and 12 samples, respectively, that were not assigned to the group corresponding to the station from which they were sampled. Analysis with the total assemblage reduced that number to 8 of the 184 samples incorrectly assigned, perhaps suggesting that the total assemblage is the best representation of modern foraminiferal assemblages. Further, analysis of the environmental controls on foraminiferal distributions showed that live, dead, and total assemblages all form elevation zones. Even though we only used four different stations across a salinity gradient (studies often use a transect across elevation), we found that elevation was the dominant controlling factor, followed by salinity, no matter which assemblage was used. These results highlight the ongoing utility of high marsh sedimentary environments and foraminifera assemblages to reconstruct sea-level change due to the consistency of high marsh foraminiferal assemblages at individual locations and their strong vertical zonation (Scott, 1976; Scott & Medioli, 1978, 1980a). We recommend the ongoing use of either dead or total foraminiferal assemblages for sea-level reconstructions due to the consistency of their assemblages through space and time, in addition to the near identical relationship between the dead and total assemblages and the measured environmental variables.

Salt-marsh foraminifera can be used as a reliable proxy to reconstruct relative sea level if their distributions and relationship to environmental controls is well understood in the modern environment.

We monitored live foraminiferal distributions from four high marsh monitoring stations in southern New Jersey over three years to evaluate small-scale spatial and temporal variability. The live foraminiferal assemblages remained consistent over space and time, even through extreme weather events, and the assemblages formed four site-specific clusters corresponding to each monitoring station. The live foraminiferal assemblages were predominantly controlled by elevation, with secondary influences from salinity and substrate.

We found that stable carbon isotopes, which can be used as a secondary proxy in relative sea-level reconstructions, also remained consistent on small-scale spatial and temporal scales. However, we note that to use carbon isotopes as a secondary proxy, it is important to fully understand the local relationship between stable carbon isotope values and tidal elevation as it may not be consistent regionally.

The live foraminiferal assemblages correspond to the previously investigated dead assemblages from the same sites (Walker et al., 2020), as well as the total (live + dead) assemblages, with the live foraminifera only occurring in smaller numbers compared to the dead. Live, dead, and total assemblages all formed four distinct clusters corresponding to the monitoring stations from which they were sampled. Additionally, the live, dead, and total assemblages each formed elevation zones. The similarities among the assemblages and consistency over space and time emphasize the reliability of salt-marsh foraminifera as a sea-level indicator.

The authors thank Isabel Hong, Ane Garcia-Artola, Andra Garner, and Tina Dura for their assistance in the field. JSW was funded by the David and Arleen McGlade Foundation, a Cushman Foundation for Foraminiferal Research Student Research Award, a Rutgers University Marine Field Station Graduate Student Award, and the U.S. National Science Foundation (award OCE-1804999). TAS and BPH were funded by the Ministry of Education Academic Research Fund MOE2019-T3-1-004, the National Research Foundation Singapore, and the Singapore Ministry of Education, under the Research Centres of Excellence initiative. Work by DB on this project was supported by the Harold F. Alderfer Fund for Environmental Studies at Bryn Mawr College. The authors acknowledge PALSEA (Palaeo-Constraints on Sea-Level Rise), a working group of the International Union for Quaternary Sciences (INQUA), and Past Global Changes (PAGES), which in turn received support from the Swiss Academy of Sciences and the Chinese Academy of Sciences. This work is a contribution to IGCP Project 725 ‘Forecasting Coastal Change’. This work is Earth Observatory of Singapore contribution 488. Appendices B and C can be found linked to the online version of this article.

Original references to the taxa identified to the species level.

• Ammoastuta inepta Cushman and McCulloch: Ammoastuta inepta Cushman and McCulloch, 1939, p. 89, pl. 7, fig. 6.

Remarks: One of the most abundant species in the live and dead assemblages from this study and Walker et al. (2020). A. inepta was found in high abundances in the low salinity, brackish high marsh-upland transition environment of Station 4 and has been found to be a dominant species in similar low salinity environments on the U.S. Atlantic coast, such as in New Jersey and North Carolina (Scott et al., 2001; Culver and Horton, 2005; Kemp et al., 2009; Kemp et al, 2013).

• Arenoparella mexicana Kornfeld: Trochammina inflata (Montagu) var. mexicana Kornfeld, 1931, p. 86, pl. 13, fig. 5.

• Balticammina pseudomacrescens Bronnimann, Lutze, Whittaker: Balticammina pseudomacrescens. Bronnimann, Lutze, Whittaker, 1989, 41: 167-177.

Remarks: One of the most abundant species in the live and dead assemblages from this study and Walker et al. (2020). Previously grouped under Trochammina macrescens, B. pseudomacrescens is differentiated from J. macrescens (de Rijk, 1995) as the distributions of the two species are spatially variable (Gehrels and van de Plassche 1999). B. pseudomacrescens is a dominant species at both low-salinity Stations 2 and 3, though Kemp et al. (2012) found lower abundances in the high marsh of other areas in New Jersey. Further, B. pseudomacrescens has not been recorded in North Carolina or Virginia (Spencer, 2000; Kemp et al., 2009; Wright et al., 2011), but has been found in New England and Newfoundland, Canada in varying abundances (de Rijk, 1995; de Rijk and Troelstra, 1997; Gehrels and van de Plassche, 1999; Edwards et al., 2004; Wright et al., 2011).

• Jadammina macrescens Brady: Trochammina inflata (Montagu) var. macrescens Brady, in Brady and Robertson, 1870, p. 47, pl. 11, figs. 5a-c.

Remarks: One of the most abundant species in the live and dead assemblages from this study and Walker et al. (2020). J. macrescens was the dominant species in the low salinity, brackish high marsh-upland transition environment of Station 4. Beginning with the work of Scott and Medioli (1978, 1980a), salt-marsh foraminifera assemblages dominated by J. macrescens have been considered the highest elevational zone and have continued to be found in maximum abundances in high marsh-upland transition environments with low salinities (e.g., Spencer, 2000; Nikitina et al., 2003; Robinson & McBride, 2006; Horton & Culver, 2008; Kemp et al., 2009).

• Milliammina fusca Brady: Quinqueloculina fusca Brady, in Brady and Robertson, 1870, p. 47, pl. 11, figs. 2–3.

• Milliammina petila Saunders: Milliammina petila Saunders, 1958, p. 88, pl. 1, fig. 15.

• Pseudothurammina limnetis Scott and Medioli: Thurammina (?) limnetisScott and Medioli, 1980b, p. 43, 44, pl. 1., figs. 1–3.

• Siphotrochammina lobata Saunders: Siphotrochammina lobata Saunders, 1957, p. 3, pl. 9, figs. 1,2.

• Tiphotrocha comprimata Cushman and Bronnimann: Trochammina comprimata Cushman and Bronnimann, 1948b, p. 41, pl. 8, figs. 1–3.

Remarks: One of the most abundant species in the live and dead assemblages from this study and Walker et al. (2020). T. comprimata was abundant at both low-salinity Stations 2 and 3. T. comprimata has previously been associated with lower salinity environments on the U.S. Atlantic coast (de Rijk, 1995; de Rijk & Troelstra, 1997; Kemp et al., 2013). However, it has also been noted as an indicator of the transition from low to high marsh as abundances decrease (Scott and Leckie, 1990) or reaching peak abundances before falling in the highest marsh environment, perhaps suggesting a preference for somewhat saline conditions (Edwards et al., 2004).

• Trochammina inflata Montagu: Nautilus inflatus Montagu, 1808, p. 81, pl. 18, fig. 3.

Remarks: One of the most abundant species in the live and dead assemblages from this study and Walker et al. (2020). T. inflata was the dominant species at Station 1, the highest salinity site. T. inflata has previously been recognized as a dominant species in the middle and high marsh at sites on the U.S. Atlantic coast (e.g. Hippensteel et al., 2000; Kemp et al., 2009; Kemp et al., 2012) and T. inflata-dominated assemblages have also been associated with higher salinity sites, such as in North Carolina, (Kemp et al., 2009).

• Trochammina ochracea Williamson: Rotalina ochracea Williamson, 1858, p. 55, pl. 4, fig. 112.

Live foraminifera counts over the three-year study period for all four monitoring stations and coefficient of variation among replicate samples. Appendix B can be found linked to the online version of this article.

Live foraminifera counts at each monitoring station with variability in local climate indices (Fig. 2), including air temperature, precipitation, and seasonal changes in mean sea level recorded by the Atlantic City tide gauge. Appendix C can be found linked to the online version of this article.