Bacterial Diversity in Submarine Groundwater along the Coasts of the Yellow Sea

Submarine groundwater (SGD) is one of the most significant pathways for the exchange of groundwater and/or source of nutrients, metals and carbon to the ocean, subsequently cause deleterious impacts on the coastal ecosystems. Microorganisms have been recognized as the important participators in the biogeochemical processes in the SGD. In this study, by utilizing 16S rRNA-based Illumina Miseq sequencing technology, we investigated bacterial diversity and distribution in both fresh well water and brackish recirculated porewater along the coasts in the Yellow Sea. The results showed that Actinobacteria and Betaproteobacteria, especially Comamonas spp. and Limnohabitans spp. were dominated in fresh well samples. Distinct patterns of bacterial communities were found among the porewater samples due to different locations, for examples, Cyanbacteria was the most abundant in the porewater samples far from the algal bloomed areas. The analysis of correlation between representative bacterial taxonomic groups and the contexture environmental parameters showed that fresh well water and brackish porewater might provide different nutrients to the coastal waters. Potential key bacterial groups such as Comamonas spp. may be excellent candidates for the bioremediation of the natural pollutants in the SGD. Our comprehensive understanding of bacterial diversity in the SGD along the coasts of the Yellow Sea will create a basis for designing the effective clean-up approach in-situ, and provide valuable information for the coastal management.


INTRODUCTION
Submarine groundwater discharge (SGD) includes both fresh meteoric groundwater to the ocean by terrestrially-driven directly and recirculated seawater by permeable sediment (Garcia-Solsona et al., 2010). More and more works showed that SGD is one of the most important pathways for the exchange of groundwater and/or source of dissolved compounds (e.g., nutrients, metals, carbon) to the ocean, which may cause negative impacts on the coastal ecosystems (Burnett et al., 2006;Moore, 2010). Nitrogen input through the SGD may contribute to the eutrophication of coastal regions (Moore, 2010). Several case studies showed that increased nutrient supply via SGD might be a key factor for initiating and fueling the persistent harmful algal blooms (HABs) (LaRoche et al., 1997;Hu et al., 2006;Lee and Kim, 2007;Smith and Swarzenski, 2012).
The Yellow Sea, which is as a shallow semi-enclosed water body of the Northwest Pacific Ocean, shows such a strong seasonality in the nutrients concentration and their components . SGD has been determined as a main nutrient source to the Yellow Sea, which is one of the largest continental shelves in the world Waska and Kim, 2011). Microorganisms have been shown as the important components mediating chemical reactions in the SGD in some concerned cases Paytan et al., 2004;Halliday and Gast, 2010). Recent studies focused on the bacterial communities in the surrounding sea water during macroalgal blooms of Ulva prolifera in the Yellow Sea (Burke et al., 2011;Liu et al., 2011;Zhang et al., 2014). For example, Analysis of nifH gene clone library revealed that heterotrophic nitrogen fixers, mostly represented by Vibrio-related Gammaproteobacteria, were dominated in both of surface waters that were covered and noncovered with massive macroalgal canopies of U. prolifera in the Yellow Sea in the summer of 2011 (Zhang et al., 2014). However, for our knowledge, up to now, there are few reports on the bacterial compositions in SGD along the coast of the Yellow Sea.
Microorganisms play pivotal roles in biogeochemical cycles in marine ecosystems (DeLong and Karl, 2005;Zehr and Kudela, 2011). It was reported that several studies using 16S rRNA (Fields et al., 2005), functional genes (Yan et al., 2003;Santoro et al., 2006Santoro et al., , 2008, microarray (Waldron et al., 2009), and metagenomics approaches (Hemme et al., 2010) to reveal bacterial communities in groundwater systems. The decreased diversity was observed in groundwater contaminated with higher levels of nitric acid and uranium waste in the Field Research Center (FRC) of the U.S. Department of Energy Environmental Remediation Science Program (Oak Ridge, TN, USA) (Fields et al., 2005). Ammonia-oxidizing bacterial and archaeal abundance as well as denitrifier communities were characterized along a nitrate and salinity gradient in groundwater samples at Huntington Beach, CA, USA (Santoro et al., 2006(Santoro et al., , 2008. Metagenomical analysis showed that a groundwater microbial community composed of clonal denitrifying γ -and β-proteobacterial populations in an extreme low-pH environment contaminated with high levels of uranium, nitric acid, technetium and organic solvents (Hemme et al., 2010). Biodegradation or biotransformation of pollutants with microorganisms is the most potential cost-effective bioremediation technology in the natural environments (Fields et al., 2005;Hwang et al., 2009). Pilot-scale uranium in situ bioremediation experiment suggested that the indigenous microbial communities can be successfully stimulated after 2 years of biostimulation with ethanol in FRC groundwater ecosystem, U(VI) levels were reduced from the initial concentration of 50 mg/L to below drinking water standard (<30 µg/L) (Xu et al., 2010). Zhang et al. (2015) reported that amendment of the slow-release polylactate hydrogen-release compound (HRC) stimulated FRC groundwater microbial communities associated with HRC degradation and reduction of NO − 3 , Cr(VI), Fe(III), and SO 2− 4 . In the present work, by using 16S rRNA gene-based Miseq Illumina sequencing approach, we focused on the bacterial compositions in both fresh well water and recirculated brackish porewater along the coast of the Yellow Sea. The goal is to explore the key bacterial groups, which can play their potential ecological roles in nitrogen and carbon flux by SGD to the coastal areas. The relationship between representative taxonomic groups and contexture environmental parameters were also investigated. Our understanding on the potential bacterial candidates for bioremediation will help us design the effective clean-up technology for coastal environmental management.

Sample Collection and Measurements of Physic-Chemical Parameters
Groundwater samples (well and pore water) in our field observation were collected from wells and beach in May, 2014. Well samples along the shore were pumped from the wells utilizing the water pump; pore water samples from offshore zone were collected from a push-point piezometer using a peristaltic pump (Charette and Allen, 2006). The temperature, salinity and pH of groundwater were measured directly in the field using a portable salinometer with multiple parameters (Germany, multi 350i). Sampling depth of each sample site ranged from 0.4 m to > 20 m ( Table 1).
Four well water samples (∼40 L) and three pore water samples (∼20 L) were collected for Ra isotope measurement. After filtrating (pore size: 0.5 µm), Ra isotopes in the water were extracted using a MnO 2 -impregnated acrylic fiber column (20 g), then the column with Mn-fiber was immediately placed in the Radium Delayed Coincidence Counter (RaDeCC) to measure the short-lived isotope 223 Ra and 224 Ra in the field (Moore and Arnold, 1996).
Corresponding water samples (∼60 mL) for nutrient analysis were collected with polyethylene bottles filtered through 0.45 µm cellulose acetate filters. Then the filtrates were poisoned with saturated HgCl 2 and stored in the dark. The nutrient concentrations (NO − 2 , NO − 3 , NH + 4 , PO 3− 4 , and SiO 2− 3 ) were then analyzed using an auto-analyzer (Model: Skalar SANplus146) (Liu et al., 2005). DOC samples were filtered via clean Nylon filter (pore size: 0.45 µM) immediately after collection and kept at −20 • C until analysis. In the lab, DOC samples were measured with a TOC analyzer (Shimadzu R TOC-LCPH).
Replicate, 300-mL of each submarine groundwater sample were collected on a 0.22 µm pore size polycarbonate filter (Nuclepore Track-Etched Membrane, Whatman). The filter was placed in a sterile 1.5-mL microcentrifuge tube and immediately store at −20 • C.

DNA Extraction, PCR, and Sequencing
Total DNA was extracted from the filter using a MoBio PowerWater R DNA Isolation Kit (MOBIO Laboratories, Carlsbad, CA, USA) according to the manufacturer's instruction. DNA from two independent extractions were combined, and their concentration and purity were measured spectrophotometrically with NanoDrop ND2000.
Minimum numbers of PCR cycles were performed and three independent PCR mixtures were pooled for each sample to decrease PCR bias. Briefly, bacterial V4-V5 hypervariable regions of 16S rRNA genes were then amplified using the specific barcoded universal primer pairs 515F (5 ′ -GTGCCAGCMGCCGCGG-3 ′ ) and 907R (5 ′ -CCGTCAATTCMTTTRAGTTT-3 ′ ) (Xiong et al., 2012). Cycling conditions were 2 min at 95 • C, followed by 25 cycles with 30 s at 95 • C, 30 s at 55 • C, and 30 s at 72 • C, and a final extension period of 5 min at 72 • C. PCR products were purified using the AxyPreDNA gel extraction kit (Axygen Biosciences, USA) following the manufacturer's protocol and then quantified by QuantiFluorTM-ST (Promega, USA). Reaction mixtures were pooled in equimolar ratios and paired-end reads were generated on an Ilumina Miseq PE250 (Majorbio Bio-Pharm Technology Co., Ltd., Shanghai, China).

Sequence Data Process, OTU Cluster, and Taxonomic Assignment
Raw fastq files were demultiplexed, quality-filtered using QIIME (version 1.17) (Caporaso et al., 2010) with the criteria as described previously (Li et al., 2014). Read data of each sample was imported in Fastq format according to index barcode sequence. Operational Taxonomic Units (OTUs) were clustered at 97% similarity using UPARSE (version 7.1 http://drive5.com/ uparse/). Chimeric sequences were identified and removed using UCHIME. By using the usearch_global command, the number of reads from each sample assigned to each OTU was generated in an "OTU table." Representative 16S rRNA gene sequence of each OTU was analyzed by RDP Classifier (http://rdp.cme.msu. edu/) against the Silva (SSU115) 16S rRNA database with 70% confidence threshold.

Phylogenetic Analyses
The sequences of the representative OTUs obtained in this study were compared to those in the National Center for Biotechnology Information nucleotide database by using BLAST searching. The closest sequences and selected reference sequences were downloaded and aligned using Clustal W. Phylogenetic tree was created in MEGA6 using the neighbor-joining algorithm with a bootstrap test of 1000 replicates and maximum composite likelihood model (Tamura et al., 2013).

Statistical Analyses
Alpha diversity metrics and coverage were calculated using the Mothur Program (Schloss et al., 2009). R program was used to plot Heatmap figure (R Development Core Team., 2013). Bray-Curtis dissimilarities were calculated based on genus/sample abundance matrix using function vegdist of the vegan package in R (Oksanen et al., 2007). Hierarchical clustering of the samples was performed with the complete method using the function hclust of stats package in R (R Development Core Team., 2013).

Nucleotide Sequence Accession Numbers
Paired end Illumina sequence data from this study were submitted to the NCBI Sequence Read Archive (SRA) under accession number SRP064231.

Site Description and Physic-chemical Characteristics
Fresh well water and brackish porewater of submarine groundwater samples were selected in this study. Wells along the shore were dug by native habitants; Porewater within Yellow Sea offshore area were collected from ∼50 to 130 cm below the sediment-water interface ( Table 1).
Based on the types of submarine groundwater and physicochemical parameters, seven samples were separated into two groups: Group W included four samples (YSGW1, YSGW3, YSGW4, and YSGW11) collected from the wells, the salinity of well waters ranged from 0.3 to 1.4 , indicating their freshwater environment; Group P included three porewater samples (YSPW2, YSPW7, and YSPW11) with relative higher salinity (∼30 ) (Table 1, Figure 1).
The other hydrological properties and chemical parameters of submarine groundwater are also summarized in Table 1. Fresh well waters had an average value of 0.269 and 0.0083 dpm/L for 224 Ra and 223 Ra. However, pore water had a higher value for 224 Ra and 223 Ra, was 6.805 and 0.144 dpm/L, respectively. Several nutrients in the specific well water samples were much higher than those in the pore water, for examples, nitrate concentrations ranged from 243.9 to 540.1 µmol/L within Group W except sample YSGW 4 (20.4 µmol/L), but still higher than those in porewater samples (1.6-14.8 µmol/L). Ammonium concentration ranged 8.09-19.40 µmol/L in porewater samples, and 0.98-1.55 µmol/L within Group W except sample YSGW 4 (9.56 µmol/L). Porewater sample YSGW 2 had the highest nitrite concentration and well water sample YSGW 1 had the highest phosphate concentration. The silicate concentrations in porewater samples were below 15 µmol/L, while the concentrations in well samples were much higher (116.2-333.3 µmol/L). No distinct pattern of DOC concentrations between well water and porewater samples, but YSGW 4 had the highest DOC concentration. The other parameters such as pH and temperature were relatively constant at all sites.

Bacterial Diversity and Cluster Analysis
We obtained a total of 87,475 high quality bacterial V4-V5 Illumina sequences, with an average read length of 396 bp from seven submarine groundwater samples.
According to 97% identity cutoff, there were 1078 OTUs in the complete data set. Among them, only 23 OTUs were found in all well water samples and 169 OTUs in porewater samples. Good's coverage was 99∼99.8% for all samples, representing almost whole range of bacterial diversity. The Chao, ACE value, and Shannon evenness all indicated that the α-diversity of bacterial community is lower in the well water samples than in porewater samples ( Table 2).
The hierarchical heatmap at the bacterial genus level showed that seven samples could be organized into two main groups: The first group was composed of four well water samples, YSGW1 and YSGW4 were clustered together, and YSGW3 and YSGW11 formed another cluster; the second group was composed three porewater samples, turbid porewater YSPW2 was clustered separately from clear porewater YSPW7 and YSPW11 (Figure 2). Figure 3 showed that there were some variations in the percentage composition of illumina sequences among the samples. Proteobacteria was the most abundant phylum in all samples (43.7∼76.7%) except YSPW11 (25.9%). Bacteriodetes, which was the most abundant phylum in YSPW11 (27.5%), was also distributed in the other samples. Actinobacteria was the other prevalent phylum in wellwater samples YSGW3 (27.2%), YSGW4 (16.7%), and YSGW11 (11.6%). Planctomycetes and Acidobacteria formed the second (17.2%) and third (12%) dominant groups after Proteobacteria within YSPW2. Cyanobacteria represented 32.8 and 22.7% of sequences within porewater samples YSGW7 and YSGW11. Interestingly, Cyanobacteria was almost absent from the other five samples.

Phylogenetic Analyses
Main Groups of Actinobacteria and Betaproteobacteria within Group W Samples OTUs related to Actinobacteria and Bataproteobacteria were frequently detected within fresh well water samples (Figure 4).

Main Groups of Bacteriodetes Detected in Both Group W and P Samples
Bacteriodetes was frequently observed in all studied submarine groundwater samples, however, the main groups of Bacteroidetes were different between Group W and P samples (Figure 6). Five Flavabacterium related OTUs were found in YSGW1 and YSGW4, indicating microbial degradation of high-molecularweight organic matter in these samples (Kirchman, 2002). Two OTUs 19 and 275 in YSGW4 matched uncultured environmental clones from Lake Zurich (HE574353 and HE574358), which may participate in N-acetyl-glucosamine uptake (Eckert et al., 2012). OTU 219 had a significantly higher abundance in the YSGW4 and 100% identity to uncultured Flectobacillus clone ZS-2-379 (FN668111) (Van den Wyngaert et al., 2011), members of the grazing-resistant Flectobacillus genus had the ability to defense filament-forming predators (Šimek et al., 2010a). OTUs 159 and FIGURE 2 | Heatmap showing the relative abundance and distribution of genus-based OTU illumina reads. The color code indicates relative abundance, ranging from blue (low abundance) to black to brown (high abundance). 700 in YSGW11 were closely affiliated with clones recovered from macroalgal surface (GU451453 and GU451530) (Lachnit et al., 2011).

Correlation Analysis between Bacterial Taxonomic Groups and Environmental Parameters
In order to investigate whether some bacterial taxonomic groups are characteristics of specific environmental conditions, we calculated the Pearson's correlation coefficient between bacterial taxonomic groups and contextual environmental parameters (Gobet et al., 2012). The selected phyla/classes were divided into two groups: rare taxonomic group is defined as the phylum/class represented by more than 1% reads only from a single station, whereas abundant taxonomic groups is defined as the phylum/class represented by more than 1% reads from at least two stations (Table S1). Within abundant taxonomic groups, we found significant Pearson's correlation coefficients between isotope 224 Ra tracers and Phyla including Cyanobacteria, Betaproteobacteria, and Gammaproteobacteria. Salinity is an important environmental parameter which significantly influenced the Betaproteobacteria (−0.849). Pearson's correlation coefficient between Betaproteobacteria and nitrate (0.76)/ammonium (−0.774) suggested increased in sequences as nitrate concentration got higher and ammonium concentration got lower, indicating ammonium oxidation by Betaproteobacteria. Cyanobacteria was positively correlated with the ammonium concentration (0.863) suggesting nitrogen fixation by Cyanobacteria. Sequences within class Cytophagia had highly negative correlation with temperature but positive with DOC concentration, suggested organic matter degradation of this low temperature favorable bacterial group. Interestingly, six out of nine rare taxonomic groups were detected in porewater YSPW2, and these six groups had highly positive correlation with nitrite concentration (Figure 7).

Distinct Bacterial Communities within Porewater Samples
In this study, different locations may result in distinct bacterial communities among the porewater samples, for example, YSPW 2 was close to the regions where the green tides frequently break out, and there were specific bacterial taxonomic groups only detected in this site, but lacked some phyla which were abundant in other two porewater samples (Figure 5). One of the negative effects of the green tides on the marine ecosystem was the production of toxic hydrogen sulfide into surrounding environments (Lyons et al., 2014). The second most sequences in turbid YSPW 2 porewater sample were affiliated with JTB255 cluster, which was considered as the putative sulfide-oxidizing bacteria (Bowman and McCuaig, 2003). It is proposed that the toxic hydrogen sulfide could be detoxified in the porewater.
The Cyanobacteria was the most abundant group with OTU 874 comprising 30.4% in YSPW 7 and 7.6% in YSPW11, but only 0.97% in YSGW 2. Sequences of OTU 874 were remotely related to the heterocystous diazotrophic cyanobacterium Calothrix genus (Zehr, 2011), and 100% identical to the environmental clone PROA52S 10 from Karenia brevis dominated surface water (Jones et al., 2010), showing this uncultured cyanobacterial group might pose the putative ecological functions to different phytoplankton blooms. In addition, the most abundant OTU 67 in YSGW 11, had 97% similarity with the non-diazotrophic picocyanobacterium Prochlorococcus marinus strain SS12, suggesting the cohabitation of two groups of cyanobacteria in this sample. Unicellular phototroph Prochlorococcus marinus strain SS12 is well known for its characteristics of the tiny cell size, autotrophic metabolism and adaption to low light niche (Dufresne et al., 2003), this lowlight-adapted ecotype may make it possible to survive in the porewater at dimmer light.
It has been shown that microalgal growth was linked to nutrients released by the nitrogen fixing cyanobacteria. For instance, following aerolian iron input, nitrogen fixer Trichodesmium produced dissolved organic nitrogen (DON) for stimulating the growth of Karenia brevis, subsequently FIGURE 4 | Neighbor-joining tree showing phylogenetic relationships among the representative OTUs obtained in this study and reference 16S rRNA sequences retrieved from the NCBI GenBank for Actinobacteria and Betaproteobacteria. These OTUs were represented by more than 1% reads from a single station and/or from multiple stations within well water samples. The OTUs obtained in this study are shown in bold type. The numbers in parentheses indicate the percentage composition of reads in each station in the following order: (YSGW1, YSGW3, YSGW4, YSGW11) (YSPW2, YWPW7, YSPW11). The scale bar represents the estimated number of nucleotide changes per sequence position. The numbers at the nodes show the bootstrap values obtained after 1000 resamplings, only values of > 50 are shown. Accession numbers for the reference sequences are shown in parentheses.
caused the dinoflagellate bloom in the Gulf of Mexico waters (Lenes et al., 2001). Jones et al. (2010) found the higher abundance of cyanobacteria in the zero/low K. brevis libraries than those of medium/high K. brevis libraries. The Cyanobacteria was also detected in the surrounding water during an early Ulva bloom . Notably, cyanobacteria usually grow better at the early stage of macroalgae blooms, this could be explained that allelopathic chemicals produced by periphyton biofilms suppress cyanobacterial growth in macroalgae dominated waters (Wu et al., 2011). This may explained lower abundant cyanobacteria detected in YSPW 2. We speculated that the newly fixed N released by the Cyanobacteria in surface water and porewater might be significant in initiating the algal blooms. Considering high abundance of cyanobacteria detected in porewater samples, SGD became the potential nutrient source for adjacent coastal areas.
Bacterial community structures associated with the surface of Ulva prolifera and in surrounding seawater in Jiaozhou Bay during the Ulva bloom phases in 2008 were reported . They found Orders Alteromonadales, Flavobacteriales, and Rhodobacterales were the common bacterial groups on the surface of U. prolifera, and Alteromonaldales was also dominant in the seawater. The genus Glaciecola within Alteromonadales was dominated not only in the algal surface but also in the seawater. In porewater sample YSGW7, sequences affiliated with three orders accounted for 42.7% of total bacteria. Ra isotope analysis showed recirculated groundwater discharge in Group P samples, it was likely that these bacteria received from the seawater. An indoor mesocosm experiment with natural plankton communities from the western Baltic Sea revealed that Glaciecola sp. was the single most abundant taxon during the phytoplankton peak in the cold water (Scheibner et al., 2014). High abundance of Glaciecola cells has also been observed in FIGURE 5 | Neighbor-joining tree showing phylogenetic relationships among the representative OTUs obtained in this study and reference 16S rRNA sequences retrieved from the NCBI GenBank for Gammaproteobacteria, Cyanobacteria, and Planctomycetes. These OTUs were represented by more than 1% reads from a single station and/or from multiple stations within porewater samples. The OTUs obtained in this study are shown in bold type. The numbers in parentheses indicate the percentage composition of reads in each station in the following order: (YSGW1, YSGW3, YSGW4, YSGW11) (YSPW2, YWPW7, YSPW11). The scale bar represents the estimated number of nucleotide changes per sequence position. The numbers at the nodes show the bootstrap values obtained after 1000 resamplings, only values of > 50 are shown. Accession numbers for the reference sequences are shown in parentheses.
temperate North Sea during a spring bloom (Teeling et al., 2012). The considerable presence of Genus Glaciecola during the blooming phase of the diatoms in different studied regions suggested its significant roles in consuming dissolved organic carbon released by phytoplankton. YSGW7 contained 13.5% Glaciecola related sequences, indicated that SGD was significant source of the carbon flux to the coastal water. Meanwhile, Gammaproteobacteria was observed as dominated bacterial group during the phytoplankton peak (Teeling et al., 2012;Scheibner et al., 2014), the coexistence of Gammaproteobacteria and Cyanobacteria in our studied porewater sample was probably intensified contributions of the nutrient and DOC by SGD to the occurring of macroalgae bloom in the coastal areas.
Several studies of bacterial communities on the surface of marine macroalgae have shown the distinctive compositions of planctomycetes could be observed in association with different marine macroalgae (Lage and Bondoso, 2011;Bondoso et al., 2014). The second most abundant OTU 1051 in YSPW 11 was closely related to Planctomycetes Rhodopirellula Baltic strain SH 1 T , this strain was also frequently found in epiphytic communities of many macroalgae in a widespread geographic locations (Lage and Bondoso, 2011). The genome of Rhodopirellula baltic strain SH 1 T contains a large number of sulfatases, which are known to be involved in carbon recycling of complex sulfated heteropolysaccharides (Glöckner et al., 2003). The ability of Planctomycetes to degrade sulfatated heteropolysaccharides provides us the efficient pathway to obtain polysaccharides with lower molecular weights in porewaters, which were proved to possess important pharmacological activities such as antioxidant activities (Qi et al., 2005;Costa et al., 2010;Li et al., 2013). FIGURE 6 | Neighbor-joining tree showing phylogenetic relationships among the representative OTUs obtained in this study and reference 16S rRNA sequences retrieved from the NCBI GenBank for Bacteriotedes. These OTUs were represented by more than 1% reads from a single station and/or from multiple stations within well water or porewater samples. The OTUs obtained in this study are shown in bold type. The numbers in parentheses indicate the percentage composition of reads in each station in the following order: (YSGW1, YSGW3, YSGW4, YSGW11) (YSPW2, YWPW7, YSPW11). The scale bar represents the estimated number of nucleotide changes per sequence position. The numbers at the nodes show the bootstrap values obtained after 1000 resamplings, only values of >50 are shown. Accession numbers for the reference sequences are shown in parentheses.

The Putative Ecological Role of Genus Limnohabitans in the Fresh Well Water
Betaproteobacteria and Actinobacteria constituted the two major groups in four fresh well water samples, in agreement with the finding that a few phylogenetic clusters are dominant in a variety of freshwater ecosystems (Allgaier and Grossart, 2006;Newton et al., 2011). Our correlation analysis confirmed Betaproteobacteria favored lower salinity condition (Figure 7).
Members of Genus Limnohabitans within Betaproteobacteria were frequently predominant in freshwater bacterioplankton communities (Šimek et al., 2010b). We reported that Bacteria of the genus Limnohabitans were the major group in the fresh well samples (Figure 3). Limnohabitans spp. are the wellstudied opportunistic bacterial groups (Šimek et al., 2005), with fast living generation time, specialization on low molecular weight dissolved organic matter, and highly vulnerability to size selective predation (Salcher, 2013). Autochthonous algalderived organic matters and products of the photolysis of dissolved materials are two main groups of substrates for the growth of members of Genus Limnohabitans (Šimek et al., 2011). Genus Limnohabitants was identified as an important player in in carbon flow to higher trophic levels since they can be consumed by food for heterotrophic nanoflagellates (Šimek et al., 2013). More recent study showed Limnohabitans Species Strains Rim28 and Rim47 had a great metabolic versatility, including photosynthesis, autotrophic carbon fixation, ammonium oxidation and sulfur oxidization (Zeng et al., 2012b). The higher nitrate concentrations in YSGW3 and YSGW 11 was in good agreement with the most abundant sequences affiliated with Limnohabitans sp. Rim 28, which may FIGURE 7 | Environmental factors associated with variations of the bacterial community structure at the phylum level. In order to obtain a higher resolution, the Proteobacteria phylum level was separated into Alpha, Beta, and Gamma-Proteobacteria classes. Phylum Bacteriodetes was separated into classes Flavabacterria and Cytophagia. The total number of sequences in each phylum is indicated in parentheses. Pearson's correlation coefficients between -1 and 1 are shown in the rectangle, which indicates correlations between phylum/class sequence abundance and selected environmental parameters. For example, a firebrick colored rectangle (0.76) between nitrate and Betaproteobacteria indicates a higher number of sequences with increasing nitrate concentration, a blue rectangle (−0.774) between Betaproteobacteria and ammonium indicates a higher number of sequences with decreasing ammonium concentration. The color code indicates Pearson's correlation coefficients, ranging from blue (−1) to white (0) to firebrick (1). The density showed the distribution of Pearson's correlation coefficients between −1 and 1. Si, silicate; P, phosphate; NO 2 , nitrite; NO 3 , nitrate; NH 4 , ammonium; 224 Ra, Radium isotope tracer; DOC, dissolved organic carbon.
serve as ammonium oxidizer in these two samples. Fresh well waters provide another potential sources for both nitrogen and carbon flux to the coastal waters.

Bacterial Candidate for Bioremediation in Well Waters
Fresh well sample YSGW 4 is located in a public park in Shandong province, garbage such as mosquito and discarded beverage bottles in the well were observed during the sampling time. Compared to the other three fresh water samples, YSGW 4 exhibited relatively lower nitrate and higher DOC concentrations. The co-occurrence of sequences closely related to Actinobacteria Rhodoluna lacicola strain MWH-EgelM2-3.2 (FJ545223), Flavobacterium granuli strain kw05 (NR_041052), and uncultured Flectobacillus clone ZS-2-379 (FN668111) (Figure 4). The presence of actinorhodopsins in Rhodoluna lacicola strain MWH-EgelM2-3.2 suggested this strain to be photoheterotrophic bacteria, surviving under low nutrient and energy conditions (Sharma et al., 2008). Flavobacterium granuli was isolated from granules used in a wastewater treatment plant (Aslam et al., 2005); Members of genus Flectobacillus usually act as the direct competitors in situ (Šimek et al., 2010a). The characteristics of these three groups indicated the hostile environment of YSGW 4 sampling site. The most abundant sequences in YSGW 4 were affiliated with C. testoteroni strain CNB-2. The complete genome of strain CNB-2 revealed its genetic versatility adaption to diverse habitats, based on its ability to metabolize a variety of compounds, such as carboxylic acids and aromatic compounds, utilize nitrate and ammonia as nitrogen source, resistant to drug and heavy metals (Ma et al., 2009). C. testosteroni strain CNB-1 has been used successfully rhizoremediation of 4-chloronitrobenzene (CNB) polluted soil (Liu et al., 2007). Other Comamonas spp. were also used for environmental applications (Tobajas et al., 2012;Wu et al., 2015). The high abundance of C. testosteroni related sequences in our samples indicated that this species could be selected as potential microbes for bioremediation in the highly polluted fresh well waters before flowing into coastal water.

Linking Bacterial Community with Environmental Variables
Based on our analysis of relationship between bacterial taxonomic groups and contexture environmental parameters, we found Betaproteobacteria in well water had highly negative correlation with ammonium, but had positive correlation with nitrate; in contrast, Cyanobacteria in porewater had highly positively correlation with ammonium. These results supported the possibilities that fresh well water and brackish porewater may provide nitrate and ammonium respectively. Only YSPW 2 contained more than 1% sequences within six rare bacterial taxonomic groups (Figure 7), which had highly positively correlation with nitrite, showing these groups in YSPW2 had the abilities to cycle nitrogen.
Four well water samples were collected from three different eco-environments impacted by heavy human activities (Table 1). We detected less diverse bacterial groups in the well water samples than those in the porewater samples ( Table 2), suggesting the stressful conditions in the well water due to anthropogenic activities. Comamonas spp. and Limnohabitans spp. were dominated in fresh well samples, but were almost absent from the brackish porewater samples, which indicates these groups favored lower salinity environment. Distinctive bacterial patterns were detected among the three porewater samples collected from the different tourist beaches, for example, six rare bacterial taxonomic groups were mainly found in YSGW 2, which may have been due to niche selections.

CONCLUSION
Our study documented bacterial community structures in the both fresh well water and brackish porewater along the coasts of the Yellow Sea. We provided perspective into the correlations between representative bacterial taxonomic groups and environmental parameters, fresh well water and brackish porewater provided the different nitrogen sources to the coastal water. Potential key bacterial groups such as Comamonas testosteroni may be excellent candidates for bioremediation of the natural pollutants in different niches. Further studies will explore the shifts of indigenous microbial community in in-situ bioremediation experiments in the submarine groundwater.

AUTHOR CONTRIBUTIONS
JD designed the experiments; JL collected the samples; QY and JL performed the experiments and analyzed the data; QY, JL, JD, and JZ wrote the manuscript.