The study site is located at the mouth of the York River Estuary at Gloucester Point (37.248884 N, 76.505324 W), VA, within the Chesapeake Bay. The details of the site description were previously reported by Beck et al. (2016). SGD rates ranged from 3.9 to 8.9 cm day^-1 depending on discharging locations (Luek and Beck, 2014; Beck et al., 2016). The GP-STE contains fine grained and permeable sediments.

A sediment core, about 1.0 m in length, was collected using a vibrocorer in October of 2015. Sediments in a single core were sub-sectioned at an about 10 cm interval, for a total of 10 subsamples (Figure 1, Table 1, and Supplementary Figure S1) were obtained and stored at -80°C until further analysis.

Groundwater at different depths of the study site was sampled in parallel with the sediment core sampling with a drive-point piezometer system called Retract-A-Tip (AMS, Inc.; Charette and Allen, 2006). The stainless-steel piezometer was driven to the depth of interest, and groundwater samples were pumped through Teflon tubing using a peristaltic pump. Samples for nutrients were collected into 30 mL acid-cleaned scintillation vials after passage through a Pall Aquaprep 0.22 μm capsule filter and stored frozen until analysis. Water properties including temperature, salinity, pH, dissolved oxygen, and Oxidation-Reduction Potential (ORP) were recorded using an YSI 600XLM multiprobe and 650MDS handheld computer. Iron (Fe) was measured using the calorimetric Ferrozine method (Viollier et al., 2000) which fixes Fe in the sample for measurement of Fe²⁺ prior to chemical oxidation for measurement of Fe³⁺ to determine total iron. Dissolved N₂O concentration was measured by collecting 30 mL of groundwater in 50 mL serum bottles with a pellet of KOH. After vigorous shaking to equilibrate dissolved gases, 10 mL of headspace gas was sampled and injected to a gas chromatograph equipped with electron capture detector (Shimadzu, Kyoto, Japan).

DNA was extracted from 0.5 g of sediment (wet weight) using a Power soil^®DNA isolation Kit (Mo-Bio Laboratories, Inc., Carlsbad, CA, United States) according to the manufacturer’s instructions. Thermo Savant Fast Prep FP 120 Cell Disrupter (Qbiogene Inc., Carlsbad, CA, United States) was used for cell disruption. PCR was conducted to amplify the V₄ hypervariable regions of 16S rRNA genes of bacteria and archaea using forward 515F and reverse 806R primers (Caporaso et al., 2011). The reverse primer was modified with 8 bp barcodes, and both forward and reverse primers were added with adapter primers corresponding to the sequencing protocol of the Ion Torrent Personal Genome Machine (PGM). Aliquots of a 25 μL PCR reaction included 12.5 μL 2 × Go-Taq master mix, 1 μL of each primer (5 μM), 1 μL template DNA and 9.5 μL nuclease-free water. The PCR reactions were performed in duplicate under the following conditions: initial denaturation at 95°C for 2 min, 25 cycles of 95°C for 30 s, 55°C for 30 s, 72°C for 30 s, and then a final extension at 72°C for 7 min. A negative control without DNA template was included in each PCR reaction to detect for any contamination of genomic DNA. The amplicons of each sample were pooled and purified with the Wizard^® SV Gel and PCR Clean-Up System.

The purified PCR products were quantified using a 2200 Tape Station instrument and D1K reagents (Agilent Technologies, Santa Clara, CA, United States) following the manufacturer’s instructions. Sequencing was performed using an Ion Torrent PGM sequencer with the pooled samples on a 316 chip, following the manufacturer instructions (Life Technologies, Grand Island, NY, United States). The sequence reads with more than one ambiguous nucleotide (N) and incomplete sequences of barcode and 16S primers were removed from a subsequent bioinformatic analysis. The filtered sequences were assigned to each sample based on the barcode.

All sequence reads with complete match of the barcodes and a single mismatch to the 16S primers were retained and then trimmed by removing the sequencing adaptor, barcodes and 16S primer sequences to obtain valid raw reads. The reads were further screened by using the following thresholds: (i) minimum average quality score of 25; (ii) minimum read length of 200 bp; (iii) sequences containing no ambiguous bases; and (iv) maximum homopolymers of 8 bp. Quality-controlled sequences were analyzed in MOTHUR (Version 1.35.1) following standard operating procedures² (Schloss et al., 2009). The Greengene database (gg-13-5-99.align) was used to align and classify the reads. The OTUs containing more than 1% of total sequences were defined as dominant OTUs. Phylogenetic analysis of the dominant OTUs (Top 50 OTUs) was conducted with the software MEGA7.0 and Evolview (He et al., 2016). A heat map was constructed in accordance with the abundance of dominant OTUs using Microsoft Excel.

Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt)³, a predictive exploratory tool, combined with the Kyoto Encyclopedia of Genes and Genomes (KEGG) ortholog classification, was used to infer the metabolisms of sediment communities and to generate predicted metagenomes based on the classified 16S rRNA reads of each sample (Langille et al., 2013; Kanehisa et al., 2014). The functional genes in methane, nitrogen and sulfur metabolisms were selected to compare potential functions of microbial communities in the sediments of GP-STE.

The q-PCR assays of the 16S rRNA gene were performed to quantify the abundance of bacteria and archaea in the sediment communities. The q-PCR reactions were performed in triplicate within a volume of 20 μL containing SYBR green using Go-Taq qPCR Master Mix (Promega Corporation, Madison, WI, United States) with the primer set of 515-F and 806-R using ABI Prism 7500 Real Time PCR System (Applied Biosystems, Carlsbad, CA, United States). Each reaction was performed in a 25 μL volume containing 2 μL of DNA template, 0.2 μL BSA (0.1%), 1 μl of each primer (20 μM, 515F, and 806R), 12.5 μL of Power SYBR Green PCR Master Mix and 9 μL of nuclease free water (Applied Biosystems, Foster City, CA, United States). The q-PCR cycle was as following: 10 min at 95°C, followed by total of 35 cycles of 95^oC for 30 s, 55^oC for 30 s, 72^oC for 30 s. A standard plasmid carrying the 16S rRNA gene was generated by amplifying 16S rRNA gene from DNA extracted from GPMD sediment and cloned into the pMD-T-18 Vector (Takara, Japan). The plasmid DNA concentration was determined on a Tape Station 2200 (Agilent) and the copy numbers of target genes were calculated directly from the concentration of the extracted plasmid DNA. Ten-fold serial dilutions of a known copy number of the plasmid DNA were subjected to q-PCR assay in triplicate to generate an external standard curve.

Chao and Shannon estimators were generated in Mothur as proxies of the alpha diversity index. Redundancy analysis (RDA) was performed to examine covariance among environmental variables and the dominant OTUs using the Canoco 5.0 software. Principal coordinate analysis (PCoA) was performed in Mothur software on relative sequence abundance data at the OTU level (Schloss et al., 2009). The vertical variation in environmental parameters and microbial abundance were analyzed using one-way analysis of variance (ANOVA), and significance level was set at α < 0.05.

The reads of the 16S rRNA genes were deposited in the ENA Short Read Archive under the submission number SRP154289.

Depth profiles of environmental characteristics [salinity, temperature, dissolved oxygen (DO), pH, N₂O, Fe (II), Fe (III), Fe (total), and ORP] and microbial abundance were shown in Figure 1. We found sharp gradients of salinity, DO, and ORP in the groundwater sampled from different depths of the sampling site in the GP-STE. The salinity increased slightly from 17.80 psu at 30 cm to 18.20 psu at 40 cm, and then declined sharply to 1.26 psu at depth of 90 cm. The spatial distribution of salinity indicated that seawater can be transported to the STE via two different ways: vertical permeation from the surface, and horizontal transportation from seawater terminal. The DO profile was measured using an oxygen probe, revealing a sharply decline from 7.3 ppm at 30 cm to 0.14 ppm at a depth of 90 cm. The temperature gradually increased along the profile from 12.70 to 17.00°C. pH value declined sharply from 7.20 at 30 cm to 6.60 at 60 cm, and then increased from 6.60 to 7.20 at 80 cm. The N₂O concentration varied from 0.06 to 0.13 ppm and peaked at the depth of 60 cm. In addition, the N₂O concentration was under detection level below depth of 60 cm, Similar to the profile of DO, the ORP value decreased with sediment depth from 157.40 to -127.00 mv. The dissolved Fe peak occurred at the depth of 70 cm where the sand layer was highly enriched with iron oxides. The profile of microbial abundance showed no significant variation (P > 0.05) in sediment core [ranging from (0.23 ± 0.03) × 10⁹ to (4.62 ± 0.08) × 10⁹ copies g^-1]. There was a remarkably high microbial abundance in the aerobic-anaerobic transition zone (AATZ, depth of 50–60 cm) than other depth sediments. The microbial abundance was lowest in the anaerobic zone (GP 8–10, depth of 80–100 cm).

The Ion Torrent PGM generated about 350,000 raw reads of 16S rRNA genes from 10 sediment samples. After removing of noise and low-quality reads, more than 270,000 reads were used for further analysis (the high quality reads in each sample were listed in Table 1). In total, 31,146 OTUs were assigned at a 97% sequences similarity threshold with the most abundant OTUs (>100 reads per OTU) accounting for 49% of all sequences. There were 30,263 singletons accounting for 52.9% of the total OTUs that might have resulted in the high microbial diversity. The number of OTUs ranged from 3,405 to 6,158 across all samples, with sample GP7 (at depth of 70 cm) sediments harboring significantly higher numbers of OTUs than other samples in the sediment core. The Shannon diversity index was also significantly higher in GP7 sediment community than others (P = 0.004). The coverage of each sample ranged from 86.0 to 96.0%, with an average value of 89.8%, indicating that the sampling efforts of all samples were sufficient to represent indigenous species. Rarefaction curves showed that the high-throughput sequencing provided sufficient bioinformation to investigate the community composition and diversity of bacteria and archaea in this study (Figure 2A). Similarly, ACE and Chao1 estimators showed clear differences in species richness across the sampling depths, which were about twofold higher than the OTU numbers (Table 1), suggesting that twice as many OTUs may exist in the sediments.

At a confidence threshold of 80%, 273,103 out of 273,789 qualified reads (count for 99.7%) could be assigned to a known phylum using the Mothur classification program, of which 59,128 belonged to archaea and 213,975 were bacteria. Relative abundance of bacteria decreased from 88.2 to 74.6%, while archaeal abundance increased correspondingly and ranged from 11.8 to 25.4% as the depth increased (Figure 2B). Four archaeal phyla, Crenarchaeota, Euryarchaeota, Parvarchaeota, and an unclassified phylum, were found in the sediment samples with varying relative abundance at different depths. Crenarchaeota was the dominant species in the archaeal communities, constituting up to 60% of the archaeal reads, followed by Euryarchaeota (30%), which were retrieved mainly from bottom sediments (80–100 cm). Parvarchaeota was a minor group, and the relative abundances in the upper samples were higher than those in the lower samples (Figure 2C).

In total, 70 bacterial phyla were recovered from our 9 samples (Supplementary Table S1). When comparing the bacterial communities in different phyla, Proteobacteria (30.70%), Chloroflexi (30.01%), Planctomycetes (8.93%), Bacteroidetes (3.12%), Acidobacteria (1.97%), Cyanobacteria (2.20%), OD1(1.85%) and an unclassified phylum (6.63%), altogether constitute up to 86.6% of the reads affiliated with bacteria (Figure 2D). However, reads belonging to Gemmatimonadetes, WS3, Chlorobi, Nitrospirae, Chlamydiae, Verrucomicrobia, TM6, OP8, Caldithrix, OP3, Spirochaetes, Firmicutes, GN02, and GN04 were found to be the minor groups. The heterogeneous distributions of some phyla along the depth profile in the sediment core were observed. For example, Euryarchaeota, Chloroflexi, and WS3 displayed an increasing trend with depth. In contrast, Proteobacteria, Bacteroidetes, Cyanobacteria and Acidobacteria showed a decreasing trend with depth.

Based on the OTU analysis using PCoA calculations, both bacterial and archaeal community structures displayed variations across the depth profile of the sediment core (Figure 2). The microbial communities at the top to 30 cm (GP1, GP2, and GP3) tightly clustered together and were different from those below 30 cm depth.

Figure 3 shows the phylogenic tree of 50 dominant OTUs and their distribution characteristic across all samples. Among the bacteria OTUs, Gammaproteobacteria (12 OTUs), Acidibacteria (OTU26), Betaproteobacteria (OTU47, OTU32, OTU48), Acidimicrobiia (OTU34), Chloroplast (OTU19 and OTU15), Flavobacteriia (OTU23), Mb-NB09 (OTU46), and Nitrospira (OTU29) were dominant in the top to 70 cm sediment layers. In contrast, Deltaproteobacteria (OTU17, OTU7, and OTU37), Dehalococciodetes (OTU17, OTU7, OTU37, OTU41, and OTU39) and Phycisphaerae (OTU39) were dominant members in the 80–90 cm depth sediment layers. Among the archaeal OTUs, Thaumarchaeota (8 OTUs) was mainly distributed in the top to 70 cm sediment layers while MCG Methannobacteria (OTU4) and MBGB (OTU6) were only found in the 80–90 cm depth sediment layers. Moreover, there was a higher variety of dominant OTUs in 49–70 cm depth of sediment layers (GP 6–8) where all of the 50 dominant OTUs were detected. The relative higher diversity in the zone was mirrored by the rarefaction analysis in Figure 2A.

Obvious vertical distribution of community composition was revealed by unweighted unifrac PCoA analysis (Figure 4A). PCoA results indicated that bottom and anoxic community assemblages (GP 8–10) were separated from others. Redundancy analysis (RDA) was used to evaluate the relationship between the community composition and the physicochemical characteristics of sediments (Figure 4B). Results showed that sampling depth, DO, ORP and salinity were key physicochemical characteristics affecting the distribution of bacteria and archaea in sediments in the GP-STE.

The metabolic functions based on the classified OTUs were inferred using the PICRUSt program, which categorizes functional gene families in nitrogen, methane, and sulfur metabolism (Table 2). The annotations were classified into 330 functional KEGG or orthology pathways at level 3 (Supplementary Table S2), which included any kind of biological reaction and regulation of gene expression.

The genes for aerobic metabolism, including those encoding ammonia monooxygenase (amoA genes), methane monooxygenase (mmoA gene), and sulfite oxidase, were predicted to have the highest abundance in the oxic surface to 70 cm depth, but the lowest in the deeper depths below 70 cm (GP 8–10). In contrast, predicted abundance of nitrous oxide reductase (nosZ) genes, nitrogen fixation (nifH) genes, methanegen (mcrA) genes, and sulfite reduction (dsrA) genes showed an opposite trend to the aerobic metabolism genes. Results also predicted that abundance of nirK and nrfA gene was higher in AATZ (50–60 cm).

In this study, the microbial community compositions of sediments at different depths from surface to 100 cm of permeable sediments were examined using high-throughput sequencing of the 16S rRNA gene. Our results revealed a distinct pattern of vertical distribution consisting of aerobic, aerobic-anaerobic transition (AATZ) and anaerobic zones. The composition of the microbial community in different zone was distinctly different. The high indexes of OTU and Shannon diversity (Table 1) in sediment samples confirmed that the permeable sediments harbor higher microbial diversity similar to estuarine sediments (Liu et al., 2014a,b). Moreover, apparent depth-related differences in measurement of diversity were found in the sediment core, which contrasts with previous studies in estuarine sediments. Furthermore, the microbial diversity index appeared to be significantly higher in the AATZ than in the shallow and deep layers. The variation along the depth profile in the sediment core with seen with PCoA analysis suggests that permeable sediments undergo predictable changes from the surface to bottom layers, which may be significantly different from the distribution of the microbial community structure in less permeable sediments. In the sediment core, bacteria were numerically dominant relative to archaea, with 2–3 orders of magnitude, while the relative abundance of archaea gradually increases from 4.4 to 22% with increasing depth (Figure 2B), suggesting archaea may play more important roles for the biogeochemical processes in the deeper zone than in the upper zone. Moreover, Crenarchaea accounted for more than 50% of the archaeal reads, suggesting that Crenarchaea may play more important roles than other archaeal groups. A similar distribution was also observed in the estuarine sediments (Liu et al., 2014a,b) and mangrove sediments (Bouchez et al., 2013).

It has been suggested that most microorganisms are restricted to specific niches, mainly due to different environmental forces (Lee et al., 2011; Qian et al., 2011; Nelson et al., 2016; Zhang et al., 2017). Our results provide plausible explanations for vertical variations of microbial communities in the permeable sediment. The temporal shift in microbial composition appeared to be associated with the vertical stratification of physiochemical parameters, including temperature, salinity, DO, corresponding significantly to changes in the microbial communities. Previous studies have shown that DO was one of the most important factors affecting the distribution of microbial communities in estuarine waters (Kan et al., 2006), marine subsurface sediments (Durbin and Teske, 2012) and in flooded Paddy Soil (Ludemann et al., 2000). Results of the present study also highlighted the importance of DO in governing the distribution of both bacterial and archaeal communities in the sandy permeable sediments. The DO profile exhibited a sharp decline from 7.3 ppm at 30 cm to 0.14 ppm at a depth of 90 cm and the microbial community distribution was closely related to the DO profile. Above 50 cm, microorganisms were mainly aerobic and facultative aerobic groups, including Gammaproteobacteria, Beltaproteobacteria, Acidobacteria, Chloroplast, Flavobacteria, Thaumarchaeata, and Mb-NB309. At a depth of AATZ (50–60 cm), anaerobic groups appeared gradually increased, showing a co-existence with facultative groups in the transition zone with the highest diversity. Below 70 cm depth, anaerobic microbial groups were dominant in the sediments, including methanogenic archaea, sulfate reducing bacteria, and reductive dehalogenating bacteria.

Previous studies have suggested that salinity is an important factor in shaping the compositions of bacteria communities in estuaries (Crump et al., 1999; Bouvier and del Giorgio, 2002). However, salinity varied with many environmental parameters across the sediment core, including DO, pH, and temperature; the role of salinity has not been obviously discriminated for the distribution of microbial community structures. Furthermore, the salinity in the sandy sediments is in a state of fluctuation due to the influence of tides (Moore, 1999; Robinson et al., 2007). Therefore, the microbial community in the sediments should be periodically impacted by different salinity. The effects of increased salinity on N-cycling microbes and their associated geochemical functions have been previously reviewed (Santoro, 2010). Salinity appeared to select for a less diverse microbial community in most cases. The activities of different functionally microbial groups have different responses to the change of salinity, for example, the rate of DNRA increased when the salinity increased to 10 psu, but the denitrification process was depressed with the same treatment (Laverman et al., 2007). Although our data showed that the salinity was an important factor that influenced the distribution of microbial groups, no clear pattern was observed for how defined groups of bacterial or archaea community shift with the salinity variation. Further study of the salinity effect of on composition and function of microbial communities will increase our understanding of microbial processes in STEs.

The surface of the core is rich in oxygen, which provides conditions for the aerobic metabolism of microorganisms. The predicted metagenome analysis showed microbial aerobic metabolism, including methane oxidization (Bender and Conrad, 1994), ammonia oxidization (Ward and Jensen, 2014), and sulfite oxidation (D’Errico et al., 2006), mainly occurring in the aerobic zone. Similar to other STE locations (Charette and Sholkovitz, 2002; Roy et al., 2013), a bright red-orange layer of precipitated iron oxide was visible in the sand cores from 55 to 70 cm. The intensity of this oxide layer and the levels of dissolved Fe appeared to fluctuate in both time and depth, which we believe was linked to hydrologic controls, similar to other locations (Michael et al., 2005; Gonneea et al., 2013). In the AATZ, almost all dominantly predicted genes related to denitrification, DNRA and sulfite reduction were observed with relatively high abundance, suggesting both aerobic and anaerobic metabolism could occur in the transition layer. Therefore, the aerobic-anaerobic interface should be a hotspot for microbial functional metabolism. Microbial anaerobic metabolic pathways were primarily found in the anaerobic zone, for example sulfate reduction and metagenesis. Overall, the predicted microbial functions also presented a typical stratification feature, and function prediction showed that the variation of potential microbial metabolic function was matched to the community shift in the sediment core, showing a unification of the distribution of microbial diversity and functional metabolism.

Oxygen plays an important role in regulating the function of microorganisms. In fact, most of the microbes were of the facultative respiration type, capable of using a variety of electron acceptors (Hogg, 2005). Under low oxygen conditions, a variety of electron acceptors can be used alternately to improve the efficiency of metabolic productivity (Jose, 2004; Canfield et al., 2005). In particular, the concentration of dissolved oxygen will fluctuate slightly, especially in the STEs where fresh ground water and sea water are mixed (Moore, 1999; Gonneea and Charette, 2014). Under these conditions, the metabolism of microorganisms will also be significantly regulated by dissolved oxygen in the pore water. Many microorganisms can alternate or simultaneously carry out aerobic and anaerobic respiration under low oxygen conditions. The higher diversity of both OTUs and predicted functions demonstrated that aerobic and anaerobic metabolism existed concurrently in the AATZ. Combined with the diversity and the predicted functions of the microbial communities, we propose that the aerobic-anaerobic interface in the highly permeable of sandy sediments is a hotspot of microbial geochemical reactions, regulating nutrient discharges to surface water.