An intertidal mudflat (the Eddy, 43.994827, and −69.648632) located within the Sheepscot River Estuary Maine, United States of America was the subject of the current study (Figure 1; Stickney, 1959). As defined by Stickney the Eddy is classified as an upper estuary mudflat with seasonal variations in salinity and temperature that receives a significant influence from brackish and/or freshwater riverine input (Stickney, 1959). From November 2015 to November 2016 sediment cores were sampled (n = 11); for porewater and solid phase geochemistry in combination with microbial community structure (n = 10) at 1 cm depth intervals up to a depth of 10 cm. Due to circumstances beyond our control cores were not collected in August and October 2016. We took care not to take monthly samples in areas disturbed by our previous sampling events. This site was not impacted by local worm and shellfish harvesters as evinced by a lack of visible digging tracks in the surface sediments (Figure 1A), or direct observations of their presence. A wooden post and a washed-up bundle of lobster traps served as a reference point for sampling throughout the study (Figure 1B). At high tide, this area of the mudflat was covered with approximately one meter of water, dependent on tidal amplitude. At low tide, we successively took samples in the same horizontal distance from shore to reduce possible longitudinal effects on the mudflat macrofaunal populations. This mudflat lies adjacent to a tidally flushed salt marsh to the east, and faces south-southwest to the main stem of the Sheepscot River. The salinity and temperature varied throughout the year, but the overlying water was never frozen in the winter months of 2015–2016. Salinity, temperature and spatial/temporal data are presented in Supplementary Table 1.

Sediment cores (approximately 15–20 cm in depth) from the Eddy were taken with a Plexiglas core barrel (inner diameter, 7.5 cm) that had predrilled 4 mm holes drilled (sealed with electrical tape before sampling) at 1 cm intervals for Rhizon CSS samplers with a 5 cm length and 0.15 μm pore size (Rhizosphere Research Products, Wageningen, Netherlands; Seeberg-Elverfeldt et al., 2005). Porewaters were extracted by inserting Rhizons into the predrilled holes starting from the top of the core and working down to the final depth interval (9–10 cm). Sterile 10 mL syringes were used to pull negative pressure on the Rhizons, and the plunger was held in place for 1–2 h with a wooden block. After ∼10 mL of porewater had been extracted, the Rhizons were removed, and the porewater was dispensed into a sterile 15 mL tube in a glove bag containing an atmosphere of N₂/H₂ (95/5%) or under ambient atmospheric conditions. Ferrous iron [Fe(II)] was immediately preserved by pipetting 250 μL of pore water into 250 μL of Ferrozine buffer (10 mM Ferrozine in 50 mM HEPES buffer, pH = 7; Stookey, 1970; McBeth et al., 2011). We did not observe substantial differences in porewater Fe(II) concentrations when rhizon-derived porewater was rapidly preserved, <10 s from rhizon collection to Ferrozine treatment under ambient air conditions, compared to working in the glove bag, thus after the April sampling all samples were collected under ambient atmospheric conditions. The Ferrozine-porewater solution was measured at 562 nm on a Multiskan MCC plate reader (Thermo Fisher) and absorbance values were calculated to concentrations by plotting them on a standard curve of known Fe(II) concentration solutions ranging from 0 to 1,000 μM prepared from a 10 mM solution of FeCl₂ in 0.1 N HCl.

Sediment cores used for porewater extraction were extruded and sliced into 1 cm sections with a brass wire and dried in an oven at 70°C for 24 h. Highly reactive iron oxides (ferrihydrite and lepidocrocite as well as iron monosulfides that were oxidized during this process) were extracted for 48 h in a 50 mL falcon tube with 10 mL of hydroxylamine-HCl (1 M) in 25% acetic acid (w/v; Poulton and Canfield, 2005). The total extracted Fe was measured with a Ferrozine solution as described above, except the samples were diluted 1:100 in distilled water due to high absorbance values of undiluted samples.

Sub-cores were sampled simultaneously with porewater extraction using 1 mL syringes with cut off tips. The 1 mL syringe corers were inserted into pre-drilled 8 mm holes at a right angle to the Rhizon sampling holes, and gently pressed into the side of the core while drawing back on the plunger to remove approximately 0.2–0.5 cc of sediment. The 1 mL sub-cores were then frozen at −80°C until DNA extraction. DNA was extracted from ∼0.1–0.2 g of wet sediment using the DNeasy PowerSoil Kit (Qiagen) with a modification of the removal of 200 μL of the bead solution and addition of 200 μL of 24:25:1 phenol:chloroform:isoamyl alcohol (Sigma Aldrich); the manufacturers protocol was then followed for the rest of the extraction, with the exception of the DNA column elution step, which was eluted with 100 μL of molecular grade water instead of the elution buffer from the kit. The DNA was quantified using a Qubit 3.0 fluorometer with the Qubit dsDNA quantification kit and protocol (Thermo Fisher Scientific).

The 16S rRNA gene was amplified from the 1 cm interval DNA extractions with the 515F/926R primer pair and the 18S rRNA gene was amplified with the E572F/E1009R primer pair at the Integrated Microbiome Resource (Dalhousie University, Halifax, NS, Canada). Paired end reads were assembled and classified using a mothur (Schloss et al., 2009) pipeline from the mothur MiSeq standard operating procedure (Kozich et al., 2013; Scott et al., 2017). Operational taxonomic units (OTUs) were clustered at either 85% or 97%, which generated 12,437 and 351,057 OTUs, respectively. The 85% identity cutoff, which corresponds roughly to the class-order taxonomic level (Yarza et al., 2014), was chosen for comparative analysis to limit the number of clusters generated for comparisons and lend clarity to data representations. While the 97% cut-off did sub-divide some clusters among the top 100 OTUs (based on relative abundance), we did not note any substantial differences that would change the functional inferences made from OTU analysis. Relative abundances were utilized to build a heatmap in RStudio (RStudio Team) of the top OTUs that were greater than 0.5% abundance in each sample.

The 18S rRNA gene data was processed using mothur (v1.43) following the mothur MiSeq SOP data analysis pipeline (Kozich et al., 2013), but referencing the pR2 database (v4.12) to assign eukaryote taxonomies (Guillou et al., 2013). Eukaryote OTUs were determined at a level of 97% sequence similarity (Hu et al., 2015; Pasulka et al., 2019), which is less conservative than some thresholds proposed for rRNA genes (Caron et al., 2009), and more conservative than others (Callahan et al., 2017; Edgar, 2018). Counts of eukaryotic lineages of interest (i.e., diatoms, polychaetes, molluscs, nematodes, and protozoa) were filtered from the taxonomic data and used to calculate relative abundances for all samples across the year.

One way to organize and find signals in complex ecological communities is to cluster them into coherent “ecostates” and then to map these ecostates back onto time and space to highlight patterns (Record et al., 2017). Each ecostate represents a consortium of taxa that clusters together consistently throughout the sampling program, and maintains the same respective relative abundances of taxa (O’Brien et al., 2016). For this analysis, samples were clustered using average linkage hierarchical clustering of the Bray–Curtis dissimilarity as implemented in the R package vegan (Oksanen et al., 2019). The clustering algorithm used bacterial and archaeal OTUs clustered at 85% identity (n = 12,437 OTUs), since, as discussed above, this identity was a compromise between capturing overall community diversity without overwhelming the analysis with a large number of OTUs. Ecostates were then numbered and mapped back onto the spatiotemporal sampling grid in R Studio (R Studio Team).

The pore water and solid phase sedimentary geochemical data can be found in the Biological and Chemical Oceanography Data Management Office (BCO-DMO) database under the DOI:10.1575/1912/bco-dmo.737962.1. The 16S and 18S rRNA gene raw read data can be found under the Bioproject ID PRJEB36098. All sedimentary geochemical and DNA data can also be found in Supplementary Table 1; all geochemical and DNA profile plots in the main text can be reproduced from this data set.

The dissolved ferrous iron [Fe(II)] in porewaters exhibited both spatial and temporal variability over the sampling period, which was approximately every month during the 2015–2016 annual cycle (Figures 2A,B). Porewater Fe(II) concentrations peaked at a maximum average concentration of approximately 120 μmol L^–1 at a depth interval between 3–4 cm, which is a consistent observation with many porewater Fe(II) profiles from coastal and continental shelf systems (Aller, 1980; Canfield, 1989; Hines et al., 1991; Thamdrup et al., 1994; Severmann et al., 2010; Antler et al., 2019), which are more often than not, bioturbated. This peak in porewater Fe(II) is consistent with the activity of iron-reducing microorganisms at this depth interval throughout the study. Throughout the year, the porewater Fe(II) concentration maxima remained in the top 5 cm of sediment, and migrated only once to a depth of 6.5 cm in September (Figure 2B). These Fe(II) porewater maxima suggest that the most active iron cycling zone occurs within the upper 5 cm of sediment. Concentrations of Fe(II) decreased to approximately a seasonal average of 100 μmol L^–1 at the maximum depth interval sampled between 9–10 cm (Figure 2A). The decline in Fe(II) with depth is likely due to the reaction with hydrogen sulfide produced by sulfate-reducing bacteria in sediments (Jørgensen, 1982) as well as around macrofaunal burrow walls (Bertics and Ziebis, 2010), which results in the precipitation of iron sulfides in sediments (Berner, 1984). There were consistent pools of porewater Fe(II) over the annual cycle, which exhibited large spatiotemporal changes (Figure 2B). The large variability over space and time could be due, in part, to the heterogeneity of natural populations of macrofauna and microorganisms in the mudflat, which we attempted to minimize by sampling cores within a consistent region over the course of the study (Figure 1). Other studies have shown consistent variability in porewater Fe(II) with space and time (Hines et al., 1991; Thamdrup et al., 1994; van de Velde et al., 2020), while another recent study showed contemporaneously collected salt marsh sediment samples show consistency in geochemical profiles despite large geographic variation (Antler et al., 2019). The abundance and activity of macrofaunal populations likely exerts a primary control on porewater Fe(II) concentrations, due to their sediment mixing, irrigating, and feeding behavior (Aller, 1994; van de Velde et al., 2020).

Generally, the porewater Fe(II) concentrations increased with rising temperature (Figure 3A). This might be indicative of elevated sedimentary carbon mineralization in summer months by iron-reducing bacteria (Canfield, 1989; Thamdrup et al., 1994). Sulfide production via sulfate reduction could also lead to iron reduction (Luther and Church, 1988); however, we did not detect hydrogen sulfide in these sediments, indicating anaerobic respiration coupled to iron reduction is, on balance, an important process. The increases in porewater Fe(II) with elevated temperatures between the months of February and July has also been observed in bioturbated sediments in the Great Bay New Hampshire, United States (Hines et al., 1984). The supply of terrigenous organic matter to this mudflat occurs mainly in the spring that corresponds with the thawing land mass and riverine run-off (Stickney, 1959). Additional organic inputs come from settling of detrital material from the water column, and benthic mats of phytoplankton that can be transported into the sediments via the bio-irrigating and bio-mixing activites of macrofauna (Christensen et al., 2000). It is reasonable that this accounts in part for the high mobility and turnover of porewater Fe(II) due to the behavior of irrigating animals influencing microbial iron cycling processes (van de Velde et al., 2020). This relationship likely explains the large variance of the porewater Fe(II) pools over the 13 month long study, which is overall influenced by seasonal variations in sediment temperature (Figure 3A) and perhaps availability of macrofaunal food sources such as phytoplankton (discussed in detail below).

Similar to measurements of Fe(II), there was a continuous and variable pool of highly reactive Fe (Fe_HR) with respect to depth over the seasonal cycle (Figures 2C,D). The concentration of Fe_HR declined slightly with depth (Figure 2C) perhaps due to the interaction with hydrogen sulfide. The most significant change in the Fe_HR pool was in the warmer months (July and September) into the summer and early fall, with a nearly two-fold increase in the Fe_HR pools compared to the other months over the year (Figure 2D). We hypothesize that this change is due to several key biogeochemical processes, which are controlled by the master variable, temperature, over the seasonal cycle. Changes in irrigation rates that likely occur throughout the year (Fang et al., 2019) will have profound effects on the exchange of solutes as well as the control of the pools of sedimentary Fe_HR—when irrigational activity increases, concomitant microbial oxygen consuming process, such as microaerobic bacterial iron oxidation and aerobic carbon mineralization increase. When irrigation activity is low, the pools of sedimentary Fe_HR are depleted in the winter months. This is likely due to interaction Fe(II) with hydrogen sulfide produced by sulfate reducing microorganisms, resulting in the formation of more iron sulfides. An increase in the pool of Fe_HR in the summer has important biogeochemical ramifications, including: the potential to enhance the iron flux from the sediment by stimulating iron reduction (Severmann et al., 2010), and by providing fresh and reactive substrates for microbial iron reduction, which may effectively compete with sulfate reduction in bioturbated sediments (Froelich et al., 1979; Canfield, 1989).

The extractable quantity of sedimentary DNA (μg g^–1 sediment) decreased from approximately 10 μg g^–1 at 1 cm to ∼7.5 μg g^–1 at 10 cm (Figure 2E). Assuming that DNA pools are approximate to microbial biomass, these values infer that microbial abundances decrease with depth, which is consistent with observations in both shallow and deep-sea sediments (Kallmeyer et al., 2012; Chen et al., 2017). Higher microbial biomass near the sediment-water interface is consistent with this being the zone where the most labile organic matter is deposited at the sediment surface, while organic matter becomes more refractory with depth. However, macrofauna may confound these interactions by actively exchanging and flushing bottom water into deeper sediments, as well as by bio-mixing, thereby introducing more labile organic matter, either dissolved or particulate, deep into the sediments. Across the year, DNA pools were distinct between fall and winter months versus summer months (Figure 2F). These variations may be caused by fall phytoplankton blooms that are stored in the sediments over the winter months, and decrease in the summer months due to the activity of microbial mineralization of detrital organic matter, including degradation and metabolism of exogenous DNA. This is evident when sedimentary DNA pools are plotted as a function of temperature, which is negatively correlated; as temperature increased, DNA pools decreased (Figure 3C).

There were OTUs that could be assigned to the iron biogeochemical cycle. The iron-oxidizing Zetaproteobacteria formed an OTU (Otu00040) that was 0.28% median seasonal abundance (range 0.01–2%), which is a typical relative abundance for this class of bacteria in marine sediments (Rubin-Blum et al., 2014; McAllister et al., 2015; Laufer et al., 2016, 2017; Beam et al., 2018). The Zetaproteobacteria relative abundances changed with respect to depth (Figure 4A), but did not exhibit pronounced seasonal variation except in one sample in September that appeared to correlate with increased Fe_HR production (Figures 2D, 4B). The relative abundance of Zetaproteobacteria ranges from 0.25–0.75%, in the upper 5 cm of sediment compared to the deeper sediment, with ranges from around 0–0.25%. The highest relative abundances were at 3–4 cm (Figure 4A), which has been previously documented from a depth profile of 16S rRNA genes from sediment and iron oxide macrofaunal burrow walls at this mudflat (Beam et al., 2018). The peak abundance of Zetaproteobacteria also coincides with the peak in porewater Fe(II) and Fe_HR from these sediments (Figures 2A,C). These Zetaproteobacteria profiles are consistent with them being most abundant in the active iron cycling zone throughout the year. It is also important to note that in this study we analyzed bulk sediment where we would expect lower overall numbers of Zetaproteobacteria, Previous studies, including one from this same site, have shown their abundance increases significantly in iron-oxide coated worm burrow walls, where they are likely involved in Fe(II) oxidation coupled to oxygen reduction and carbon dioxide fixation (McAllister et al., 2015; Beam et al., 2018). Consistent with these observations is that increased dark CO₂ fixation with concomitant Fe(II) oxidation has been observed in burrow walls of the polychaete Neries diversicolor (Reichardt, 1986), which was the most abundant polychaete at the Eddy. Thus, dark CO₂ fixation by iron-oxidizing Zetaproteobacteria may contribute to the autochthonous organic matter pools in coastal sediments, similar to that observed in autotrophic sulfur-oxidizing Gammaproteobacteria (Dyksma et al., 2016).

Microorganisms that catalyze the reduction of iron are phylogenetically diverse (Melton et al., 2014) and, particularly in marine sediments, are more difficult to identify with 16S rRNA gene profiling. However, culture-based studies from marine sediments have suggested that members of the Deltaproteobacteria family Desulfuromonadaceae can reduce iron, as well as other some other Deltaproteobacteria lineages (Vandieken et al., 2006), which have been subsequently identified in coastal sediments (Reyes et al., 2017; Otte et al., 2018; Buongiorno et al., 2019). We identified one OTU containing members of the Desulfuromonadaceae family (Otu00028) that was relatively abundant seasonally accounting for a median relative OTU abundance of 0.80% (range = 0.40–1.7%; Figures 4C,D). The Desulfuromonadaceae OTU was also most abundant in the active iron cycling zone (Figure 6C), which appears to mirror the porewater and highly reactive iron profiles (Figures 2A,C). The Desulfuromonadaceae abundances fluctuated temporally (Figure 4D) perhaps in response to organic matter and/or iron oxide availability. Correlations of Desulfuromonadaceae with porewater and highly reactive iron geochemistry were not indicative of a strong relationship (not shown), but this does not preclude a causal relationship. Canonical iron-reducing bacterial genera such as Geobacter and Shewanella-like OTUs were extremely rare or below detection in these sediments, similar to findings from coastal sediments near Denmark and Sweden (Reyes et al., 2017; Otte et al., 2018; Buongiorno et al., 2019).

The sedimentary microbial communities were dominated by bacteria (total OTU abundance = 97.5%), with minor contributions by archaea (total OTU abundance = 2.5%)—consistent with the observations of Chen et al. (2017) and Kim et al. (2008) that bacteria dominate over archaea in bioturbated sediments. Two of the three most abundant bacterial OTUs throughout the seasons were related to the Deltaproteobacteria, and one was associated with the Gammaproteobacteria (Figure 5). Deltaproteobacteria and Gammaproteobacteria are commonly observed as abundant members of coastal sediment microbial communities (Kim et al., 2008; Lasher et al., 2009; Wang et al., 2012; Chen et al., 2017; Cleary et al., 2017; Qiao et al., 2018). The two Deltaproteobacteria OTUs were related to the families Desulfobulbaceae (Otu00001; median seasonal abundance = 6.4%; range = 3.0–8.8%) and the genus Desulfococcus in the family Desulfobacteraceae (Otu00002; median seasonal abundance = 5.7%; range = 1.6–9.8%), which are likely involved in sulfur oxidation and sulfate reduction, respectively, in these sediments. The Desulfobulbaceae is a diverse family within the Deltaproteobacteria that contains the so-called “cable bacteria” that transport electrons from sulfur oxidation over centimeter scale distances and catalyze oxygen reduction at sediment-water interfaces (Nielsen et al., 2010; Meysman et al., 2015; Nielsen and Risgaard-Petersen, 2015). We have observed on many occasions the presence of cable bacteria-like filaments from this field site (Figure 1F), and hypothesize that they are functional in these bioturbated sediments, which has been recently demonstrated (Aller et al., 2019), and may impact iron dynamics (Sulu-Gambari et al., 2016). Members of the Desulfobulbaceae are also involved in sulfur disproportionation, which could also be an important process in these sediments, and would impact the iron cycle (Finster, 2011).

To gain a more holistic understanding of the microbial community dynamics, each discrete spatiotemporal sample (n = 100) was clustered based on the relative abundance of OTUs at 85% 16S rRNA gene identity. The composition of the sedimentary microbial community did not exhibit large fluctuations with the changing seasons in this bioturbated mudflat (Figures 5, 6). Although some individual OTUs exhibit spatial and temporal variability, the community was generally comprised of approximately 20–25 dominant OTUs (Figure 5). Seven main clusters predominated the various communities with respect to time and space (Figure 5). These clusters comprise different “ecostates” that were identified at various depths and times during the season. Ecostates are operationally defined as a collection of inferred taxonomic groups from distinct ecological samples, and they can be used to track changes in ecological communities through space and time (O’Brien et al., 2016). There is a strong depth stratification of ecostates that persists temporally (Figure 6). Ecostate 2 typically dominates all depth intervals <5 cm at all months. We speculate that this ecostate, which has greater abundances of putatively aerobic microbes (Supplementary Table 3), experiences more oxic conditions, partly due to diffusion of O₂ from the sediment air/water interface, but primarily driven by biodiffusive transport by animals that create sporadic “oxygen oasis” in these otherwise anoxic sediments (Glud, 2008; Volkenborn et al., 2010). Ecostate 1 was usually confined to deeper sediments >5 cm where lower oxygen conditions would be prevalent (Glud, 2008; Figure 6). The presumed presence of labile organic matter in the upper 5 cm of sediment also likely contributes to the separation of ecostates 1 and 2 as has been observed between other surface and subsurface coastal sediments (Cleary et al., 2017; Qiao et al., 2018). These ecostates are relatively stable throughout the seasons (Figure 6), and appear to be unaffected by observed seasonal changes in environmental conditions and other stochastic events likely occurring in the intertidal zone. More than half of the samples at the sediment-water interface were classified as ecostate 3, which suggests some unique properties of this interface that contributes to the observed community signal (Figures 5, 6). This consistent pattern of stratification was only disrupted once, in the April profile. This profile differed markedly, with three unique ecostates (4, 5, and 6) that were not observed in any other month. These ecostates contained a cluster of high relative abundances of OTUs, described below, that only comprise a minor community component of the other ecostates (Figure 6). The occurrence of unique ecostates in April could reflect a community response to ecosystem disturbance or a stochastic assembly of horizontal spatial heterogeneity, which would both imply resilience of the system.

Approximately 50% of the taxa assigned to each of the different ecostates, either at the order, family, or genus level, could be associated to known physiological representatives that are presumed to be either predominately anaerobic or aerobic. For example, all known relatives of the family Desulfobulbaceae, or genus Desulfococcus, predominant in Ecostates 1 and 2 are known to be anaerobes; while known relatives of the family Flavobacteriaceae, and family Marcinicellaceae (both abundant in ecostate 3) are obligate aerobes. From this we surmise that the overall dominant ecostates 1 and 2 had relatively high abundances of putative anaerobes, compared to putative aerobes (Supplementary Table 2), while ecostate 3, found primarily in the top 1 cm (Figure 5) had higher relative abundances of putative aerobes, and corroborates the availability of oxygen declining with depth (Glud, 2008). In general, we found members of the phyla Bacteroidetes more prevalent in ecostates 2, 3, and 6, all associated the upper 5 cm of the sediments than in deeper ecostates 1 and 7. The Bacteroidetes are known for their involvement in the degradation of more complex organic matter, either aerobically or anaerobically (Kirchman, 2002; Hahnke et al., 2016). This is consistent with the most labile organic matter being present in the upper layers of sediment, in part driven by spring and fall phytoplankton blooms. Ecostates 4, 5, and 6 were only present in April, and OTUs that increased in abundance were Flavobacteriaceae (Otu00006), Rhodobacteraceae (Otu00022), Vibrionales (Otu00051), Helicobacteraceae (Otu00052), unclassified Bacteroidetes (Otu00009), Clostridiales (Otu00066), and Alteromonadales (Otu00012) OTUs (Supplementary Table 1). These OTUs may be involved in the degradation of phytoplankton-derived organic matter—relatives of some of these OTUs have been observed to change abundance in response to phytoplankton-derived organic matter as well as plant-derived organic matter (Raulf et al., 2015). OTUs that were elevated in April were similar to the elevated OTUs in ecostate 3, which further supports the hypothesis that these OTUs are potentially involved in organic matter-specific responses such as carbon mineralization. Although these microbial responses to phytoplankton blooms are tentative, it would be worth exploring in other coastal environments with respect to seasonal variation.

An important ancillary question related to our focus on iron-cycling communities, is the role for pelagic-derived phytoplankton, such as diatoms, to be transported to these intertidal sediments by sedimentation and/or bio-mixing activity of suspension feeding macrofauna. These phytoplankton may act as an important source of labile organic matter in sediments at different times of the year (Kristensen, 2001; Kristensen and Mikkelsen, 2003). To address this, we used 18S rRNA gene sequencing to determine the presence and relative abundance of phytoplankton and other eukaryotes in these sediments over the seasonal cycle. We recovered 18S rRNA gene amplicons from 99/100 samples over the different seasons from depths of 1–10 cm. However, the total number of quality filtered sequences from each sample varied significantly with a median of 2,636 sequences per sample and a range of 1–24,484 sequences. One explanation for this variability is that it reflects biological reality, in that there are simply few eukaryotes present in those low-yield samples, thus yielding few 18S rRNA gene sequences. Pasulka et al. (2019) utilized a similar approach to investigate benthic microbial eukaryotes at Guaymas Basin hydrothermal vents, and also recovered a highly variable numbers of sequences among samples. Further sequencing of samples from marine benthic habitats will be needed to understand the relationship between sequence yields, the biomass of microbial eukaryotes or meio-fauna in unprocessed samples, DNA template concentrations in extracted samples, and possible PCR inhibition in these complex and heterogeneous environmental DNA samples. It cannot be ruled out that some sequences are sourced from relict, extracellular DNA bound to sediment particles, and are not indicative of the conditions that prevailed during the course of this study. However, a recent paper looking at marine sediments found little influence of extracellular DNA on amplicon studies similar to these samples (Ramírez et al., 2018).

With these caveats in mind, we found a total of 119 eukaryotic OTUs present at a relative abundance of ≥0.1% of total sequences (range = 0.10–3.87%; Supplementary Table 4). Of these OTUs, 34 were identified as belonging to the two most abundant classes of eukaryotic marine phytoplankton, including the diatoms (Bacillariophyta = 18 OTUs) and dinoflagellates (n = 16 OTUs). For illustrative purposes, the spatiotemporal distribution of diatoms are shown in Figure 7. While there may be temporal patterns to the diatom distribution, further sampling and correlation to additional parameters like quantification of chlorophyll and cell counts in the overlaying water will be necessary to make these associations. The dinoflagellates (data in Supplementary Table 4) showed similar patterns. The fact that significant diatom abundance was found in the depth range of 2–5 cm, i.e., below the water-sediment interface, Figure 7, indicates the transport of diatoms, or diatom-DNA, into the subsurface sediments presumably via bio-mixing by marine worms. This would make an important source of labile carbon that would help sustain biogeochemical cycling of iron and sulfur, and could help explain the profiles of Fe(II) and Fe_HR discussed above.

Many of the most abundant phytoplankton OTUs were associated with benthic habitats (Supplementary Table 4). Diatoms and dinoflagellates were identified as four of the ten most abundant OTUs and included a dinoflagellate from the order Peridiniales (OTU0004, 3.87%), a raphid-pennate diatom (OTU0006, 3.24%), the diatom Psammodictyon sp. (OTU0009, 2.85%) and the dinoflagellate Islandinium minutum (OTU0010, 2.72%). The Peridiniales include a diverse group of dinoflagellates rendering the taxonomic designation of OTU0004 somewhat ambiguous. Nonetheless, OTU0004 is notable with respect to its prevalence across depth, with greater proportions found deeper in the sediments (Supplementary Table 4). Raphid-pennate diatoms include both pelagic and benthic species, but the motile benthic species are thought to be substantially more diverse than pelagic types (Mann et al., 2017), and are likely the type associated with OTU0006. The diatom Psammodictyon (OTU0009) is also associated with benthic habitats, living at the sediment water interface (Round et al., 1990). Islandinium minutum (OTU0010) is primarily associated with a particular “spiny” morphology of dinoflagellate cysts that are common in high latitude sediments (Potvin et al., 2013). The extent to which phytoplankton distributions in the sediments can be used as indicators of subsurface biomixing and ventilation remains to be tested under more controlled experimental conditions, but the possibility exists that these organisms could serve as useful biomarkers for this important biogeochemical process.

The 18S rRNA sequences also revealed the presence of Nereis (Hediste) diversicolor, which was the presumed dominant polychaete worm in this mudflat. This was the most abundant polychaete 18S rRNA gene sequence, accounting for 3.7% of the total 18S rRNA gene sequences (n = 14,149 sequences). Mollusca 18S rRNA gene sequences were not a very abundant component of any samples, although we detected 0.04% of the total 18S rRNA gene sequences (n = 148 sequences) related to the small, intertidal-dwelling clam, Macoma balthica, which is common to the Sheepscot Estuary and to North Atlantic intertidal sediments (Stickney, 1959). We have observed small clams similar to M. balthica from the Eddy on numerous occasions.

Nematodes can comprise a large fraction of the biomass of marine sediments (Heip et al., 1985) and are bewilderingly diverse (Bik et al., 2010). The 18S rRNA gene libraries revealed that nematodes comprised a median relative abundance of 21.9% with a range of 0–67.2% of the total 18S rRNA gene sequences. Their relative abundance was highest near the sediment-water interface, and abundance declined with depth (Figure 8). This finding corroborates observations primarily based from direct species counts and identification (Montagna et al., 1989). Nematode abundance and grazing activity have also been observed to increase bacterial activity and diversity in terrestrial sediments (Jiang et al., 2017). Nematodes play an important role in grazing of bacteria and other microbial eukaryotes in sediments (Pascal et al., 2008), and potentially play an important role in controlling iron biogeochemical dynamics.