Seasonal Dynamics of Microbial Diversity at a Sandy High Energy Beach Reveal a Resilient Core Community

Introduction

One third of earth’s ice free coastline are sandy beaches (Luijendijk et al., 2018). They provide a wide range of ecosystem services like nutrient cycling, water purification and protection from storm events (Nel et al., 2014). In general, sandy sediments are characterized by a high permeability (Huettel et al., 1998; Riedel et al., 2010) and low organic matter concentrations (Anschutz et al., 2009). Their high permeability and physical forcing by waves and tides lead to advective porewater exchange between beach sediments and the coastal ocean. The magnitude of exchange is higher than in less permeable sediments, where porewater transport is driven by slow diffusion (Webster and Taylor, 1992; Reckhardt et al., 2015). Beaches receive nitrogen and phosphorous from terrestrial runoff (Burnett et al., 2003), while seawater infiltration introduces labile organic matter, oxygen, sulfate (Huettel et al., 2014) and seasonally nitrate (Ahrens et al., 2020b). Moreover, beaches around the globe are affected by submarine groundwater discharge (SGD) (Moore, 1996; Paerl, 1997; Santos et al., 2009; Rodellas et al., 2014; Anschutz et al., 2016; Cho et al., 2018). SGD includes fresh terrestrial SGD as well as recirculated sea water. The volumetric flux has been estimated to exceed all riverine input into the world’s oceans (Moore et al., 2008; Kwon et al., 2014). Mixing of meteoric groundwater and tide-driven infiltrating seawater forms a subterranean estuary (STE). Steadily increasing consumption of electron acceptors with increasing groundwater age results in a horizontal redox zonation toward the discharge area at the low water line (LWL) (Reckhardt et al., 2015). Reduced chemical components such as ammonium or dissolved iron either end up in the coastal ocean or are re-oxidized and partly precipitated in surficial beach sediments and the nearshore water column (Riedel et al., 2011; Reckhardt et al., 2015; Beck et al., 2017; Kim et al., 2017). In this transition zone, advective mixing of terrestrial and marine constituents leads to high microbial activities (Anschutz et al., 2009; Gao et al., 2012; Charbonnier et al., 2013). As benthic microbial communities alter the nutrient and carbon pool, subterranean estuaries were described as natural bioreactors (Anschutz et al., 2009; Seidel et al., 2015). Remineralization products are either discharged and end up in coastal oceans, or are sequestered within beach sediments (Webster and Taylor, 1992; Moore, 1999; Kwon et al., 2014; Beck et al., 2017).

The global impact of SGD on a variety of biogeochemical cycles has led to increased scientific attention to investigate microbial diversity in STEs of sandy beaches. Most of these studies targeted sheltered, microtidal beaches, exhibiting microbial communities structured along horizontal and vertical redox gradients (McAllister et al., 2015; Hong et al., 2018). Others intensively analyzed microbial diversity of sandy beaches, yet often lack data for seasonality, geochemical parameters or spatial variability within the intertidal zone (Boehm et al., 2014). One study at a non-sheltered beach of Jeju Island by Lee et al. (2017) investigated the effect of changes in groundwater discharge flux on microbial diversity, yet only porewater samples were analyzed. They report different community structures between flood and ebb tide samples as a result of detachment and transport of subsurface bacteria from the aquifer, caused by high discharge velocities. Thus, they found periodical shifts in the community structure due to tidal cycles and SGD rate. Their findings support the hypothesis that SGD under high energy conditions has different effects on microbial diversity than at sheltered beaches.

In our study we complement the previous results reported by other studies, and explore further the specifics of microbial communities in high-energy beach subterranean estuaries, by using a multidisciplinary approach. We investigated microbial diversity patterns in combination with geochemical measurements over three seasons, covering the whole intertidal area down to a depth of 1 m at a sandy, high energy beach of Spiekeroog Island, Germany. For every sediment sample taken, porewater from the respective site and depth was analyzed, in order to corroborate geochemical and microbial diversity patterns. Simultaneously to our sampling campaigns, complementary studies were conducted on spatial and seasonal variations in porewater geochemistry. Waska et al. (2019b) found gradients in porewater chemistry to be strongest in the cross shore direction. Physico-chemically, the high water line (HWL) and the LWL were most distinct, possibly due to relatively stable in- and ex-filtrating conditions, respectively. Intermediate sites showed strong spatial and temporal heterogeneity, which was hypothesized to be an effect of intense morphological changes throughout the year. Ahrens et al. (2020b) found seasonal variations in seawater composition affecting the pathways of microbial organic matter degradation. In colder seasons, oxic and nitrate reducing conditions dominated the porewaters, while during warmer seasons, manganese and iron reduction became more important. This lead to a shift of the intertidal sediments from being a sink for dissolved inorganic nitrogen (DIN) during March to a net source in August. Thus, the strong dynamics altering the geochemistry of this environment may shape the microbial diversity to an unknown degree. So far, it is not known how the changes in salinity and chemical composition in the SGD influence the microbial communities. While previous studies on the SGD observed a distinct microbial community in the mixing zone at low energy beaches, we hypothesize, that those observations might not be valid for the high energy beach at Spiekeroog Island due to the strong physical disturbances.

Materials and Methods

Site Description

The study site is a sandy high energy beach at the seaward side of the barrier island Spiekeroog, located in the Southern North Sea (Figure 1A). The island has a freshwater lens which partly discharges within the intertidal area, exporting nutrients into the North Sea (Beck et al., 2017). The groundwater flow paths in the intertidal zone of Spiekeroog beach are influenced by characteristic topographic features (Figure 1B). These are the HWL and an intertidal ridge, which are zones of porewater infiltration due to their elevation above the mean sea-level. Located below the mean sea-level, the LWL and a runnel are porewater exfiltration sites (Waska et al., 2019b). These morphological features are not stationary, but are relocated and reshaped within the intertidal zone in a timescale of days to months by longshore currents and wave action (Anthony et al., 2005).

Figure 1. Location of the sampling site at Spiekeroog Island modified after Ahrens et al. (2020b). (A) The dashed line indicates the dimensions of the fresh water lens. The time series station in the tidal inlet (indicated by a star) monitors various sea water constituents. (B) The lower panel displays a cross-section through the beach, including hydrological flow paths and beach morphology. T1 and T2 indicate transects sampled across the intertidal area, the squares mark the sampling sites. MHWL, mean high water line; MSL, mean sea-level; MLWL, mean low water line.

Sediment and Porewater Sampling

Sampling was conducted during three campaigns in October 2016, March 2017, and August 2017 at the seaward beach of Spiekeroog Island (Figure 1A). Sediment cores were recovered with aluminum core liners along two replicate transects (distance: 60 m) between the HWL and LWL at fixed GPS positions. Four cores were taken at each transect at the specific topographical features: HWL, runnel, ridge and LWL. In October, an additional core (LWL b) was recovered from the LWL. Due to different tidal ranges this site could not be reached in later sampling campaigns. The cores were subsampled at 0, 10, 30, 50, and 100 cm depth using sterile cut-off syringes. Samples were stored frozen at −20°C until further processing. Porewater samples were obtained using stainless steel push-point samplers and pre-rinsed polyethylene syringes. Samples were processed according to the methods described in Ahrens et al. (2020b): O₂ was determined optically on-site with a hand-held flow cell. Porewater conductivity was measured on-site with hand-held multimeters. Both parameters were determined immediately after sampling. For nutrient, trace metal, and dissolved organic carbon (DOC) analyses, samples were filtered, preserved, and stored as outlined in Ahrens et al. (2020b) and Waska et al. (2019a). Whereas NO₃^–, NO₂^–, NH₄⁺, and Si were determined photometrically, Fe was analyzed with inductively coupled plasma mass spectrometry (ICP-MS). PO₄^3– was determined photometrically, following the methods of Itaya and Ui (1966) for samples below 2.1 μM and Laskov et al. (2007) for samples exceeding 2.1 μM PO₄^3–. DOC samples were treated and processed according to Waska et al. (2019b). The complete data set of O₂, NO₃^–, NO₂^–, NH₄⁺ and Fe data are freely available at the data publisher PANGAEA (Ahrens et al., 2020a). Data used for this study are compiled in the Supplementary Table S1.

DNA-Extraction, Amplification and Sequencing

Genomic DNA was extracted by using a modified phenol-chloroform protocol (Gabor et al., 2003; Lueders et al., 2004). PCR was performed in 50 μl setups. The master mix contained 30.4 μl PCR grade water, 10 μl Phusion-Polymerase Buffer, 1 μl (10 pmol/μl) of each primer published by Klindworth et al. (2013); forward: S-D-Bact-0341-b-S-17 (5′-CCT ACG GGN GGC WGC AG-3′), reverse: S-D-Bact-0785-a-A-21 (5′-GAC TAC HVG GGT ATC TAA TCC-3′), 0.2 μL MgCl₂ (50 mM), 1 μl BSA (25 mg/ml), 2.5 μl DMSO, 0.4 μl dNTPs (10 mM each) and 0.5 μl Phusion Polymerase (1 U/μl) per reaction. Reactions were performed in 96 well plates (Eppendorf) in triplicates using 3 μl template. Cycler settings included an activation step at 98°C for 1 min, followed by 32 cycles of denaturation at 98°C or 45 s, annealing at 65°C for 45 s, elongation at 72°C for 30 s and a final elongation step at 72°C for 5 min. Sequencing was performed at the Göttingen Genomics Laboratory (Germany) on an Illumina MiSeq System as previously described by Pohlner et al. (2017).

Bioinformatic Analysis of Sequence Data

Generated datasets were processed as follows: Timmomatic 0.36 (Bolger et al., 2014) was used to truncate low quality read ends if the average quality dropped below 15. Primer sequences were removed from amplicon sequences using bbduk¹. All datasets were subsequently processed with USEARCH v.10.0.240 (Edgar, 2010). Sequences were merged and low quality sequences were discarded (shorter than 300 base pairs (bp), accumulative sequencing error rate ≥1%) resulting in 5,336,558 high quality (HQ) amplicon sequences. HQ sequences from all samples were pooled and truncated to equal length of 400 bp before dereplication. Subsequently, dereplicated sequences were sorted by abundance. Chimeric sequences were removed and remaining sequences clustered into zero-radius operational taxonomic units (ZOTUs) using the unoise3 algorithm with a minimum unique sequence abundance of eight across all samples. ZOTUs were taxonomically classified by alignment employing the usearch algorithm against the SILVA SSU Nr99 Ref database (release 132) (Quast et al., 2013) with an e-value cut-off of 1e-10 with and minimum sequence identity of 90% and maxaccepts/maxrejects option disabled. An abundance table was created by mapping HQ sequences of each sample to the ZOTUs. The sequences were deposited at the European Nucleotide Archive under the accession number PRJEB39926.

Statistical Analysis

Statistical analyses were performed using R v3.5.3 (R Core Team, 2019) with the packages vegan v2.5-6, ape v5.3, drc v3.0-1, and picante v1.8 (Kembel et al., 2010; Ritz et al., 2015; Paradis and Schliep, 2019). Environmental parameters were normalized by z-scoring and analyzed by principal component analysis (PCA), excluding samples with incomplete environmental parameters. Relative abundances of read abundances were calculated and used in further analysis. Bray–Curtis distances of relative bacterial abundances at each station were visualized by non-metric multidimensional scaling (NMDS) (k = 2; 999 permutations). Environmental parameters were fitted to the ordination using vegans envfit-function with 9999 permutations and removal of unavailable data enabled. A permutational multivariate analysis of variance (PERMANOVA) with 9999 permutations was performed to determine the dataset variance explained by environmental parameters, excluding samples with unavailable environmental data. Core communities for sample depth, season and transect position were calculated using a category wide accumulative minimal abundance filter of 0.01%. To account for varying sequencing depth, count tables were repeatedly rarefied to 2831 sequences per sample (100 times). Subsequently, richness and Shannon entropy were calculated for each iteration and the mean value was used for further analysis. Effective Number of Species (ENoS) was calculated according to Jost (2006). ENoS and richness of sample depth, season and sampling site was compared using a Pairwise Wilcox test. Results with p ≤ 0.05 were considered significant. Also, Kruskal–Wallis rank sum test was performed to analyze differences of Bray-Curtis-Distances between within group and out group samples from August, March, and October. All comparisons were always performed between groups of samples, since no replicate cores could be recovered for the individual transects. Thus, seasonal comparisons were made with n = 40–50, station wise comparisons with n = 30–40, and different depths with n = 26 samples. Details can be found in Supplementary Table S1. Grouping was performed with the factor variables Month (October, March, August), site class (hwl, runnel, ridge, lwl) and depth (0, 10, 30, 50, 100). This approach was used to counterbalance the lack of individual sample replication.

Results

Morphodynamic Variations Govern the Hydrologic Zonation of the Sampling Sites

In all seasons, samples were taken along two transects at the same GPS coordinates and depth intervals of 0, 10, 30, 50, and 100 cm. The individual sampling sites exhibited topographical differences of up to 2.5 m. Their elevation above or below mean sea-level defined whether they were located in a seawater infiltration, or porewater exfiltration zone. Generally, cores taken at the HWL were always located in infiltration zones (Figure 2). The beach topography continuously dropped below mean sea-level between the HWL and the LWL, resulting in exfiltrating conditions for most sites. However, the ridge is considered to be an additional infiltration zone, which was especially pronounced at transect one. Topographical differences along transect two were less pronounced, resulting in a less steep beach profile. The elevation of the individual sampling sites over the three investigated seasons displayed the high morphodynamics of the intertidal zone. For instance, the elevation of the HWL of transect two dropped by almost two meters between the October and August sampling campaign. Due to intense erosion and deposition, the upper meter of sediment sampled at the same coordinates, but at different seasons, did not represent the same sediment material.

Figure 2. Relative height to mean sea-level of the individual sediment cores during the different sampling campaigns in October 2016, March 2017, and August 2017. The squares represent the samples at 0, 10, 30, 50, and 100 cm depth, respectively. Displayed are both transects spanning the whole intertidal zone, with the sites high water line (HWL), runnel, ridge and low water line (LWL) sampled always at the same GPS coordinates for each season. Site LWL b, sampled in October 2016, was in the subtidal, and therefore not accessible, in March and August.

Geochemical Porewater Composition Reflects the Heterogeneity of the Intertidal Area

A principle component analysis of the physico-chemical parameters showed no distinct clustering of sampling sites, depths or seasons, yet some trends could be observed (Figure 3). Combined, the first two principle components explain 51.92% of the variance within our dataset. Here, a clear geochemical depth zonation, as typical for sheltered marine sediments, was not evident. Focusing on site-specific differences, samples from the HWL and LWL were coupled to redox-sensitive parameters. A combined porewater analysis over all sampling campaigns showed higher concentrations of oxygenated compounds such as nitrate at the HWL, which is not as clear for the individual seasons (Supplementary Figure S1). Furthermore, salinity values were typical for seawater in these samples and pointed toward seawater infiltration. Samples from the LWL were characterized by higher concentrations of reduced compounds like dissolved iron and ammonium. Additionally, they exhibited higher concentrations of silicate, which is enriched in the island’s freshwater lens. As silicate serves as a proxy for fresh water influence (Oehler et al., 2019), it is indicative for the exfiltration of brackish porewater at the beach. The site specific orientation within the PCA is best defined for the October campaign. The samples from the infiltration zones (HWL and ridge) exhibited DOC, total dissolved nitrogen (TDN) and low temperatures, with low the HWL additionally showing an imprint of more oxygenated conditions. The LWL samples, in turn were characterized by an imprint of the SGD in terms of elevated concentrations of silicate and reduced compounds. The samples from March showed a station independent orientation toward high nitrate and oxygen concentrations, but also to low DOC, TDN and low temperature. Seasonal variations between March and August were mainly related to higher temperatures, DOC and TDN concentrations. An in depth analysis of the physico-chemical conditions at the sampling site was published previously (Waska et al., 2019b; Ahrens et al., 2020b).

Figure 3. Spatial and temporal variance of physico-chemical parameters. Porewater samples from all seasons, sites and depths were analyzed by a principle component analysis. Seasons are indicated by different symbols with varying sizes for the respective depths. The sampling sites are differentiated by color coding.

Analysis of the Microbial Diversity Reveals a Large and Resilient Core Community

After quality control, all 130 samples comprised 5,336,558 reads, resulting in 33,846 ZOTUS with an assigned taxonomy. Rare fraction analysis indicated a coverage of 85 ± 4% on species level. For further analysis, ZOTUs were summarized at species level. The diversity within the dataset was analyzed by including all species with a relative abundance of >0.01%. True diversity was calculated by determining the ENoS. Comparing the ENoS of the different samples revealed no significant correlations with depth or sampling site. In contrast, ENoS between the investigated seasons varied significantly (March–August, p = 0.017; March–October, p < 0.001; August–October, p = 0.020), while richness was only significantly different between March and August (p = 0.023), and March and October (p = 0.001).

By comparing the different sediment layers, 80% of all species occurred at all depths, indicating no depth zonation (Figure 4). This large overlap was also observed throughout the entire intertidal area with 72% individual species being detected at all sites. Among the four sampling sites, only the HWL harbored a substantial amount of individual species. Apart from the general overlap, the HWL and LWL shared the least similar communities. Over 90% of all species detected were present throughout the year. We define this group of microorganisms as a “core community,” showing no significant temporal or spatial variability. Only few community members were either exclusively present at one season or found during only two of the sampling campaigns. The overall large core community spanning all depths of the intertidal area throughout the year was probably composed of generalists, adapted to the heterogeneous and changing environmental conditions.

Figure 4. Distribution of individual microbial species within sediment samples of Spiekeroog beach. Venn diagrams show comparisons of the sampling depths (n = 26), sites (n = 30, LWL n = 40) and campaigns (n = 40, October n = 50). Diversity is based on ZOTUs with assigned taxonomy using the SILVA database.

The detailed sequence analysis on genus level displays the presence of a large core community (Figure 5). The high similarity among all sites and depths is most pronounced for both transects of the October sampling campaign. The bulk community was comprised of members of the Acidobacteria, Actinobacteria, Chloroflexi, Planctomycetes, Proteobacteria, Verrucomicrobia and a minor fraction of Bacteroidetes, Latescibacteria, and Firmicutes. Further 20–40% of the community were composed of unknowns and taxa present with relative abundances of <1%, at 90% of all stations. Those were combined and are referred to as “others.” Each of the phyla mentioned above was represented by two to five genera except the Planctomycetes and the Proteobacteria with a high internal diversity of nine genera each. Also, most of the phyla detected included a proportion of so far uncultured representatives. The community structures in the March samples revealed a significant shift in the relative abundances of the genera Gillisia as well as Planococcus and Planomicrobium, both Planococcaceae. While in October, these members of the Bacteroidetes and Firmicutes were only present in low relative abundances, their relative proportion increased at every station, yet not at every depth. Such a temporal increase in relative abundances of a few taxa is comparable to a typical bloom situation, known from pelagic environments and hereafter referred as “subsurface bloom.” However, this subsurface bloom was less pronounced in transect two. In general, deeper sediment layers were more stable, exhibiting higher similarities to the core communities identified in October. Especially at the HWL, the shift was limited to the surface sample. Additionally, an increasing relative abundance of Bacillus with sediment depth was observed at the HWL of transect one, sampled in March. Interestingly, increased relative abundances of Bacillus were also detected in August, but here at the LWL. Generally, the community structures of the August samples appeared to be a hybrid of those described for October and March. Elevated relative abundances of Planococcus and Planocmicrobium were only found at the HWL and runnel, and of Gillisia only at the HWL. At deeper depths of these stations as well as at the ridge and LWL, the microbial communities bore resemblance to the composition observed in October. Comparisons of total Bray-Curtis distances showed only small differences between August and October (Kruskal–Wallis rank sum test, chi² = 18.6 and 73.4, p < 0.001). However, both were considerably different from March samples (chi² = 1419 to 2199, p < 0.001). In general, we interpret this backshift as an indication for the reestablishment of a core community.

Figure 5. Community composition of transect one and two based on Illumina-sequencing of 16S rRNA genes. Taxonomy is based on ZOTUs, assigned using the SILVA database. The community composition is presented on genus level. Genera with an abundance of <1%, at 90% of all stations, were combined and are referred to as “others.” From left to right every block represents a sediment core of the individual sampling sites across the intertidal zone, with each bar representing one depth descending from 0 cm down to 100 cm. In October, an additional core (LWL b) was recovered from the LWL. Due to different tidal ranges this site could not be reached in later sampling campaigns.

Combined Analysis of Physico-Chemical and Microbial Diversity Data Highlights the Impact of the Beach Topography

We investigated the microbial community variation in different seasons by using an NMDS analysis and fitting environmental porewater geochemistry data to the ordination (Figure 6). Some subsets of physico-chemical data were not available for every sample and therefore excluded to avoid artificial correlations. The majority of non-HWL samples exhibit a high degree of similarity. A PERMANOVA revealed that the most influential geochemical factor is the DOC concentration, significantly explaining 11% of the difference between all samples (p ≤ 0.05) (Supplementary Figure S2). The impact of the DOC concentration on the community composition is especially high for the samples located at or slightly below sea-level. The strongest dissimilarity is related to the absolute depth above or below the sea-level, which is most pronounced for samples of the HWL. This leads to a clear separation of HWL samples from the other sites, independent of the sediment depth. The separation correlates to the response of microbial diversity to oxygenated and reduced porewater constituents, such as NO₃, NO₂, or NH₄ and Fe²⁺. Furthermore, the samples are also oriented according to salinity and silicate, which are proxies for marine and freshwater influence on porewater chemistry, respectively. The comparably weak correlations are possibly due to an overprinting of the biogeochemical depth zonation by the high physical forcing at the system.

Figure 6. The influence of porewater chemistry on microbial diversity was analyzed using NMDS analysis. Sampling sites are indicated by color coding, while the different symbols indicate the season. Grid lines indicate elevation above or below sea level. Data points closer to the center of the vector matrix are more similar to each other and influenced less by the parameters chosen.

Discussion

The aim of this research was to determine the effects of SGD on microbial diversity in the intertidal zone of a sandy high energy beach. We included seasonal variations, cross shore diversity and a wide range of geochemical parameters to obtain a holistic picture on the investigated system.

High Energy Conditions Overprint Effects of SGD on Microbial Communities

The mesotidal high-energy beach of Spiekeroog Island showed a horizontal redox zonation along the flow path of the STE. This is in accordance to a previous study by McAllister et al. (2015) which was performed on a micro-tidal beach. The authors investigated a STE at a wave-sheltered beach in Delaware over several seasons. They reported a distinct horizontal zonation of redox sensitive chemicals along the flow path. The availability and seasonal variations of chemical constituents structured the microbial community compositions. For instance, a stimulation of potentially iron oxidizing bacteria, predominantly Candidate division OP3, was observed in Fe(II) rich fresh groundwater. Sulfide was mostly produced by Desulfobacterales within the STE and transported into the intertidal mixing zone, leading to a transformation of FeOOH to Fe(II) sulfides. In contrast, the dynamics of infiltration and exfiltration zones at Spiekeroog beach governed porewater chemistry and thus overprinted a clear structuring of the microbial communities. High oxygen penetration depth at infiltration sites enhances the transport of oxidized compounds to greater depths. Physical reworking of the sediments leads to shifts in elevation, altering the groundwater flow path within the STE. Due to changes in the location of exfiltration zones, the SGD at Spiekeroog beach fluctuates along the LWL and the runnel (Waska et al., 2019b) and is much patchier than observed at low energy beaches. Additionally, the dilution of porewaters by oxygenated seawater also leads to enhanced reoxidation of the reduced compounds transported within the STE. Even though we detected areas with brackish conditions enriched in silicate and reduced components like ammonium and Fe(II), the community composition was not affected (Figure 4). Hong et al. (2018) found a vertical zonation of metabolic diversity in the upper meter of one sediment core taken from a STE in Virginia, applying molecular biological techniques. They detected a maximum of microbial diversity and abundance at the oxic-anoxic transition zone. In contrast, microorganisms present in Spiekeroog beach sediments are adapted to a broad range of changing physical and chemical conditions rather than to a specific niche within the STE.

The Formation of a Homogeneous Core Community Is Related to the Dynamics of Sandy High Energy Beaches

The strong dynamics of the system’s hydrogeochemistry force benthic microorganisms to adapt. This results in the formation of a generalist core community capable of sustaining populations under changing environmental conditions. Likewise, Boehm et al. (2014) identified a group of cosmopolitan microorganisms in supratidal sandy surface sediments from the Californian coast. They found similarities between beaches with comparable physico-chemical properties. Another comprehensive investigation compared sediment communities from eleven different beaches throughout the northern hemisphere at and above the HWL, including those of freshwater lakes (Staley and Sadowsky, 2016). The authors report that depth influences diversity less than distance to the shoreline. Both studies give a solid insight into microbial community patterns of sandy beaches above the HWL, also suggesting the presence of cosmopolitan core communities. However, studies sampling only surficial sediments at the HWL might not have captured the full diversity of the remaining intertidal area. Furthermore, seasonal changes might have been overlooked since we detected temporal disturbances of the core community composition. While the homogeneous community structure was omnipresent in October, it was still present in the deep and HWL samples in March, where the community was likely less affected by sediment reallocation and subsequent biogeochemical changes. August represented an intermediate state. Here, the communities of several depths layers were approaching a composition more similar to that observed in October, while others still showed an imprint of the bloom. Perhaps this bloom is a reoccurring event which is linked to seasonality, but this can only be confirmed when multiple years are analyzed.

Deeper Penetration of O₂ and NO₃ Lead to a Subsurface Bloom of Bacteroidetes and Firmicutes

The bloom was composed of Gillisia, a genus within the Bacteroidetes, and the two Firmicute genera Planococcus and Planomicrobium, both Planococcaceae. Gillisia have previously been identified as cosmopolitan members of beach communities along hundreds of kilometers of the Californian coast (Boehm et al., 2014). In the same study, members of the family Planococcaceae have been detected in high abundances in samples were Gillisia was also present. In general, Planococcaceae are found in a great variety of habitats ranging from soils (Luo et al., 2014) to Antarctic freshwater ponds (Reddy et al., 2002). While Planococcus has previously been found at a North Sea beach in England (Engelhardt et al., 2001), Planomicrobium was reported to be present in marine sediments in China (Dai et al., 2005). They are aerobic or facultative anaerobic with some being able to reduce nitrate (Shivaji et al., 2014).

The genera dominating the bloom might have reacted quickest to higher and deeper oxygen and nitrate availability, leading to a strong increase in their relative abundances in March. During this sampling campaign, Ahrens et al. (2020b) measured highest oxygen and nitrate concentrations within the sea- and porewater, accompanied by a deeper penetration depth of both constituents. In spring, lower temperatures resulted in higher oxygen solubility and lower respiration rates within the sediments. Further, coastal seawater nitrate was not yet exhausted by phytoplankton assimilation or potential denitrification (Grunwald et al., 2010) and thus supplied to benthic communities. The deeper penetration depth of electron acceptors at the HWL is again related to the hydrology and topography of the study site. The hydrogeological features also explain why the microbial communities of transect one showed a stronger bloom than those of transect two. Here, the steeper slope between the HWL and the LWL (Figure 2) led to stronger tidal pumping and thus to a deeper recharge of oxygen and nitrate (Ahrens et al., 2020b).

Conclusion

We could not deduce a clear correlation of the microbial community composition with STE salinity gradients along this high-energy beach. As main findings, we identified an equilibrium state of the microbial diversity all along the intertidal area down to a depth of one meter, which we defined as the core community. It has probably established due to the dynamic redox zonations caused by advection driven porewater flow and the fast physical reworking of the sediments by tide and wave action. These dynamics support the prevalence of generalist microorganisms, adapted to quickly changing environmental conditions. Over the investigated time period, the homogeneous community structure was disturbed by a subsurface bloom of individual genera, yet showing a strong resilience by returning to an equilibrium state.

Data Availability Statement

16S sequencing data have been deposited publically accessible at the European Nucleotide Archive (ENA) under the project accession number PRJEB39926. Geochemical data used for this study are available at PANGAEA under https://doi.org/10.1594/PANGAEA.905928.

Author Contributions

JD performed the research and wrote the manuscript. LD did the bioinformatics and statistics. JA and HW conducted the geochemical analysis. HW and MB revised the geochemical data and coordinated the sampling campaigns. BE was involved in sampling, writing, and revision of the manuscript. All authors read and approved the manuscript.

Funding

This study was financed by the Niedersächsisches Ministerium für Wissenschaft und Kultur (MWK) in the scope of project “BIME” (ZN3184).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We want to thank the BIME research consortium and especially Frank Meyerjürgens, Franziska Gluderer, and Manoj Mandal for their help during sampling. We would also like to thank Helmo Nikolai and Gerrit Behrends for providing logistics and equipment. Likewise, we would like to acknowledge the Niedersächsischer Landesbetrieb für Wasserwirtschaft, Küsten- und Naturschutz (NLWKN), Svante Pitters, the Umweltzentrum Wittbülten, and the Quellerdünen hostel, for on-site logistics and infrastructure. We are especially grateful to Verona Vandieken, for her assistance during sampling and valuable scientific discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.573570/full#supplementary-material

Supplementary Figure 1 | Values of oxygen, nitrate, ammonium, iron, and manganese for the infiltration zone (HWL and Ridge) and the exfiltration zone (LWL and Runnel). Sampling depths were combined to “shallow” (0–30 cm) and “deep” (50–100 cm). (A) combined over all sampling campaigns, (B) October 2016, (C) March 2017, (D) August 2017. For the October campaign (B), oxygen values were not available.

Supplementary Figure 2 | PERMANOVA of Bray-Curtis dissimilarities and physicochemical parameters; n indicates the number of samples for the respective parameter.

Supplementary Table 1 | Combined physico-chemical and diversity data used for this study.

Footnotes

1. ^ https://sourceforge.net/projects/bbmap/

References