Distinct Seasonal Patterns of Bacterioplankton Abundance and Dominance of Phyla α-Proteobacteria and Cyanobacteria in Qinhuangdao Coastal Waters Off the Bohai Sea.

Qinhuangdao coastal waters in northern China are heavily impacted by anthropogenic and natural activities, and we anticipate a direct influence of the impact on the bacterioplankton abundance and diversity inhabiting the adjacent coastal areas. To ascertain the anthropogenic influences, we first evaluated the seasonal abundance patterns and diversity of bacterioplankton in the coastal areas with varied levels of natural and anthropogenic activities and then analyzed the environmental factors which influenced the abundance patterns. Results indicated distinct patterns in bacterioplankton abundance across the warm and cold seasons in all stations. Total bacterial abundance in the stations ranged from 8.67 × 104 to 2.08 × 106 cells/mL and had significant (p < 0.01) positive correlation with total phosphorus (TP), which indicated TP as the key monitoring parameter for anthropogenic impact on nutrients cycling. Proteobacteria and Cyanobacteria were the most abundant phyla in the Qinhuangdao coastal waters. Redundancy analysis revealed significant (p < 0.01) influence of temperature, dissolved oxygen and chlorophyll a on the spatiotemporal abundance pattern of α-Proteobacteria and Cyanobacteria groups. Among the 19 identified bacterioplankton subgroups, α-Proteobacteria (phylum Proteobacteria) was the dominant one followed by Family II (phylum Cyanobacteria), representing 19.1-55.2% and 2.3-54.2% of total sequences, respectively. An inverse relationship (r = -0.82) was observed between the two dominant subgroups, α-Proteobacteria and Family II. A wide range of inverse Simpson index (10.2 to 105) revealed spatial heterogeneity of bacterioplankton diversity likely resulting from the varied anthropogenic and natural influences. Overall, our results suggested that seasonal variations impose substantial influence on shaping bacterioplankton abundance patterns. In addition, the predominance of only a few cosmopolitan species in the Qinhuangdao coastal wasters was probably an indication of their competitive advantage over other bacterioplankton groups in the degradation of anthropogenic inputs. The results provided an evidence of their ecological significance in coastal waters impacted by seasonal inputs of the natural and anthropogenic matter. In conclusion, the findings anticipate future development of effective indicators of coastal health monitoring and subsequent management strategies to control the anthropogenic inputs in the Qinhuangdao coastal waters.

Introduction

Marine and coastal environments are of major concern nowadays on account of the rapid developments and growth in anthropogenic activities. Most commonly the anthropogenic activities include terrestrial pollution, aquaculture, urban development, tourism, maritime transport, agricultural and industrial activities, oil refineries, and mining, etc. (Islam and Tanaka, 2004). The nutrient inputs from these activities are the major drivers of eutrophication in most of the coastal areas, ultimately affecting the marine food web components (Vaquer-Sunyer et al., 2015). In addition, rivers and allochthonous biological processes also enrich coastal water with dissolved compounds (Hiriart-Baer et al., 2008; Barrera-Alba et al., 2009; Vargas et al., 2011; Gustafsson et al., 2014; Lafon et al., 2014). Thus, organic matter mineralization and nutrient cycling are crucial for maintenance of coastal ecosystem health. In this context, the assemblage and abundance of bacterioplankton involved in the marine microbial loop (Farooq and Francesca, 2007) are canonical indicators of ecosystem status and the extent of anthropogenic impact on the coastal waters (Zhang et al., 2014). The spatial and temporal dynamics of bacterioplankton community is driven by multiple environmental factors such as latitudinal gradient (Fuhrman et al., 2008; Milici et al., 2016), temperature (Kim et al., 2016), predation (Oguz et al., 2013; Karus et al., 2014), nutrients (Kim et al., 2016), succession (Chen H. et al., 2016), seasons (Gilbert et al., 2012), and timescales (Fuhrman et al., 2015). Apart from natural environmental factors, anthropogenic activities also shape the bacterioplankton assemblage and abundance pattern (Zhang et al., 2009; Thiyagarajan et al., 2010; Zhou et al., 2011; Sauret et al., 2012; Jeffries et al., 2016; Meziti et al., 2016). For instance, the abundances of the bacterioplankton groups that increased and decreased in the impacted sites were significantly correlated with nutrients enrichment (Fodelianakis et al., 2014). Also, there is a direct evidence of specific metal (Cd)-induced patterns in bacterioplankton communities in coastal systems (Wang et al., 2015). In oligotrophic coastal water microcosm with nitrate perturbation, the bacterioplankton community composition was greatly influenced by nitrate loading mode, indicative of nitrate loading impact on marine environments (Dong et al., 2017). Thus, nutrients and toxic pollutants enrichment seem to be the important drivers of the bacterioplankton composition in sites receiving direct human impact (Aguilo-Ferretjans et al., 2008).

A closer inspection of the bacterioplankton community dynamics reveals many interesting relationships and patterns in the taxonomic diversity and abundance. A high concentration of Synechococcus and low concentration of Prochlorococcus has been observed in sites impacted by human activities with inputs from the land (Bouvy et al., 2012). In a microcosm study, it was found that anthropogenic impacts do not necessarily influence the abundance pattern of rare species (<0.1% of relative abundance), and in fact, nutrients enrichment increased the relative abundance of only abundant species (>1% of relative abundance) Polaribacter, Tenacibaculum, and Rhodobacteraceae (Baltar et al., 2015). By the use of 16S rRNA gene cloning and sequencing, the bacterioplankton diversity of oligotrophic marine environments is relatively well studied than coastal environments (Zhou et al., 2011; Ortega-Retuerta et al., 2014; Hartmann et al., 2016; Milici et al., 2016; Dong et al., 2017). Bacterioplankton diversity and abundance patterns on a regional-scale and also for local sites under anthropogenic impact remains a hot topic and several reports are available (Simonato et al., 2010; Niu et al., 2011; Peierls and Paerl, 2011; Zhou et al., 2011; Fodelianakis et al., 2014; Laas et al., 2014; Almutairi, 2015; Baltar et al., 2015; Meziti et al., 2015; Xiong et al., 2015; Jeffries et al., 2016). Overall, the available reports on bacterioplankton community composition at regional-scale or in local sites indicate diverse patterns and complex relationships between species composition and environmental factors. In addition, phylogenetic analysis and abundance patterns of bacterioplankton assemblages provide significant insights into the bacterioplankton dynamics in the impacted coastal areas which further leads us to the development of coastal environmental monitoring and management strategies (Demarcq et al., 2012). This calls for more research to identify the bacterioplankton species that show strong linkages with the implicit environmental factors for developing effective descriptors of the anthropogenic impacts.

Qinhuangdao is a famous port city, an important energy export port, and coastal tourist spot on the northwest coast of the Bohai Sea in North China with an area of approximately 7812 km2 and a population of 2.9 × 106 (Gu et al., 2016). The Bohai Sea coastline experiences the strong influence of riverine systems, aquaculture, municipal sewage, industrial wastewaters, and agricultural runoffs (Zhu et al., 2014). Land-use surrounding the Qinhuangdao coastal area shifted drastically from agriculture to industrialization in the last 30 years (Xu et al., 2010). The area has undergone development of several petrochemical, steel, and other projects, besides land reclamation and construction of embankments (Zhu et al., 2014). The resulting land-use change has overall influenced the coastal seawater nutrients cycling and associated bacterioplankton, phytoplankton, and other biotic factors (Xu et al., 2010; Liu et al., 2011). Results from other and our previous studies have shown a distinct pattern of bacterioplankton abundance and diversity as a consequence of environmental influence such as estuaries, mariculture, and rainfall (Li J. et al., 2011; Li et al., 2013; Shang et al., 2016; Zhou et al., 2016). These reports further encourage us to investigate closely the relationship between the bacterioplankton community and nutrient levels across different seasons in the near-shore and off-shore stations near Qinhuangdao coast and address the following questions. What are the patterns in abundance of major taxonomic groups and their spatial diversity? What drives the spatiotemporal variations of bacterioplankton in the Qinhuangdao coastal waters?

The objectives of the present study were to: (1) analyze the abundance of total bacteria and major bacterioplankton groups across different seasons, (2) identify the key environmental factors that influence the bacterioplankton abundance, and (3) assess the bacterioplankton distribution and diversity along the near-shore and off-shore coastal areas of Qinhuangdao.

Materials and Methods

Study Area and Water Collection

Qinhuangdao, one of the highly urbanized regions, is a port city on the coast of China located in the northeastern Hebei province, and it is situated about 280 km east of Beijing on the northwest Coast of Bohai Sea, the innermost gulf of Yellow Sea. The Qinhuangdao coastal area has a humid continental climate with four distinct seasons. It has a monsoon-influenced humid continental climate with the highest precipitation (152–189 mm) during July–August period. From October month the temperature starts to fall (<5°C) until April, and a very low (<10 mm) precipitation occurs during November to March period (Gu et al., 2016).

Six stations (W1–W6) were selected that had the distinct influence of the local environment (Figure ). The W1 station is located outside the port of Qinhuangdao mainly polluted by the shipping industry and municipal sewage, the W3 station is located in the Xinkai river estuary near the Beidaihe Forest Wetland and Qinhuangdao Wildlife Park, and the W5 station is located between the Yang River Estuary and Dai River Estuary. The Yang River and Dai River are severely polluted by microorganisms and the total number of bacteria exceeds the standard levels. W3 and W5 stations are also located in the tourist area and near the bathing beach. W2, W4, and W6 are off-shore stations with aquaculture facility near W6 station. Water samples were collected in the month of March, July, October, and December (2014) at the stations W1, W2, W3, W4, W5, W6 (Figure ). For each month one water sample (1 L) at 1 m depth from the surface was collected using a sample hydrophore. A plexiglass hydrophore (JC-800, Juchuang, China) of capacity 5 L, diameter 15 cm, and height 37.5 cm was used as the sample hydrophore. All samples were stored in sterile glass bottles at 4°C and then transported to the laboratory. A 500 mL volume of water sample was filtered using 0.22 μm membrane filter (AmeriTech Inc., United States). The membrane filters coated with microorganisms were stored at -80°C.

Map of sampling stations.

Water Quality and Nutrient Analysis

The pH, temperature, dissolved oxygen (DO), and salinity of the water samples at six stations were measured using a portable YSI Pro Plus Multiparameter instrument (YSI Inc., United States). The contents of total nitrogen (TN), total phosphorus (TP), ammonia nitrogen (NH4-N), nitrate nitrogen (NO3-N), nitrite nitrogen (NO2-N) and Chlorophyll a (Chl a) were determined in water samples in accordance with the “Specification for oceanographic survey” (GB/T 12763.4-2007). Briefly, NH4-N was measured by sodium bromate oxidation method which involves oxidization of ammonium salt to nitrite in alkaline medium, followed by determination of the TN content using diazo spectrophotometry. After deducting the concentration of NO2-N, the concentration of NH4-N was estimated. NO3-N was measured by zinc cadmium reduction method where nitrate is quantitatively reduced to nitrite in water by zinc cadmium reduction method following determination of total nitrite by the diazo and azo method, the original nitrite was then corrected and the nitrate content was calculated. NO2-N was measured by a diazo azo method which involves the reaction of nitrite with sulfanilamide in acidic medium following product reaction with hydrochloric acid and synthesis of a red azo dye whose absorbance was measured at 543 nm. Chl a measurement was done by an extraction-fluorimetric method which involves excitation by blue light of the acetone extract of Chl a to produce red fluorescence (Holm-Hansen and Bo, 1978). The fluorescence value of the extract was determined at 685 nm before and after acidification.

The abundance of total bacteria at the 6 stations near Qinhuangdao coastal area over warm and cold seasons were monitored using fluorescence microscopy to assess the influence of anthropogenic and natural activities in the near-shore and off-shore waters. The total bacterial count in water samples (formaldehyde-fixed) at each station was estimated by fluorescence microscopy (Eclipse Ni-U, Nikon Instruments Inc., United States). Sample (2 mL water) filtration was done at a pressure of 150 mm Hg through a polycarbonate filter (0.22 μm aperture, 25 mm diameter) (AmeriTech Inc., United States). The membrane filter was stained with 2X SYBR-Gold solution (300 μl) and was visualized at 480–485 nm wavelength (blue light excitation). Data reported are the means of counts in 20 randomly selected fields per sample. For each station, triplicate samples were processed and for each replicate sample, bacterial numbers were counted in 20 randomly selected fields. The data reported are the means of counts in triplicate samples.

DNA Isolation and Sequence Analysis of Bacterioplankton Assemblage

Five hundred milliliter of surface water samples collected in July 2014 from 6 stations were filtered onto 0.2 μm polycarbonate membrane filters (AmeriTech Inc., United States). The resulting filters were stored in ultra-low temperature refrigerator (-80°C) until DNA isolation. The total DNA was extracted using E.Z.N.A.TM Water DNA Kit (Omega Bio-Tek Inc., United States) following the manufacturer’s instructions. The isolated DNA was used as template for bacterial 16S rRNA gene amplification using universal primers 27F (5′-AGAGTTTGATCCTGGCTCAG-3′) and 1492R (5′-GGTTACCTTGTTACGACTT-3′) on Bio-Rad S1000 thermal cycler (Bio-Rad Laboratories Inc., United States) using the following protocol: 95°C 5 min; 94°C 45 s, 55°C 45 s, 72°C 1 min; 35 cycles; 72°C 10 min. PCR amplification for each sample was performed in triplicates and the triplicate PCR products were mixed to obtain the final sample for clone library construction. All PCR products were gel-purified, cloned into a pMD18-T Simple Vector (Sino Biological Inc., China), and then transformed into Escherichia coli DH5α competent cells by traditional heat shock treatment. 16S rRNA clone libraries of a total of six water samples (one sample per station) for the month of July were constructed and 90 positive clones were randomly picked from each sample plate for sequence analysis (Beijing Genomics Institution, China). The coverage of the clone libraries ranged from 31.1 to 56.7% and the impact of low coverage on diversity analysis was minimized by calculating Simpson diversity estimate (Haegeman et al., 2013) in Mothur software package (Schloss et al., 2009). The sequences have been submitted to NCBI Genebank database with accession number MF498066 – MF498485.

All sequence analyses were performed within versions 1.35.1 of the Mothur software package (Schloss et al., 2009). The FASTA format sequences were aligned with SILVA alignments and the distance was calculated at cutoff 0.03 (97% sequence homology) using align.seqs and dist.seqs commands, respectively. Chimeras were identified using chimera.uchime command with default parameters. Sequences were classified using classify.seqs command with reference to Ribosomal Database Project (RDP) as the taxonomy file (silva.bacteria.rdp.tax) and the remove.lineage command was used to remove the sequences that belong to mitochondria, chloroplast, eukaryote, archaea, and unknown lineages. Cluster command was used to assign sequences to OTUs with opticlust clustering method. α-diversity was calculated by counting the number of OTUs and using the reciprocal of Simpson Index (invsimpson). The invsimpson calculator is preferred to other measures of α-diversity as it indicates the richness in a community with a uniform evenness that would have the same level of diversity. Moreover, the inverse of the Simpson index has some biological interpretation and do not tend to be affected by sampling effort because it is independent of abundance distributions (Haegeman et al., 2013).

Quantitative PCR Analysis of Bacterioplankton Groups

To assess the phylum/class level abundance distribution of total bacteria in warm and cold seasons, qPCR experiments with phylum and class specific primers (Supplementary Table S1) were conducted for all 6 stations. qPCR method was used to detect and quantitate the abundance of the major phylogenetic groups of bacterioplankton in the water samples from Qinhuangdao coastal area. The advantage of qPCR method in our study was its ability to quantitate specific bacterioplankton groups by using phylum/class specific primers. In addition, qPCR quantitation is fast, reliable, and accurate. Based on the results of clone library analysis and identification, the following dominant classes were selected to carry out quantitative PCR experiments: α-Proteobacteria, β-Proteobacteria, Cyanobacteria, Actinobacteria, Firmicutes, and Bacteroidetes. PCR was carried out in a 20 μL reaction mixture containing 1 μL DNA template (2.8-10.5 ng/μL), 10 μL SYBR Select Master Mix (Thermo Fisher Scientific Inc.), 1 μL forward primer (10 nmol/mL), 1 μL reverse primer (10 nmol/mL), and 7 μL nuclease-free molecular-grade water. The qPCR cycling condition was: 95°C 55 s; 95°C 20 s, 60°C 35 s, 72°C 30 s; 40 cycles; 95°C 1 min, 60°C 1 min, heating-up from 75 to 95°C, the increment of 0.5°C for 8 s. The annealing temperature of each primer is shown in Supplementary Table S1. The standard curves were constructed with known amounts of plasmid DNA containing the sequence insert. For the positive plasmid preparation, the PCR products were ligated to pMD18-T Vector by pMD18-T Vector Cloning Kit (Takara Biotechnology (Dalian) Co., Ltd.) and then transformed into E. coli competent cells (DH5α) followed by a screening of positive clones. Plasmid DNA was extracted using E.Z.N.A® Plasmid Midi Kit (Omega Bio-Tek Inc., United States) and the concentration of the extracted DNA was measured using NanoDrop ND-1000 spectrophotometer (Nanodrop Technology). The gene copies were calculated using the mean mass of pMD18-T Vector ligated with the target DNA sequence. DNA of each sample and standard plasmid DNA dilution (101 and 105) were used as templates for qPCR. Each sample and standard were analyzed in triplicates. The standard curve was constructed to estimate the copy number of each sample according to the dilution ratio and copy number of the standard plasmid template.

Statistical Analysis

Principal component analysis (PCA) was done to reveal the relationships between the environmental parameters and the distribution of the sampling stations in an ordination plot. Pearson’s correlation analysis was performed to assess the correlation between the abundance of total bacteria, major bacterioplankton groups, and environmental parameters. Redundancy analysis (RDA) was carried out to identify the explanatory environmental parameters that significantly (p < 0.01) constrain the bacterioplankton abundance data by forward selection in Canoco version 5.02. RDA was performed using abundance (qPCR data) and environmental data. The PCA of environmental parameters was done using normalized (z-score) data in Canoco version 5.02. For RDA, the qPCR data were log(x+1) transformed and the environmental data was normalized by z-score.

Results

High Seasonal Variations of Environmental Parameters

The environmental parameters across 6 stations near Qinhuangdao coastal area appeared highly heterogeneous spatiotemporally. In the warm season (July–October), the temperature of coastal waters in all stations was within 21–25°C, whereas in the cold season (December–March) the temperature was significantly lower ranging from -0.3 to 6.8°C. pH was notably higher (11.6 ± 0.09) in October at all stations. Salinity levels were similar (∼29 ppt) in all stations and also during warm and cold seasons. The DO levels were relatively higher in cold period than warm as a result of temperature effect on oxygen solubility in water. The average values of environmental parameters (nutrients) with their variation over different periods at the 6 stations near Qinhuangdao coastal area are shown in Table . The NH4-N level was high in W5 station than that of the others. In October the NH4-N was relatively high at all stations (0.05–0.17 mg/L). No significant variation in NO2-N and NO3-N content was observed at the 6 stations. TN was detected in all stations (0.35–0.6 mg/L) with relatively large fluctuations (0.6–1.8mg/L) in the warm season (July) and was lowest in the cold season (December). In July and March, the TP was found to be comparatively higher than in October and December. Furthermore, the average TP content was higher in W1, W2, W3, W5 (0.12–0.13 mg/L) than that in W4 and W6 stations (0.07–0.08 mg/L). W1 and W5 stations showed highest average Chl a level (11.6–16.0 mg/m3) followed by W3 and W6 (5.6–5.9 mg/m3) and lowest levels at W2 and W4 (2.6-3.2 mg/m3). The variation in Chl a level was substantially less in W2 and W4 stations suggesting a low abundance of phytoplankton at these stations than at the polluted stations (W1, W3, W5, and W6). In the Bohai Sea, silicates and Si/dissolved inorganic N ratio have been associated with phytoplankton abundance (Chen Y.H. et al., 2016). In addition, climate and mariculture activity have a significant correlation with the Chl a trends (Fu et al., 2016). River discharge and suspended sediment also influence Chl a in Bohai Sea coast. These factors might have played a role in the spatial variation of Chl a observed in the present study. Also, Chl a mostly was higher in July in all stations probably due to more growth of phytoplankton resulting from optimal light and temperature conditions.

By analyzing all the measured environmental variables in combination, the resultant ordination plot showed distinct partitioning of coastal water samples across warm and cold seasons rather than across near-shore and off-shore stations (Figure ). The ordination plot could explain 59.6% of total variation in the environmental data and revealed a linear positive correlation between TN, temperature, and Chl a, and a negative correlation of DO with these variables. These parameters were most crucial in the partitioning of samples from the warm and cold seasons. Similarly, NH4-N and NO3-N had a positive correlation with each other, and salinity, NO2-N, and TP did not show any notable correlations to any of the other variables. Overall, the PCA biplot clearly showed that seasonal variation of the measured environmental parameters was more than spatial variation. Thus, seasonal loads of anthropogenic inputs seem to have a significant impact on the nutrients conditions in the coastal areas of Qinhuangdao.

Biplot of the principal component analysis (PCA) for environmental parameters in the sampling stations near Qinhuangdao coastal area. The two principal components (PC1 and PC2) explained 59.64% of total variation in environmental data, and shows clear partitioning of July and October samples from December and March along the PC1, and partitioning of July from October along PC2. PC1 and PC2 strongly correlated to variables NH4-N and NO3-N, and total nitrogen (TN), temperature (T) and Chlorophyll-a (Chl a), respectively. Total phosphate (TP), NO2-N, and salinity were not correlated to PC1 and PC2. Blue arrows represent environmental parameters, and circles in color represents sampling points.

Bacterioplankton Abundance Patterns and Influencing Factors

The results of fluorescence microscopy based detection of total bacterial abundance in the near-shore and off-shore stations are shown in Figure . In July, all the stations mostly showed highest abundance (5.32 × 105 – 1.71 × 106 cells/mL) of bacteria compared to other months. Overall, the total bacterial abundance in each of the 6 stations over the 4 months was in between 8.67 × 104 and 2.08 × 106 cells/mL. Within the other 3 months, March recorded higher abundance than October and December in all the stations. The stations W1, W3, W5, and W6 typically had relatively more abundance than the W2 and W4 perhaps due to nutrients overload in these stations as a consequence of higher anthropogenic and natural activities as evident from their higher level of Chl a. A significant (p < 0.01) positive correlation of total bacterial abundance with TP was seen (Figure ).

Fluorescence microscopy based quantification of total bacteria at the near-shore and off-shore stations near the Qinhuangdao coastal area of Bohai Sea. ‘W1’ to ‘W6’ indicate the different sampling stations.

Significant (p < 0.01) correlation between total bacterial abundance and total phosphorous. Values on Y-axis are shown in log-scale.

The qPCR based abundance estimates of 6 major phylum/class i.e., α-Proteobacteria, β-Proteobacteria, Cyanobacteria, Actinobacteria, Firmicutes, and Bacteroidetes are shown in Figure . α-Proteobacteria, Actinobacteria, Firmicutes, and Bacteroidetes showed a fairly similar pattern of abundance distribution in the 6 stations. Notably, in the warm season, these 4 groups exhibited relatively higher abundances than in cold season. However, β-Proteobacteria and Cyanobacteria showed distinct trends, wherein the former one showed a similar pattern in all stations over the 4 months, i.e., higher abundance in cold than warm season. Bacteroidetes abundance at W1 and W6 stations remained high in October, and in W5 its abundance increased only in July and December. Near-shore and off-shore differences in the abundance of the major phylotypes were also apparent. In the warm season, near-shore stations possessed more Cyanobacteria and Bacteroidetes than off-shore stations. Alternatively, in cold months, near-shore stations possessed a relatively higher abundance of α-Proteobacteria, Firmicutes, Actinobacteria, and β-Proteobacteria than off-shore stations.

16S rRNA gene abundance (qPCR quantification) of major phylogenetic groups from coastal bacterioplankton assemblages in 6 stations near Qinhuangdao coastal area. The symbols and lines in black, red, green, blue, cyan, and magenta represent stations W1, W2, W3, W4, W5, and W6. Values on Y-axis are plotted in log-scale.

Table shows the Pearson’s correlation coefficients for the relationship between bacterioplankton abundance and some environmental parameters. α-Proteobacteria showed significant (p < 0.05) positive correlation with TN, Chl a, DO, and temperature. Furthermore, from the ordination plot (Figure ) based on redundancy analysis, temperature and Chl a were the major factors that positively correlated with the first principal component and were the significant explanatory parameters (p < 0.01) that constrained the abundance of α-Proteobacteria and Cyanobacteria phylotypes. On the contrary, DO showed a negative influence (p < 0.01) and salinity was non-significant in explaining any variation of the abundance of bacterioplankton groups. Overall, samples collected in warm July or October had a much higher abundance of α-Proteobacteria, Firmicutes, and Actinobacteria than that of cold March or December. In contrast, β-Proteobacteria, Cyanobacteria, and Bacteroidetes did not show any distinct seasonal influence which likely suggests their diverse metabolic potential and ability to utilize organic matter even at low concentrations that allow them to attain invariable growth throughout the year.

Biplot of redundancy analysis (RDA) showing relationship between the environmental variables and 16S rRNA gene abundance (qPCR quantification) of major phylogenetic groups in six stations near Qinhuangdao coastal area. The two RDA axes (RDA1 and RDA2) explained 70.35% of total variation in abundance data, and the significant (p < 0.01) explanatory variables were temperature (temp), dissolved oxygen (DO), and Chlorophyll a (Chl_a). Blue arrows represent different phylotypes, gray arrows represent environmental parameters, stars represent near-shore samples and circles represent off-shore samples. Different colors of the star or circle symbols represent different sampling time, with March samples in green, July samples in red, October samples in yellow and December samples in purple.

Spatial Distribution and Diversity of Bacterioplankton Assemblages

The assemblage of bacterioplankton communities featuring diversity and distribution at the 6 stations in the warm season was assessed and the results are reported in Table and Figure . Based on the environmental data, July exhibited greater levels and fluctuations in the measured parameters perhaps due to more nutrient inputs that resulted from anthropogenic and natural activities in warm season. Thus, further characterization of bacterioplankton phylogeny, composition, and diversity was conducted on water samples collected in July. Out of the selected clones analyzed, 62, 62, 68, 44, 47, and 71 OTUs were identified in W1, W2, W3, W4, W5, and W6 stations, respectively (Table ). A broad phylogenetic distribution of the bacterioplankton community at the 6 stations in warm season is shown in Figure . The analysis of 16S rRNA gene sequences from these stations revealed 7 distinct phyla (Figure ). The majority of the sequences generated from the 6 bacterioplankton clone libraries were affiliated to phyla Proteobacteria (21.7–75.9%), Cyanobacteria (4.6–54.2%), Bacteroidetes (9–21.2%), Actinobacteria (1.1–5.6%), Firmicutes (0–4.4%), Verrucomicrobia (0–3.5%), Planctomycetes (0–3.3%). The unclassified group represented sequences in the range of 9.2–26.7%. At the phylum level Proteobacteria (49% of total sequences), Cyanobacteria (24% of total sequences), and Bacteroidetes (17% of total sequences) were the dominant groups in all the stations. There were 19 bacterioplankton classes in the total sequences acquired from the coastal waters of Qinhuangdao (Figure ). The predominant phylum Proteobacteria comprised of five classes namely α-Proteobacteria, γ-Proteobacteria, 𝜀-Proteobacteria, β-Proteobacteria, and δ-Proteobacteria accounting to ca. 82, 9, 3.6, 2, and 1.6% of Proteobacteria, respectively. Further taxonomic resolution of the assemblages revealed nearly 30 different families that were assigned to the identified OTUs from the bacterioplankton assemblages (Supplementary Figure S1). The OTUs within Proteobacteria phylum were assigned to 15 families: Campylobacteraceae, Comamonadaceae, Methylophilaceae, Pseudobacteriovoracaceae, Burkholderiaceae, Bacteriovoracaceae, Moraxellaceae, Halomonadaceae, Oceanospirillaceae, Enterobacteriaceae, Rhodospirillaceae, SAR11, Hyphomicrobiaceae, Rhodobiaceae, and Rhodobacteraceae. The second major phylum Cyanobacteria comprised OTUs which were all assigned to a single family — Family II. The OTUs representing phylum Bacteroidetes belonged to 4 families: Flavobacteriaceae, Cryomorphaceae, Saprospiraceae, and Chitinophagaceae. Rhodobacteraceae represented major proportion (>10% of total 16S r RNA sequences) in all the stations followed by Family II especially in W1 and W4 stations. Only three OTUs assigned to Rhodobacteraceae, Family II, and Cryomorphaceae were cosmopolitans, i.e., they were found in all the stations. Overall, α-Proteobacteria (Rhodobacteraceae) (19.1–55.2%) and Cyanobacteria (Family II) (2.3–54.2%) were the most dominant cosmopolitan groups. Interestingly, it was noted that these two dominant groups had an inverse relationship (r = -0.82).

Composition and distribution of bacterioplankton assemblages at the 6 stations near Qinhuangdao coastal area. Percentage distribution of the 16S rRNA gene sequences at the (A) phylum level and (B) class level indicates higher relative abundance of subgroups Proteobacteria, Cyanobacteria, Bacteroidetes, and Actinobacteria.

Based on inverse Simpson index (1/D), the ascending order of diversity at the 6 stations was W5 < W4 < W1 < W6 < W2 < W3 ranging from 10.2 to 105 (Table ). The results suggested a wide difference in the bacterioplankton diversity, which likely resulted from the presumed shifts in the environmental conditions of the coastal waters of Qinhuangdao. The inverse Simpson index (1/D) may serve as a sensitive indicator of the anthropogenic impact in the Bohai Sea and can accurately monitor the extent of human activities in Qinhuangdao coastal region.

Discussion

Spatiotemporal Variations and Factors Influencing Bacterioplankton Abundance

Most global-scale studies on the biogeography of bacterioplankton communities have revealed intriguing findings and have largely shown that diversity and abundance pattern are closely linked to latitudinal gradient (Fuhrman et al., 2008), seasons (Cuevas et al., 2004; Fuhrman et al., 2015), temperature and nutrient levels (Pommier et al., 2007; Vichi et al., 2007). In addition, previous studies have well established that bacterioplankton abundance is potentially linked to elevated anthropogenic nutrient loads in aquatic ecosystems. Overall, nutrient contents (TP, TN, Chl a, and DO) and gradient of the natural parameter (temperature) displayed significant differences across seasons and sampling stations in our study. In the warm season, we observed a greater concentration of nutrients (TN, TP, and Chl a) and decreased DO which is consistent with other studies (Lee et al., 2011; Shang et al., 2016). The heavy rainfall in the warm season seems one of the main causes which overload organic matter and pollutants through the estuaries and run-offs. Freshwater inflows and aquaculture activities also deliver nutrients and contaminants besides nutrients from fertilizers and anthropogenic activities, which may also cause the nutrient concentration variations, and shape the bacterioplankton abundance pattern (Jeffries et al., 2016). Generally, low salinity in coastal water indicates that the area has received river discharge, however, in the present study any significant lowering of salinity in the coastal waters was not observed. This perhaps is an indication that Qinhuangdao coastal area does not experience the heavy impact of riverine system round the year. In addition to nutrient inputs, the direct input of pathogenic and fastidious microorganisms from wastewater treatment, wildlife sanctuary, and sewage sources, which are known to influence the bacterioplankton abundance (Vandewalle et al., 2012; Liu et al., 2015) also might affect bacterioplankton abundance and diversity in Qinhuangdao coastal waters.

Our results demonstrated that bacterioplankton abundance in the Qinhuangdao coastal waters exhibit substantial shifts over time and relatively less in space. Furthermore, the strong correlation in total phosphorous and abundance agree well with previous findings that describe inorganic phosphorous as an essential nutrient and that it serves a vital role in cellular energy storage and transformation and is also a growth-limiting nutrient (Brandsma, 2016). Although bacterioplankton and phytoplankton abundance (Chl a) in many studies are closely linked, our study show no relationship between them, an uncoupling phenomenon often observed when inorganic phosphorus is available (Robarts and Carr, 2009). The considerable reduction in bacterioplankton abundance in October perhaps is an indication of anthropogenic inputs or short-term physical forcing events that bring in nutrients which offer phytoplankton an opportunity to outgrow bacterioplankton. The higher abundance of bacterioplankton in warm July and cold March is probably associated with dissolved organic carbon of algal exudates from algal bloom and allochthonous carbon, respectively. Overall, this study extends upon the mechanism of bacterioplankton dynamics and explains the consistency of abundance pattern with the utilization of specific substrates. Nevertheless, the tight connection of abundance with TP bespeaks high implication in coastal pollution monitoring.

Bacterioplankton abundance in the near-shore area is larger than that in the open ocean, especially in the area of the Gulf and the estuarine area. In our study, generally the near-shore abundance was higher than off-shore, perhaps caused by the uneven distribution of organic matter in the coastal stations that are affected directly by run-offs, sewage, human activities, etc. than off-shore stations. As the seawater gets polluted by industrial wastewater and/or domestic sewage, the bacterioplankton may go in a state of unstable nutrition which subsequently affects the total bacterioplankton abundance, transiently. In our study, the abundance was relatively higher in the stations adjacent to river estuaries than that of other stations, while within the near-shore stations, the difference was not significant. Interestingly, in station W5 the bacterioplankton abundance fluctuated dramatically that possibly indicated a strong influence of nearby river estuaries heavily polluted with invasive microorganisms. The results showed that the total abundance in the stations near Qinhuangdao coastal area fluctuated between 8.67 × 104 (W1 in March) to 2.08 × 106 (W4 in October) cells/ml which is in good agreement with the estimate projected by Whitman et al. (Whitman et al., 1998) and empirical data (Al-Rifaie et al., 2008; Robarts and Carr, 2009; Li H.-Y. et al., 2011).

Although total abundance pattern provides an indication of anthropogenic impacts on coastal waters to a certain extent, knowing the pattern of variation in the composition of bacterioplankton allows further understanding of how the growth of different groups of bacterioplankton is controlled by sporadic nutrient inputs. In addition, different phylogenetic groups of bacterioplankton play a specialized role in their consumption of nutrients, and thus dissecting their abundance fluctuations would provide greater insight into the extent of nutrient loading and their recycling by functional groups. The distinct pattern of seasonal variations in the abundance of α-Proteobacteria, Actinobacteria, Firmicutes, and Bacteroidetes groups most likely explains their role in organic matter transformation as chemoheterotrophs. These groups exhibited higher abundance mostly in warm season when there were intense anthropogenic activities and subsequent enrichment of seawater with pollutants. In fact, it was shown that warming promotes bacterioplankton to benefit from the phytoplankton bloom and increase their abundance ultimately shifting toward a more heterotrophic system (Scheibner et al., 2014). These groups perhaps are directly involved in the decomposition of anthropogenic chemicals and algal exudates via a microbial loop as evident from their strong relationship with TN and Chl a. In contrast, β-Proteobacteria was found to be tightly coupled with Chl a, which probably indicated that this less known bacterioplankton group feeds mostly on the algal exudates and may not be directly involved in the decomposition of complex polymers. Bacteroidetes also take part in organic matter degradation and have been associated with the occurrence of algal bloom. Some members of this group use proteorhodopsin to capture light energy for growth making them photoheterotrophic (Gomez-Consarnau et al., 2007). The strong correlation between Bacteroidetes group and temperature could likely be a direct upshot of more light intensity in the warm season.

In our study, temperature, TN, Chl a, and DO were identified as key drivers of spatiotemporal dissimilarity between bacterioplankton abundance patterns in the coastal waters of Qinhuangdao area. Salinity, temperature, and pH are naturally varying parameters, while nutrient concentrations and DO levels indicate the extent of anthropogenic activities besides allochthonous inputs. Temperature and salinity are strong drivers of variations in microbial assemblages of estuaries (Jeffries et al., 2016), but in coastal water salinity changes are marginal, and therefore non-significant in explaining community variations. Previous studies have reported strong association of TN, Chl a, temperature, and TP with the microbial community composition variations and perhaps are the major drivers of community dynamics in aquatic ecosystems (Fodelianakis et al., 2014; Wang et al., 2015; Jeffries et al., 2016; Dong et al., 2017). However, in the present study, these key parameters were tenuous in partitioning the near-shore and off-shore waters near Qinhuangdao area, perhaps indicating the influence of other factors such as Si/dissolved inorganic N ratio, river discharge, short-term physical forcing events, and suspended sediments. In conclusion, this study shows temporal variations in the environmental parameters were more pronounced than spatial in shaping the abundance patterns of bacterioplankton inhabiting the coastal waters of Qinhuangdao. This further connotes the consequence of sporadic anthropogenic activities and their chemical inputs in influencing the bacterioplankton assemblages in the harbor area. Our study provided the first insight into the bacterioplankton ecology of an anthropogenically impacted regional-scale coastal water of Qinhuangdao, by combining bacterioplankton abundance patterns with physicochemical characteristics.

Spatial Patterns of Bacterioplankton Composition and Diversity

The regional and global ocean sampling expeditions have highlighted the importance of bacterioplankton community dynamics and their assembly into functional communities by elucidating their community structure inhabiting different marine locations (Yutin et al., 2007; Lauro et al., 2014). Yet, the information on bacterioplankton diversity within littoral areas adjoining harbors is quite limited perhaps due to complex interactions of multitudinous factors that underline their structure and distribution. As far as the littoral area of Bohai Sea is concerned, this is one of the first studies to our knowledge which reports the diversity of the bacterioplankton assemblage nearby harbor and the factors that shape their structure and composition. Proteobacteria, Cyanobacteria, Actinobacteria, and Bacteroidetes represented the major phylotypes at the 6 stations, and is in agreement with other reports which have found these bacterioplankton groups to be the major inhabitant of marine habitats (Li J. et al., 2011; Li et al., 2013; Smedile et al., 2014; Xiong et al., 2015; Jeffries et al., 2016). Proteobacteria and Cyanobacteria exhibited higher dynamics in their composition among others. Although it is quite challenging to relate phylogenetic diversity and ecological function in marine microbial diversity, yet it is hypothesized that these major groups efficiently metabolize exudates and dissolved organic matter composed of carbohydrates and/or amino acids derived from both growing and senescent phytoplankton (Teira et al., 2008). In addition, these groups are also known to degrade hydrocarbons and inorganic compounds from allochthonous sources (Robarts and Carr, 2009).

Within phylum Proteobacteria, the class α-Proteobacteria represented the majority of the 16S rRNA gene sequences and was the predominant group in all stations. α-proteobacteria are usually abundant in coastal waters and they can form the predominant surface- and particle-colonizing group (Dang and Lovell, 2002a,b). Rhodobacteraceae (also called the marine Roseobacter clade) is usually the largest subgroup in the α-proteobacteria and they are the pioneering and dominant microorganisms on submerged surfaces and organic particles (Dang and Lovell, 2000; Dang and Lovell, 2016). The class γ-Proteobacteria constituted only a minor fraction in our molecular-based study, although members of this class are best known for their rapid growth and also represented as a major group in global surveys (Pommier et al., 2007). OTUs assigned to α-Proteobacteria and Cyanobacteria phyla were cosmopolitans in our study, and this finding is consistent with previous reports (Pommier et al., 2007; Smedile et al., 2014). Cyanobacteria dominated at W1 and W4 stations and most OTUs belonged to the genus Synechococcus. They are mainly free-living and are the main participants and contributors to primary productivity of the global carbon cycle and are the common bacteria inhabiting coastal, temperate, and cold environments (Brown et al., 2014; Xia et al., 2015). From genome analysis, Synechococcus sp. WH8102 was found to be able to use some new organic compounds as sources of nitrogen and phosphorus (Palenik et al., 2003). In addition, Cyanobacteria (Synechococcus) are typically more abundant in higher nutrients or dynamic conditions and are often described as generalist better suited to exploit the fluctuating environments (Palenik et al., 2003). In contrast, Rhodobacteraceae comprise mainly aerobic photo- and chemo-heterotrophs besides purple non-sulfur bacteria which perform photosynthesis in anaerobic environments. Members of this subgroup are metabolically heterogeneous and are capable of anoxygenic phototrophy, aromatic and organosulfur degradation (Lenk et al., 2012). They are the main groups involved in the demethylation of dimethylsulfoniopropionate in the water column (Curson et al., 2011). They are largely involved in sulfur and carbon biogeochemical cycling and symbiosis with other micro- and macro-organisms (Pujalte et al., 2014). The dominance of Rhodobacteraceae in the coastal waters of Qinhuangdao port area epitomizes the degradation of organic and inorganic compounds containing sulfur. The inverse relationship between the two most predominant taxonomic groups (Rhodobacteraceae and Family II) may be either due to competition, niche partitioning, or another mechanism. Further investigation is needed to ascertain the exact cause of the inverse relationship. The low proportion of γ- and δ-Proteobacteria at the 6 stations perhaps indicate the absence of chronic polyaromatic hydrocarbons (PAHs) influence in the coastal waters near Qinhuangdao port, because these groups are often seen associated with high levels of chronic PAHs (Jeanbille et al., 2016). Moreover, the SAR11 considered as highly diverse and abundant at different depths and across habitat types within euphotic and mesopelagic zones was detected only in few stations (W1, W2, and W5). The absence of SAR11 clade at other stations perhaps indicates their rarity and thus this clade was unable to represent in our clone libraries. Nevertheless, the use of inverse Simpson Index which does not tend to be affected by rare species could predict diversity more accurately than other estimators of diversity, and thus diversity analysis in our study was relatively robust and could give a reliable and meaningful interpretation. The wide range and high variation in diversity estimate (1/D) indicate environmental heterogeneity of the Qinhuangdao coastal waters due to varied nutrient types and their concentrations that directly influence the relative abundance of bacterioplankton functional groups.

In the present study, the richness (species number) of the bacterioplankton community was not estimated taking into account the high rate of estimation error with low coverage libraries (Haegeman et al., 2013). Instead, we implemented the estimation of Simpson diversity index which is a robust estimator and provides an accurate estimation of the species diversity. It is independent of any species abundance distributions (SADs) unlike the richness estimators (Haegeman et al., 2013). The use of Simpson index in our study assured accurate estimation of bacterioplankton diversity even with low coverage clone libraries. Also, as explained in an earlier study, coverage does not correlate with library size and the authors clearly mentioned that “a large library size did not guarantee a stable estimate of phylotype richness, nor was a small library necessarily too small” (Kemp and Aller, 2004). The authors also showed that the average size of 220 bacterial libraries was 81 ± 10, ranging in size from 5 to 417 clones. Therefore, our study with 90 clones per library is in good agreement with the average size 81 ± 10 (Kemp and Aller, 2004).

Overall, our study on the bacterioplankton diversity at stations near Qinhuangdao port suggests prominence of only two cosmopolitan OTUs with several endemic OTUs. The predominance of few species is probably an indication of selective advantage of these species to overcome a load of chemical inputs despite scatting grazing and predation by virus and macro-organisms. Also, the effect of anthropogenic and allochthonous inputs was evident from the wide range of inverse Simpson diversity estimates for bacterioplankton assemblages at the stations near Qinhuangdao port. This is consistent with the notion that disturbed marine and coastal environments possess high variations of diversity as a result of complex and unstable nutrient resources, in contrary to oligotrophic marine habitats and deep oceans. The overall findings of this study are expected to serve the basis for Qinhuangdao coastal water monitoring and development of strategies for anthropogenically impacted coastal areas.

Author Contributions

YH has contributed to data acquisition, analysis and drafting the work. BS has contributed to data interpretation, analysis, drafting the work and critical revision. SZ, YZ, and JZ have contributed to data acquisition. NX helped in data analysis. GW contributed to conception of the work, data interpretation and manuscript revision. All authors have final approval to the published version and accountable for all aspects of the work.

Conflict of Interest Statement

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.

Table 1

The data of environmental parameters at near-shore (W1, W3, and W5) and off-shore (W2, W4, and W6) stations near Qinhuangdao coastal area.

Stations	NH4-N (mg/L)	NO2-N (mg/L)	NO3-N (mg/L)	TN (mg/L)	TP (mg/L)	Chl a (mg/m3)
W1	0.067 (0.037)	0.035 (0.024)	0.337 (0.196)	0.793 (0.199)	0.131 (0.080)	11.604 (6.334)
W3	0.061 (0.032)	0.024 (0.017)	0.355 (0.419)	1.103 (0.542)	0.126 (0.077)	5.603 (5.389)
W5	0.104 (0.048)	0.028 (0.018)	0.337 (0.248)	0.959 (0.447)	0.128 (0.019)	16.004 (10.559)
W2	0.053 (0.012)	0.024 (0.017)	0.435 (0.431)	1.158 (0.506)	0.121 (0.079)	3.209 (1.181)
W4	0.041 (0.015)	0.024 (0.022)	0.362 (0.296)	0.851 (0.334)	0.079 (0.054)	2.651 (1.826)
W6	0.054 (0.017)	0.024 (0.019)	0.331 (0.189)	0.937 (0.461)	0.069 (0.063)	5.873 (8.444)

Table 2

Pearson’s correlation analysis showing correlations between abundance of total bacteria, bacterioplankton groups and some environmental parameters.

Bacterial groups	Total nitrogen (mg/L)	Total phosphate (mg/L)	Chlorophyll-a (mg/m3)	Dissolved oxygen (mg/L)	Temperature (°C)
Total bacteria	NS	0.719∗∗(0.00)	NS	NS	NS
α-Proteobacteria	0.490∗(0.015)	NS	0.461∗(0.023)	NS	0.974∗∗(0.00)
β-Proteobacteria	NS	NS	0.532∗∗(0.007)	NS	NS
Cyanobacteria	NS	NS	NS	NS	NS
Actinobacteria	NS	NS	NS	NS	NS
Bacteroidetes	NS	NS	NS	NS	0.410∗(0.046)
Firmicutes	NS	NS	NS	NS	NS

Table 3

The number of 16S rRNA sequences, operational taxonomic units (OTUs), and inverse Simpson diversity (1/D) estimates for the bacterioplankton clone libraries.

Stations	16S rRNA Sequence number	Operational taxonomic units (OTUs)	Inverse Simpson diversity index (1/D)
W1	87	62	34.9
W2	85	62	83.0
W3	80	68	105.0
W4	82	44	16.3
W5	85	47	10.2
W6	89	71	70.3