Interestingly, S5.3, a minor group in marine environments, was widely detected from the bottom layer of stratified stations. S5.3 has at least six clades and shows depth partitioning (Huang et al., 2012). S5.3-I, represented by RCC307, is mainly present in surface water layer, while S5.3-II, -IV, -V, and –VI prevail in the medium to low light layer (Huang et al., 2012). Based on the rpoC1 gene sequence, we found that S5.3 in the Pearl River estuary was not as diverse as in the open ocean and was abundant in the bottom layer (Supplementary Figure S2). Their distribution was significantly positively related to salinity while negatively correlated with temperature, NH4 + and NO3 − (**Figure 5**). This is in agreement with Hashimoto et al.'s (2012) observation that S5.3 mainly occurs in deep waters. Apart from S5.3, we observed that clade I also widely occurred in the bottom layer where temperature could exceed 23◦C. This is in contrast with the conclusion of previous studies that clade I is restricted in high latitude cold water (Zwirglmaier et al., 2008; Huang et al., 2012; Sohm et al., 2016). A recent study reported that clade I contains at least six subclades with different thermal preferences (Xia et al., 2017). Consistently, only warm water subclades, IA and IC (see Figure 8 in Xia et al., 2017), were detected in the Pearl River estuary (Supplementary Figure S3). Besides these two subclades, OTU25, the most abundant clade I OTU, did not cluster with all reported subclades (Xia et al., 2017), but formed another novel subclade (subclade IG) (Supplementary Figure S3). The fact that subclade IG, defined by this study, was mainly distributed in deep water may be the reason why this subclade has not previously been detected. Huang et al. also detected clade I in the South China Sea at relatively deep layers of 75 and 100 m depth with relatively high abundance by sequencing 16S-23S rRNA internal transcribed spacer (ITS) (Huang et al., 2012). This suggests that clade I may be globally distributed and some subclades are specifically distributed in the deep water of tropical/subtropical region.

#### CONCLUSION

The river-estuary-shelf continuum is a highly complex system, which provides a wide array of niches for a highly diverse Synechococcus assemblage ranging from freshwater Synechococcus to euryhaline and strictly marine Synechococcus. Our data suggest that Synechococcus lineages have markedly different abilities to deal with environmental variations. In the estuary, salinity is an important factor influencing the distribution of Synechococcus groups. More studies are needed to reveal the mechanisms involved in salinity tolerance. The fact that high abundance of clade III occurs in the brackish coastal water may revise our previous understanding that clade III prefers oligotrophic oceanic water. Our results further reveal that clade I and S5.3 contain subgroups that have different niches. Further studies should focus on isolation of Synechococcus strains from the studied area and the physiological traits of clades I, III, and S5.3 strains. Moreover, to uncover more details about the distribution of Synechococcus in the salt wedge estuary, high resolution sampling (both vertical and horizontal) need to be conducted in future studies.

#### AUTHOR CONTRIBUTIONS

HL designed the experiment. XX and WG performed the experiments. Data were analyzed by XX in collaboration with WG and HL. XX and HL wrote the manuscript. ST attended the cruise and collected FM and DNA samples. All authors reviewed and approved the final version of the manuscript.

## FUNDING

This study was funded by National Natural Science Foundation of China (NSFC)<PHONE_NUMBER>1), and partially supported by a grant from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. T21/602/16). HL also acknowledge the support of RGC GRF 661813.

## ACKNOWLEDGMENTS

We are grateful to Ms. Candy Lee for analyzing flow cytometric data. Yanping Xu is greatly acknowledged for measuring the nutrients. We also thank Prof. Minhan Dai from Xiamen University for providing us opportunities to collect samples.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2017.01254/full#supplementary-material

## REFERENCES

fmicb-08-01254 July 4, 2017 Time: 16:4 # 11

cyanobacterium Synechococcus sp. WH8102. ISME J. 3, 835–849. doi: 10.1038/ ismej.2009.31

**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.

Copyright © 2017 Xia, Guo, Tan and Liu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

fmicb-08-01254 July 4, 2017 Time: 16:4 # 12

# Patterns and Processes in Marine Microeukaryotic Community Biogeography from Xiamen Coastal Waters and Intertidal Sediments, Southeast China

#### Weidong Chen1,2, Yongbo Pan1 , Lingyu Yu1,2, Jun Yang2 \* and Wenjing Zhang1 \*

1 State Key Laboratory of Marine Environmental Science, Marine Biodiversity and Global Change Research Center, College of Ocean and Earth Sciences, Xiamen University, Xiamen, China, 2 Aquatic EcoHealth Group, Key Laboratory of Urban Environment and Health, Institute of Urban Environment, Chinese Academy of Sciences, Xiamen, China

#### Edited by:

Martin G. Klotz, Washington State University Tri-Cities, United States

#### Reviewed by:

Adélaïde Roguet, University of Wisconsin–Milwaukee, United States Stefan Bertilsson, Uppsala University, Sweden

#### \*Correspondence:

Wenjing Zhang<EMAIL_ADDRESS>Jun Yang<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Aquatic Microbiology, a section of the journal Frontiers in Microbiology

Received: 05 May 2017 Accepted: 20 September 2017 Published: 12 October 2017

#### Citation:

Chen W, Pan Y, Yu L, Yang J and Zhang W (2017) Patterns and Processes in Marine Microeukaryotic Community Biogeography from Xiamen Coastal Waters and Intertidal Sediments, Southeast China. Front. Microbiol. 8:1912. doi: 10.3389/fmicb.2017.01912 Microeukaryotes play key roles in the structure and functioning of marine ecosystems. Little is known about the relative importance of the processes that drive planktonic and benthic microeukaryotic biogeography in subtropical offshore areas. This study compares the microeukaryotic community compositions (MCCs) from offshore waters (n = 12) and intertidal sediments (n = 12) around Xiamen Island, southern China, using high-throughput sequencing of 18S rDNA. This work further quantifies the relative contributions of spatial and environmental variables on the distribution of marine MCCs (including total, dominant, rare and conditionally rare taxa). Our results showed that planktonic and benthic MCCs were significantly different, and the benthic richness (6627 OTUs) was much higher than that for plankton (4044 OTUs) with the same sequencing effort. Further, we found that benthic MCCs exhibited a significant distancedecay relationship, whereas the planktonic communities did not. After removing two unique sites (N2 and N3), however, 72% variation in planktonic community was explained well by stochastic processes. More importantly, both the environmental and spatial factors played significant roles in influencing the biogeography of total and dominant planktonic and benthic microeukaryotic communities, although their relative effects on these community variations were different. However, a high proportion of unexplained variation in the rare taxa (78.1–97.4%) and conditionally rare taxa (49.0– 81.0%) indicated that more complex mechanisms may influence the assembly of the rare subcommunity. These results demonstrate that patterns and processes in marine microeukaryotic community assembly differ among the different habitats (coastal water vs. intertidal sediment) and different communities (total, dominant, rare and conditionally rare microeukaryotes), and provide novel insight on the microeukaryotic biogeography and ecological mechanisms in coastal waters and intertidal sediments at local scale.

Keywords: biogeography, community assembly, plankton, benthos, variation partitioning analysis, neutral model, deterministic processes, stochastic processes

## INTRODUCTION

fmicb-08-01912 October 12, 2017 Time: 17:34 # 2

Microeukaryotes are found in almost all environments on Earth and cover a wide spectrum of cell sizes, shapes and taxonomic affiliations (Schaechter, 2012; de Vargas et al., 2015). Microorganisms, such as algae, protozoa and fungi play a variety of crucial roles in marine ecosystems from primary producers, predators, decomposers to parasites (Sherr et al., 2007). Microeukaryotic community compositions (MCCs) can vary among marine ecosystems with different environmental conditions (Wang et al., 2014; Massana et al., 2015). Sediment and water column form the two most different and important components of marine ecosystems, each with unique environmental conditions and microbial community structure (Feng et al., 2009; Yang et al., 2016). Such differences may lead to different patterns of microbial biogeography. To date, most studies have separately investigated either planktonic or benthic microbial communities (Gong et al., 2015; Yu et al., 2015). Unfortunately, studies comparing planktonic and benthic MCCs, and the relative influence of environmental/spatial factors on the distributions of planktonic and benthic MCCs in marine ecosystems are still extremely limited (Massana et al., 2015).

Revealing the assembly mechanisms that determine the microbial community composition, and how microbes respond to environmental and spatial change, are major challenges in microbial ecology (Hanson et al., 2012; Liu et al., 2015). The growing efforts to understand these mechanisms mainly focus on two major categories (Sloan et al., 2006; Logares et al., 2013; Liu et al., 2015; Wang et al., 2015a; Yang et al., 2016; Liao et al., 2017). The first one is the deterministic mechanism, which considers that environmental filtering determines the biogeographical pattern of microbes (Gilbert et al., 2012). In general, deterministic mechanism predicts that environmental factors will influence microbial community composition (Dumbrell et al., 2010). A large number of studies have confirmed that microbial community composition might be influenced by a variety of environmental variables including nutrients, salinity, pH, temperature and biotic interactions (e.g., predators) (Fernandez-Leborans et al., 2007; Wang et al., 2014, 2015b; Yu et al., 2015). To fully understand the importance of deterministic processes, the sampling scale should be large enough to incorporate multiscale environmental variables (Yeh et al., 2015). The second category is the stochastic mechanism, which suggests that microbial community assembly is influenced by stochastic processes (e.g., drift). Some studies have quantified the importance of stochastic processes using the neutral model (NM) published by Sloan et al. (2006) since it can correctly explain MCC in diverse environments (Logares et al., 2013; Roguet et al., 2015; Burns et al., 2016). Further, the spatial factors (e.g., distance decay), which are part of the stochastic processes, also cause widespread concern in the study of microbial community assembly. However, some literature suggests that both environmental filtering and stochastic processes simultaneously shape microbial biogeography (Liu et al., 2015; Wang et al., 2015c; Yang et al., 2016), therefore both of these two processes need to be considered when investigating the community assembly of plankton or benthos.

Recently, high-throughput sequencing (HTS) has revealed the biogeography of marine rare biosphere community for bacteria or archaea (Galand et al., 2009; Hugoni et al., 2013), which is essential for understanding overall microbial diversity (Telford et al., 2006). Some HTS studies of microbial diversity have provided insight into distribution patterns, and influencing factors, of both abundant and rare microbes in marine and freshwater systems (Logares et al., 2014; Gong et al., 2015; Liu et al., 2015). For example, Liu et al. (2015) revealed that local environmental factors (e.g., water temperature) significantly affected the distribution of rare bacterial taxa in lakes and reservoirs, whereas spatial factors predominately influenced the distribution of abundant taxa (AT). However, it is still unclear how different environmental and spatial factors affect the distribution of planktonic and benthic abundant/rare microeukaryotic taxa in marine ecosystems (Logares et al., 2014; Lynch and Neufeld, 2015). The abundant taxa (AT) and rare taxa (RT) have significantly different richness and community composition, play different ecological functions, and may have discrepant ecological niches (Liu et al., 2015). Thus, they are likely to have distinct responses to environmental and spatial changes (Pedrós-Alió, 2012; Logares et al., 2013), and we hypothesize that abundant and rare microeukaryotes have different assembly mechanisms between water and sediment habitats.

In this study, we compared the community composition of marine planktonic and benthic microeukaryotes, and their interactions with the environmental and spatial factors. Further, we assessed the relative importance of environmental and spatial factors that structure the distribution of microeukaryotic communities (including total, dominant, rare and conditionally rare taxa) around Xiamen Island, southeast China. We aimed to address the following main questions: (1) How different are the MCCs found in the Xiamen coastal waters and intertidal sandy sediments? (2) How are the distribution of planktonic and benthic MCCs (including total, dominant, rare and conditionally rare taxa) influenced by the environmental and spatial variables? (3) Do stochastic processes explain the community variation of plankton and benthos?

## MATERIALS AND METHODS

#### Study Area, Sampling and Environmental Factors

Sampling was conducted during two field cruises around Xiamen offshore area (118◦ 010–118◦ 140 E, 24◦ 250–24◦ 340 N, **Figure 1**). Three biological replicates were selected in water sample stations in the north (N), east (E), south (S), west (W) and sediment stations A, B, C, D to represent of local ecosystem conditions around Xiamen Island (**Figure 1**). The intertidal sandy sediment and coastal water sampling sites were not at the exact same locations, which were ascribed to the intertidal sandy sediment are only distributed in East and South of Xiamen Island. For water samples, the surface (<0.5 m) water samples were transported to the laboratory and processed immediately. For planktonic microeukaryotic community analyses, 800 ml of

surface seawater samples were pre-filtered by a 200 µm poresize sieve to remove debris, large metazoans and grains, and next the water samples with microeukaryotes (smaller than 200 µm) were filtered through 0.22 µm pore polycarbonate membrane (Millipore, Billerica, MA, United States). For sediment samples, the surface (0–5 cm) intertidal sandy sediment samples were transported to the laboratory and processed within 4 h of sampling. The microeukaryotes were mechanically shaken and separated from the sand with sterile seawater over five times, and mixture waters with microeukaryotes were consecutively pre-filtered by a 200 µm pore-size sieve (to remove debris, large metazoans and grains) and filtered/concentrated through 0.22 µm pore polycarbonate membrane (Millipore, Billerica, MA, United States). The membranes were stored at −80◦C until DNA extraction. To limit DNA contamination, the sterile seawater was pre-filtered through a nominal 0.22 µm pore-size membrane (Millipore, Billerica, MA, United States). Filtration equipment was rinsed with sterile water before each sample filtering.

At each sampling site, we measured environmental variables (temperature, pH, salinity, turbidity, chlorophyll-a, and dissolved oxygen) using a Hydrolab DS5 multiparameter water quality meter (Hach Company, Loveland, CO, United States). Total carbon (TC) and total nitrogen (TN) were analyzed with a TOC/TN-VCPH analyzer (Shimadzu, Tokyo, Japan) and total phosphorus (TP) was determined using spectrophotometry according to established standard methods (Wang et al., 2015b). In addition, a Lachat QC8500 Flow Injection Autoanalyzer (Lachat Instruments, Hach Company, Loveland, CO, United States) was used to determine nitrite and nitrate nitrogen (NOX-N), and phosphate phosphorus (PO4-P) concentrations.

#### DNA Extraction, PCR and Illumina Sequencing

Planktonic and benthic DNAs were extracted by E.Z.N.A. DNA Kit (Omega Bio-Tek Inc., Norcross, GA, United States) and FastDNA SPIN kit (MP Biomedicals, Santa Ana, CA, United States) following the manufacturers' instructions,

respectively. The hyper-variable V9 region of eukaryotic 18S rDNA was amplified using the primers 1380F and 1510R (Amaral-Zettler et al., 2009). Each DNA sample was individually PCR-amplified. The 30 µl PCR mixture contained 15 µL of Phusion Master Mix (New England Biolabs, Beverly, MA, United States), 0.2 µM of forward and reverse primers and 10 ng of the sample DNA. PCR reactions including 1 min initial denaturation at 98◦C, followed by 30 cycles of denaturation at 98◦C for 10 s, annealing at 50◦C for 30 s, extension at 72◦C for 60 s. Finally, the amplicons were subjected to 10 min extension at 72◦C. Triplicate PCR products for each sample were conducted and purified using GeneJET Gel Extraction Kit (Thermo Scientific, Hudson, NH, United States). Sequencing libraries were generated using NEB Next Ultra DNA Library Prep Kit for Illumina (New England Biolabs, Beverly, MA, United States) according to manufacturer's instructions, and index codes were added. The library quality was evaluated using the Agilent Bioanalyzer 2100 system (Agilent Technologies, Palo Alto, CA, United States) and Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, United States). The barcoded amplicons from each sample were mixed in equimolar amounts and were then sequenced using an Illumina Miseq platform (Illumina Inc., San Diego, CA, United States) following a paired-end (2 × 250 bp) approach (Caporaso et al., 2012).

#### Bioinformatics

Sequenced paired-end reads were merged with FLASH (Magocˇ and Salzberg, 2011). Raw data were processed and analyzed using the QIIME v.1.8.0 to remove low-quality reads (Caporaso et al., 2010). Sequences were quality controlled using the following settings: maximum number of consecutive low-quality base = 3; minimum of continuous high-quality base = 75% of total read length; ambiguous bases >0 were removed, last quality score = 3 (Liu et al., 2017). Chimeras were identified and removed using UCHIME before the downstream analyses (Edgar et al., 2011). After that, sequences were clustered into OTUs using UPARSE (Edgar, 2013) with the 97% sequence similarity cut-off. Representative sequences in each OTU were selected and blasted against the SILVA 115 reference database (Quast et al., 2013). Unassigned OTUs (sequence similarity to a reference sequence is <80%) and singletons were discarded prior to further analyses, resulting in a mean of 72048 sequences recovered per sample (range from 33996 to 116717). To minimize biases associated with sequencing coverage and allow for comparison between samples, we randomly selected 33996 sequences from each sample.

#### Definition of Abundant and Rare Taxa

We artificially defined thresholds as 0.01% for rare taxa (RT) and 1% for abundant taxa (AT), and classified all OTUs into six categories based on their relative abundance as previously described (Liu et al., 2015; Dai et al., 2016). (1) AT were defined as the OTUs with abundance ≥1% in all samples; (2) RT were defined as the OTUs with abundance <0.01% in all samples; (3) moderate taxa (MT), OTUs with abundance between 0.01 and 1% in all samples; (4) conditionally rare taxa (CRT), with abundance below 1% in all samples and <0.01% in some samples; (5) conditionally abundant taxa (CAT), taxa with abundance ≥0.01% in all samples and ≥1% in some samples but never rare (<0.01%); (6) conditionally rare and abundant taxa (CRAT), OTUs with abundance varying from rare (<0.01%) to abundant (≥1%) (Dai et al., 2016). In this study, we combined AT, CAT, and CRAT as AT to perform further analyses with these three pooled categories taxa (AT, CAT, and CRAT) referred to as 'dominant taxa' to avoid confusion. The ecological studies usually consider 'rarity' to be a continuous variable, therefore, there is always an effect of arbitrariness when defining a cutoff point for rarity (Gaston, 1994). To reduce the effect of arbitrary definition of dominant and rare OTUs, we performed the Multivariate Cutoff Level Analysis (MultiCoLA) to assess the impact of various abundance or rarity cutoff levels on our resulting data set structure and on the consistency of the further ecological interpretation (Gobet et al., 2010). MultiCoLA is an effective method to systematically explore how large community data sets are affected by different definitions of rarity.

## Analyses of Community Diversity

Rarefaction curve, Venn diagram and alpha diversity indices including OTU richness, ACE (abundance based coverage estimator), Chao 1, Shannon-Wiener, Pielou's evenness indices and Simpson index of diversity (1-D) were calculated for each sample or entire samples in vegan 2.4-1 with R software (version 3.2.3) (R Core Team, 2015). Good's coverage was performed in MOTHUR v.1.33.3 software (Schloss et al., 2009). We compared alpha-diversity using one-way ANOVA and Student's t-test.

Bray–Curtis similarity matrix is considered to be one of the most robust similarity coefficients for ecological studies (Kent, 2012) and was applied to our microeukaryotic community dataset. The non-metric multidimensional scaling analysis (NMDS) was employed for detecting possible differences in microeukaryotic planktonic and benthic communities using PRIMER v.7.0 package (PRIMER-E, Plymouth, United Kingdom) (Clarke and Gorley, 2015). Significant difference (P < 0.01) between groups was assessed by analysis of similarities (ANOSIM). The analysis of similarities statistic global R represents separation degree of between-group and within-group mean rank similarities. R = 0 indicates no separation, whereas R = 1 indicates complete separation (Clarke and Gorley, 2015).

#### Relationships between Community Composition, Environmental Variables, and Geographical Distance

To explore the potential controlling factors for the composition of the microeukaryotic community, we used standard and partial Mantel tests to reveal the correlations between the community similarity and environmental factors. Environmental factors, except pH, were square-root transformed and Euclidean distances between samples were calculated. A geographical distance matrix was calculated based on the longitude and latitude coordinates of each sampling site.

Spearman's rank correlation coefficients were calculated to explore the relationships between the Bray–Curtis similarity of microeukaryotic communities and the geographical distance/environmental factors, and the correlation between the geographical distance and the Euclidean distance of all environmental factors.

We quantified the relative effects of environmental and spatial factors in shaping MCC with variation partitioning analysis (VPA) based on redundancy analysis (RDA), as previously described (Wang et al., 2015c). First, a set of spatial variables were generated through the method of principal coordinates of neighbor matrices (PCNM) analysis (Borcard and Legendre, 2002), based on the longitude and latitude coordinates of sampling sites. Subsequently, variance inflation factors (VIFs) were calculated to check the presence of collinearities among environmental and PCNM variables using the "vegan" package and variables with VIF >10 were removed to avoid the impact of collinearity. In addition, to provide unbiased estimates of the variation partitioning based on RDA, microeukaryotic data were Hellinger-transformed prior to the analyses (Legendre and Gallagher, 2001). Forward selection procedure was used to select environmental and spatial variables (Blanchet et al., 2008). Finally, VPA was performed using the "varpart" function of the package vegan. For the seawater habitat, we artificially removed sites N2 and N3 due to their distinct microeukaryotic plankton communities compositions with 10 other sites based on the results of NMDS analyses, and also performed both Mantel tests and VPA using the 10 other water samples.

#### Neutral Community Model for Microeukaryotes

We evaluated the fit of the Sloan neutral community model for MCCs to determine the potential importance of stochastic processes (Sloan et al., 2006, 2007), which considers OTU occurrence frequency in a set of local communities, and their regional relative abundance across the wider metacommunity. The model used here is an adaptation of the neutral theory (Hubbell, 2001) adjusted to be appropriate for large microbial populations. In this model, Nm is an estimate of dispersal between communities. The parameter Nm determines the correlation between occurrence frequency and regional relative abundance, with N describing the metacommunity size and m being immigration rate (Sloan et al., 2006). In this study, we separately used the following data sets – seawater, sediment, and all 24 samples. Further, for seawater habitat, we assessed the fit of the NM for all 12 water samples and 10 water samples (without sites N2 and N3), respectively. All computations were performed in R (version 3.2.3) (R Core Team, 2015).

#### Accession Number

All raw sequences from this study have been submitted to the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) database under the BioProject number PRJNA342297 and the accession number SRP089752.

## RESULTS

#### Comparison of Environmental Factors between Water and Sediment Pore Water

The environmental factors of the studied sites were summarized in Supplementary Table S1. In general, the dynamics of environmental factors from surface waters were more complex than those from intertidal sediment pore waters. The temperature in sediment sites (31.5–33.9◦C) was generally higher than water sites (28.1–30.2◦C) except N2 and N3 (31.5 and 31.7◦C, respectively). Mean values of salinity and pH in the water sites (26.4 psu and 9.35) were higher than sediment sites (25.0 psu and 7.97). However, both TN and TP showed higher mean values in sediment sites (3.07 and 0.112 mg l−1 ) than water sites (0.57 and 0.018 mg l−1 ).

#### Comparison of Richness and Diversity between Plankton and Benthos

Benthic communities showed a significantly higher species richness compared to planktonic communities (Supplementary Figures S1–S3), although Shannon-Wiener, Simpson index of diversity and Pielou's evenness indices did not exhibited significant differences. Rarefaction curves for each sample was seldom saturated (Supplementary Figure S1A). However, the global rarefaction curve indicated that sequencing of ∼0.81 million V9 rDNA sequences from 24 samples was sufficient to approach saturation of microeukaryotic richness (Supplementary Figure S1B). Further, the total number of microeukaryotic OTUs (9003) was roughly equivalent to the number estimated by abundance based richness estimators such as Chao 1 (9005 ± 1.74) and ACE (9026 ± 47) (Supplementary Table S2). Good's coverage ranged from 96.62 to 98.92% in each sample and the index of all 24 samples combined was 99.99%, respectively (Supplementary Table S2). The rarefaction curves, extrapolated species richness indices (Chao 1 and ACE) and Good's coverage indices indicated that the majority of the microeukaryotic taxa had been recovered from the studied metacommunity (Supplementary Figure S1 and Table S2).

A total of 9003 microeukaryotic OTUs were identified from 815 904 high-quality reads at 97% identity level for the 24 samples. Further, Supplementary Table S3 summarized the contribution of each microbial taxa category to microeukaryotic community. MultiCoLA showed that when the structure patterns of community data were compared between the truncated and the original matrices, little variation in data structure was observed after removing up to 60% of the rare types, indicating consistent ecological patterns between the two different matrices. In addition, when the increasing amount of rare types was >5%, consistent ecological patterns were also maintained between the truncated and the original matrices (Supplementary Figure S4). These results indicate that our definitions of dominant taxa (47.54%) and RT (1.88%) were reasonable and objective in this study.

The planktonic microeukaryotic metacommunity generated 4044 OTUs, while the benthic microeukaryotic metacommunity generated 6627 OTUs with the same sequencing depth

(Supplementary Figure S3). It appeared that the absolute number of OTUs was higher in benthic than planktonic metacommunities. The most abundant OTUs that had mean relative abundance >1% represented 29.88% and 21.39% sequences in 12 water and 12 sediment samples, respectively. All of these OTUs were classified as either Stramenopiles or Animalia. In addition, one Arthropoda OTU (OTU\_3) and one Diatomea OTU (OTU\_7452) were present in both water

and sediment habitats with mean relative abundance ≥1%

## Comparison of Microeukaryotic Community Compositions

(Supplementary Table S4).

Overall, the most diverse and abundant OTUs in both water and sediment habitats were assigned to the groups of Opisthokonta (mean relative abundance: 15.64% OTUs and 31.81% sequences for plankton vs. 19.63% OTUs and 40.98% sequences for benthos, respectively) and Stramenopiles (21.13% OTUs and 33.94% sequences for plankton vs. 19.69% OTUs and 27.17% sequences for benthos, respectively) (**Figure 2A**). Some groups such as Alveolata (21.84% OTUs and 12.82% sequences for plankton vs. 14.77% OTUs and 7.11% sequences for benthos, respectively) and Diatomea (14.13% OTUs and 31.88% sequences for plankton vs. 8.37% OTUs and 21.06% sequences for benthos, respectively) exhibited significantly higher abundance and diversity in the water habitat compared with the sediment habitat. However, some groups such as Excavata (1.19% OTUs and 0.33% sequences for plankton vs. 4.28% OTUs and 1.97% sequences for benthos, respectively) exhibited an inverse trend, which were more diverse and abundant in benthic habitat (**Figure 2A**). Further, Venn diagrams indicated that benthic microeukaryotes were more diverse than plankton (including total, dominant taxa, RT, and CRT), since more OTUs were only presented in sediment than in water habitats (**Figure 2B**).

There was a clear difference between planktonic and benthic communities with high global R-values in the ANOSIM tests (**Figure 3**). Interestingly, planktonic communities were further separated into two subgroups, the first subgroup comprised sampling sites N2 and N3, while the second contained 10 other sample sites. Further, the total, dominant, rare and conditionally rare communities showed similar patterns in all NMDS plots.

#### Relationships between Environmental, Spatial Factors and Microeukaryotic Community

The Mantel tests showed that three physicochemical variables (salinity, TP, and TN) were significantly related to the changes of all of the four categories (total, dominant taxa, RT, and CRT) of benthic communities (**Table 1**). Further, both PO4-P and NOX were found to be significant determinants of the total, dominant and CRT subcommunities (P <0.05). Temperature was significantly related to the variation of dominant and CRT for the benthic community. However, the Mantel tests for the planktonic communities (including 12 water sites and 10 water sites) showed that no environmental factor was significantly related to the variation of planktonic community composition.

The Bray–Curtis similarity in the total, dominant, rare and conditionally rare benthic communities between any pair of samples significantly decreased with the increasing of geographical distance (P < 0.01, **Table 1**), indicating that sites in sandy sediments which were geographically distant had less similar communities – presumably because of dispersal limitation. However, no significant distance-decay relationship was found for planktonic communities (including both of 12 sites and 10 sites, **Table 1**). Further, both planktonic and benthic samples that were closer to each other showed more similar environmental conditions (P < 0.01, Supplementary Figure S5). However, benthic sites exhibited a stronger environment-distance relationship than planktonic sites (Spearman's correlation coefficient 0.378 for planktonic sites and 0.839 for benthic sites, respectively).

The VPA indicated that both environmental and spatial factors play significant roles in structuring microeukaryotic metacommunity in water and sediment habitats (**Table 2**). For total and dominant taxa, the spatial factors exhibited a slightly higher contribution to the community variation than environmental factors, however, a high proportion of unexplained variation was observed in the RT (78.1–97.4%) and CRT (49.0–81.0%) (**Table 2**). For total MCC, the variance explained by purely spatial factors was slightly higher (27.1% vs. 24.2% for plankton and 30.2% vs. 19.1% for benthos, respectively) than that of purely environmental factors in both habitats. For the dominant subcommunity, spatial factors alone also explained higher community variations (28.5% vs. 23.2% for plankton and 33.0% vs. 18.5% for benthos, respectively) than the environmental factors. Interestingly, for the planktonic samples after removing sites N2 and N3, spatial factors exhibited a much higher contribution to the community variations than environmental factors among all of total, dominant taxa, RT, and CRT. Moreover, the contribution of environmental factors to the planktonic MCC variations exhibited a rapid decline among all of total (3.4%), dominant (11.0%), rare (0%), and conditionally rare (4.2%) planktonic taxa. For the benthic microeukaryotic community, the unexplained variation of the RT (78.1%) and CRT (49.0%) was higher than that of the total (46.3%) and dominant (40.4%) communities. Similarity, in the water sites (12 sites and 10 sites, respectively), the unexplained variation of the RT (87.1 and 97.4%, respectively) and CRT (53.3 and 81.0%, respectively) was also higher than that of the total (44.6 and 76.2%, respectively) and dominant (43.2 and 60.3%, respectively) communities.

#### Fit to the Neutral Model of Community Assembly

The NM estimated a low fraction of the variation in the frequency of occurrence of different OTUs in planktonic communities (**Figure 4A**, R 2 = 0.25). The benthic microeukaryotic community showed a larger fraction of the variation in the frequency of occurrence than planktonic microeukaryotic community (**Figure 4B**, R 2 = 0.43), indicating a medium fit to the NM, although the Nm-value was higher for planktonic communities (Nm = 10103) than for benthic communities (Nm = 9338). We

and abundant taxa; RT, rare taxa; CRT, conditionally rare taxa.

also tested the fit to the NM for 10 planktonic communities excluding sites N2 and N3, and found the highest Nmvalue (Nm = 34403). Interestingly, this model explained the largest proportion of the variability in occurrence frequency of planktonic community without sites N2 and N3 among the all different data-sets (**Figure 4C**, R 2 = 0.72). Further, we evaluated the fit to neutral community model on all 24 planktonic and benthic communities, and found that the observed data exhibited no fit to the neutral curve (Supplementary Figure S6, R 2 = −0.09, note that negative R 2 -values can occur when there is no fit to the model), indicating a strong role for habitat filtering relative to stochastic processes.

## DISCUSSION

## Microeukaryotic Composition and Diversity Patterns in Water and Sediment

Our DNA metabarcoding results validated a well-known observation that planktonic and benthic microeukaryotic communities composition and structure are significantly different (Zinger et al., 2011; Massana et al., 2015; Forster et al., 2016), which was evident in the alpha-diversity (e.g., OTU richness and ACE index), beta-diversity and taxonomic composition (e.g., Diatomea) (**Figures 2**, **3** and Supplementary Figures S2, S3). The comparison of MCCs between water and sediment habitats provided additional insight into the diversity and distribution of subtropical marine microeukaryotes.

The microeukaryotic community differences between planktonic and benthic taxa held at both supergroup and OTU levels, although 18.53% OTUs were common to both habitats (**Figure 2**). First, the dominant overlap between benthic and planktonic taxa in our data was related to Alveolata, Stramenopiles (represented mainly by Diatomea) and Opisthokonta (represented mainly by Arthropoda). Although a 200 µm pore-size sieve was used to remove large metazoans, we found that the Opisthokonta was an abundant group, and they were dominated by Arthropoda, which accounted for 23.62% of planktonic sequences and 11.76% of benthic sequences. One explanation for this is that some large-sized organisms can get through the 200-µm pores and/or presence of eggs, spores, or larvae of the large-sized organisms (Liu et al., 2017). Another alternative explanation for this is that our DNA-based molecular approach cannot exclude the potential effect by "free DNA in the water" or "animal debris" passing the 200 micron filter on the perceived microeukaryote communities (Thomsen and Willerslev, 2015; Liu et al., 2017). Both Alveolata and Stramenopiles are dominant benthic microeukaryotic groups, as is also shown by previous studies (Gong et al., 2015; Forster et al., 2016). However, our finding showed that Alveolata was more dominant in water habitat. Further, another study (Mann and Evans, 2007) showed that Diatomea were more often detected in the benthos compared with plankton among phototrophic protists, which is in contrast to our data. This difference might be ascribed to different sampling areas in these various studies. Excavata is a major eukaryotic supergroup, however, only a few studies have focused on it (Hampl et al., 2009) compared to other groups (e.g., Alveolata and Diatomea). Second, our study showed that benthic microbial communities exhibited higher richness than planktonic communities although there was a little overlap between the two communities (**Figure 2B**). Studies in the Gulf of Maine and Long Island Sound in North America also revealed little overlap between genetic signatures of the water and sediment microeukaryote assemblages (Doherty et al., 2010). One explanation for this observation might be that benthic habitats exhibited stronger horizontal and vertical gradients

TABLE 1 | Mantel tests for the correlations between geographical distance, environmental distance, individual environmental factors, and microeukaryotic community similarity in water and sediment using Spearman's coefficients.

Geo\_distance, pairwise geographical distances between sampling sites; Env\_distance, Euclidean distance of all environmental variables between sampling sites; Community similarity was based on the Bray–Curtis similarity. The significances are tested based on 999 permutations, and bold values indicate statistical significance (P < 0.05). Water samples with 10 sites indicated removing sites N2 and N3. All, all taxa; AT, abundant taxa; CAT, conditionally abundant taxa; CRAT, conditionally rare and abundant taxa; RT, rare taxa; CRT, conditionally rare taxa. TN, total nitrogen; TP, total phosphorus; DO, dissolved oxygen; TC, total carbon; NOX-N, nitrate and nitrite nitrogen; PO4-P, phosphate phosphorus.

in both physical and chemical characteristics. Environmental heterogeneity is much more pronounced in the benthic than in the planktonic sites, and higher environmental heterogeneity likely promotes the existence of highly specialized organisms, thereby probably driving species-richness patterns (Hortal et al., 2009).

Although further investigations are necessary to characterize the effects of DNA extraction protocols and PCR methods on MCCs, the robust patterns in microeukaryotic community were distinct between coastal waters and intertidal sediments. Recently, a study of amplicon-based characterization of microbial community structure revealed that using different DNA extraction kits did not significantly affect alpha diversity/richness, beta diversity and dominant members of the community of samples (Staley et al., 2015). Therefore, we consider that different DNA extraction kits did not significantly affect the overall biological conclusions drawn in this study. Our NMDS analysis showed that communities from the same habitat cluster together (**Figure 3**), which indicated a distinct composition among planktonic and benthic communities. Further, in the water habitat, the differences between sites N2 and N3 compared to the 10 other water sites have been showed by a previous study (Yu et al., 2015).

#### Controlling Factors Shaping the Microeukaryotic Community Compositions

The distinct distribution patterns of MCCs between water and sediment could be attributed to two types of mechanisms: environmental filtering and stochastic processes (Liu et al., 2015; Yang et al., 2016; Liao et al., 2017).

TABLE 2 | Variation partitioning of the microeukaryotic communities among environmental and spatial variables in water and sediment.

Environment, community variation explained by pure environmental factors; Spatial, community variation explained by pure spatial factors; Shared, community variation explained by joint effect of environmental and spatial factors. Bold values indicate statistical significance (P < 0.05). Negative value of adjusted coefficients of determination (Adjusted R2 ) was converted to zero. Water (10 sites) indicated removing sites N2 and N3.

First, environmental factors such as nutrient concentrations can influence MCCs because they are essential for the growth and development of microorganisms, and different microorganisms are adapted to their optimal growth concentrations (Kneip et al., 2007; Wang et al., 2015b). MCCs can be both directly and indirectly affected by nutrient concentrations, which can influence the photosynthesis of autotrophs, and heterotrophic microeukaryotes which can prey on autotrophs (Hecky and Kilham, 1988). In this study, the Mantel tests revealed that benthic MCCs (total, dominant taxa, RT, and CRT) were significantly correlated with salinity, TP and TN in sediments (**Table 1**). These results are consistent with studies of bacterioplanktonic or microeukaryotic communities responses to TP or TN concentrations (Liu et al., 2011; Wang et al., 2015a). Further, salinity is also an important environmental factor affecting the survival and growth of marine microeukaryotes. Changes in salinity directly affect the osmoregulation and metabolism of microorganism, leading to effects on the life activities of microorganism (Decamp et al., 2003). However, our results indicate that there is a clear difference in correlation between plankton and benthos with environmental factors (**Table 1**). Interestingly, no measured environmental variable was significantly associated with planktonic MCC (including 12 water sites and 10 water sites). There are at least two reasons that might possibly explain the observed non-significant relationship between plankton and environmental variables. The first possibility is that plankton communities are extremely dynamic, and snapshot or even regular but low resolution sampling includes a large fraction of noise, which may mask main ecological patterns (Özkan et al., 2014). Second, variation in planktonic communities not explained by environmental (e.g., abiotic) controls might indicate strong biotic interactions (Wang et al., 2015b). In addition, there are some potential limitations that merit further discussion. Some unmeasured environmental

variables, such as irradiance and ocean currents, are important in shaping the distributions of MCCs (Yeh et al., 2015). The low statistical power for detecting patterns was perhaps largely due to small sample size in field metacommunity studies, and the sampling scale should be considered when designing sampling schemes (Martiny et al., 2011). That is, the sampling size should be large enough to include most important environmental gradients and dispersal factors.

Second, an increasing body of literature in microbial community ecology indicates that spatial factors are another major process shaping biogeographical patterns besides environmental selection (Vyverman et al., 2007; Cermeno and Falkowski, 2009; Hanson et al., 2012). Our results showed that spatial effects were significant in shaping benthic MCCs (total, dominant taxa, RT, and CRT), whereas marine planktonic MCCs (including 12 water sites and 10 water sites) do not adhere to significant distance-decay relationship (**Table 1**). Recently, Gong et al. (2015) surveyed microeukaryotes of coastal sediment sites in the Yellow Sea and revealed limited dispersal of benthic microeukaryotes, which showed a similar distribution pattern to our results. Comparable and similar results have been noted for benthic foraminifera in deep sea and benthic diatoms in lakes (Telford et al., 2006). However, marine planktonic picoeukaryote such as MAST-4 did not associate with significant distance-decay pattern, but the environmental selection (e.g., temperature) drove MAST-4 distribution (Rodríguez-Martínez et al., 2013). Another study on marine planktonic diatom assemblages also showed that these eukaryotic microbes were not limited by dispersal (Cermeno and Falkowski, 2009). Overall, these results suggest that benthic and planktonic microeukaryotes might have fundamentally different biogeographical patterns or ecological mechanisms. There are large differences in benthos compared to plankton including life environment, food-web and body structure (Schaechter, 2012). In sediments, movements of microeukaryotes are obviously more limited relative to water columns, where planktonic microbes can disperse passively to distant locations because of ocean currents (Yeh et al., 2015). These differences can lead to distinct biogeographical patterns between plankton and benthos at a regional scale.

More importantly, our results show that the rare and conditionally rare benthic subcommunities follow the similar and general biogeographical distribution patterns with dominant taxa, because all of the three subcommunities exhibited significant distance-decay patterns (**Table 1**). This could be because of dispersal limitation and/or the fact that benthic samples probably have similar environmental conditions, if they are close to each other (Supplementary Figure S5; Hanson et al., 2012). Recent study of bacteria using HTS in the lakes and reservoirs indicated that rare bacteria showed a similar spatial pattern with the abundant bacterial subcommunity (Liu et al., 2015). These results indicate that similar structuring processes (e.g., dispersal) can influence benthic dominant, rare and conditionally rare subcommunity compositions (Hanson et al., 2012). However, although there were similar distribution patterns of dominant taxa, RT, and CRT, we found that dominant taxa have a weaker distance-decay relationship (r = −0.400) than RT and CRT (r = −0.579 and r = −0.798, respectively, **Table 1**). This suggests that the benthic dominant microeukaryotic community has a higher capacity to dispersal than RT, hence leading to a more cosmopolitan distribution (Liu et al., 2015).

#### Relative Importance of Environmental and Spatial Factors Influencing Microeukaryotic Distribution

Our results indicate that both the environmental and spatial factors play significant roles in structuring the microeukaryotic metacommunity, although the relative effects of the two types of influencing factors on community variations were different (**Table 2**). The mechanisms underlying the assembly processes are important question in the field of microbial ecology (Hanson et al., 2012). The different relative effects of environmental and spatial factors influencing MCCs may be ascribed to the dominant taxa, RT, and CRT which had different community compositions and properties (Liu et al., 2017). In aquatic ecosystems, the responses of microbial communities to environmental and spatial changes are mediated by their properties; such as physiological tolerance, dispersal capacity, taxonomic and functional diversity (Gong et al., 2015; Liu et al., 2015). Our result was consistent with previous study on aquatic and sedimentary bacteria communities in 16 lakes from western China (Yang et al., 2016). However, our study was not in agreement with another study in Antarctica coastal lakes (Logares et al., 2013), which suggested that MCC was strongly influenced by environmental factors (e.g., salinity) and weakly correlated with geographical distance. This inconsistency maybe attributed to the spatial distance differences among the studied area and different environmental gradients (e.g., salinity ranged from 0 to 100 psu in coastal lakes in Antarctica vs. 12.6 to 30.9 psu in water habitat and 23.2 to 30.0 psu in sediment habitat in Xiamen) (Logares et al., 2013). In addition, our VPA confirmed that a large proportion of the observed MCCs variations could not be explained, especially for the RT (78.1–97.4% of variance unexplained) and CRT (49.0–81.0%), indicating that more complex mechanisms may generate and maintain the rare biosphere diversity in the coastal waters and intertidal sandy sediments. The unexplained variation may be due to unmeasured environmental and ecological factors, or interactions between taxa and species' own vital rate. For example, Barberán and Casamayor (2010) found that species sorting solely dominated microbial community composition. Evolutionary drift (stochastic genetic diversification), and ecological drift (stochastic processes of birth, death, colonization, and extinction) could also contribute to spatial effects and unexplained variation (Hanson et al., 2012). Further, another studies indicated that the co-occurrence correlations among microbes should also be responsible for the community structure (Lima-Mendez et al., 2015; Wei et al., 2016). These findings suggested that extending the number of sampling sites, and integrating more environmental factors, along with conducting manipulative experiments across space and time, are necessary for better understanding of the influence and mechanisms of environmental selection and stochastic processes on the biogeography of planktonic and benthic microeukaryotes

(Yeh et al., 2015). Therefore, the relative importance of spatial and environmental factors on microeukaryotic distribution awaits further investigation given the unexplained variation.

After removing sites N2 and N3, interestingly, the contribution of environmental factors to the planktonic MCC variations exhibited a rapid decline among all of total, dominant, rare, and conditionally rare planktonic taxa, whereas the contribution of spatial factors did not show significant variation. This result confirmed that environmental factors had a stronger influence on microbial biogeography than spatial factors in these 12 water sites than for the 10 water sites in this study (**Table 2**). More details of sites N2 and N3 were reported in the previous study (Yu et al., 2015).

Finally, we fitted the NM of community assembly for determining whether stochastic processes could explain variation of MCCs (**Figure 4** and Supplementary Figure S6). The NM was very powerful in explaining microbial community structure (R 2 = 0.76) within aquatic environments in a previous study (Roguet et al., 2015). In this study, the fit to the NM showed that stochastic processes appeared to have more influence in total benthic MCC (R 2 = 0.43, 12 sediment sites) compared with total planktonic MCC (R 2 = 0.25, 12 water sites). After removing sites N2 and N3, however, this model explained 72% of the variation in planktonic MCC (R 2 = 0.72, 10 water sites, **Figure 4**). These results indicate that stochastic processes contributed the most impact in planktonic MCC without sites N2 and N3 (Sloan et al., 2006), and this finding may be ascribed to high and random dispersal rate of microeukaryotes in surface water. Their movements are more restricted in sediment than in water habitat, where they can disperse randomly or passively with ocean currents at regional or local scales. The NM findings (the R 2 -value of 12 water sites was 0.25 and 12 sediment sites was 0.43, respectively) were almost fully consistent with results revealed by VPA (the community variation explained by spatial factors for plankton and benthos), although the contributions of spatial factors to the MCC variations differed between the two statistical methods. The difference may be attributed to the MCC variation explained by VPA only involve the influence caused by geographical distance of sampling sites (dispersal). However, various stochastic-related processes (drift and other mechanisms) could contribute major influence to the MCC variation in NM (Sloan et al., 2006, 2007).

#### CONCLUSION

Our data demonstrated that MCCs were significantly different between water and sediment habitats, and benthos exhibited significantly higher species richness than plankton largely due to high proportions of RT and CRT. The proportions of OTUs and sequences of Excavata were significantly higher in the benthos than in the plankton, whereas Alveolata and Diatomea had significantly higher richness and abundance in planktonic than benthic communities. Mantel tests revealed that environmental factors (e.g., salinity, TP, and TN) were significantly related to the variation of benthic MCCs. However, the lack of significant relationship between environmental factors and planktonic MCCs might indicate plankton communities were affected by more complex mechanisms (such as biotic interactions, stochastic processes and other unmeasured environmental variables). Further, benthic MCCs (including total, dominant, rare, and conditionally rare taxa) exhibited a significant distancedecay relationship, whereas the planktonic communities did not. However, more than 70% variation in planktonic community was well explained by stochastic processes when removing two unique sites (N2 and N3). Both the environmental and spatial factors play significant roles in structuring total and dominant microeukaryotic metacommunities although their relative effects on the planktonic and benthic community variations were different. However, unlike the higher explained variation of the total and dominant taxa, the high proportion of unexplained variation in the RT and CRT indicated that more complex mechanisms may influence the rare subcommunity assembly in the coastal waters and intertidal sediments. Altogether, our results indicate that patterns and processes in marine microeukaryotic community assembly may differ among the different habitats (coastal water and intertidal sediment) and different MCCs (total, dominant, rare and conditionally rare microeukaryotes).

## AUTHOR CONTRIBUTIONS

WZ conceived the idea and designed the comparison experiments. YP and LY collected the samples and performed the experiments. WC, JY, and WZ analyzed the data. JY and WZ contributed reagents/materials/analysis tools. WC, JY, and WZ wrote the paper.

#### FUNDING

This work was supported by the National Natural Science Foundation of China (41276133), the Natural Science Foundation of Fujian Province (2014J01163) and the Fundamental Research Funds for the Central Universities<PHONE_NUMBER>5).

#### ACKNOWLEDGMENTS

The authors thank Lemian Liu, Yuanyuan Xue, and Min Liu for data analysis assistance, and thank Prof. David M. Wilkinson and Dr. Alain Isabwe for constructive comments on the manuscript. The authors also thank two reviewers for insightful comments that helped in improving the clarity of this paper.

## SUPPLEMENTARY MATERIAL

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

## REFERENCES

fmicb-08-01912 October 12, 2017 Time: 17:34 # 13

**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.

Copyright © 2017 Chen, Pan, Yu, Yang and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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

#### Yaodong He1† , Biswarup Sen1† , Shuangyan Zhou1 , Ningdong Xie1 , Yongfeng Zhang2 , Jianle Zhang2 and Guangyi Wang1,3 \*

1 Center for Marine Environmental Ecology, School of Environmental Science and Engineering, Tianjin University, Tianjin, China, 2 Qinhuangdao Marine Environmental Monitoring Central Station, State Oceanic Administration, Qinhuangdao, China, 3 Key Laboratory of Systems Bioengineering (Ministry of Education), Tianjin University, Tianjin, China

#### Edited by:

Hongyue Dang, Xiamen University, China

#### Reviewed by:

Mohd Ikram Ansari, King Abdullah University of Science and Technology, Saudi Arabia Jérôme Hamelin, Institut National de la Recherche Agronomique (INRA), France

\*Correspondence:

Guangyi Wang<EMAIL_ADDRESS>†These authors have contributed equally to this work.

#### Specialty section:

This article was submitted to Aquatic Microbiology, a section of the journal Frontiers in Microbiology

Received: 17 February 2017 Accepted: 03 August 2017 Published: 18 August 2017

#### Citation:

He Y, Sen B, Zhou S, Xie N, Zhang Y, Zhang J and Wang G (2017) Distinct Seasonal Patterns of Bacterioplankton Abundance and Dominance of Phyla α-Proteobacteria and Cyanobacteria in Qinhuangdao Coastal Waters Off the Bohai Sea. Front. Microbiol. 8:1579. doi: 10.3389/fmicb.2017.01579 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

**89**

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.

Keywords: anthropogenic impacts, environmental variables, phylogenetic diversity, bacterioplankton abundance, seasonal variations, redundancy analysis

#### INTRODUCTION

fmicb-08-01579 August 16, 2017 Time: 15:58 # 2

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).   Overall, the available reports on bacterioplankton community composition at regionalscale 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 1**). 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 offshore 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 1**). 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.

## 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

Frontiers in Microbiology | www.frontiersin.org

universal primers 27F (50 -AGAGTTTGATCCTGGCTCAG-30 ) and 1492R (50 -GGTTACCTTGTTACGACTT-30 ) 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

fmicb-08-01579 August 16, 2017 Time: 15:58 # 4

**92**

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 R 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 1**. The NH4-N level was high in W5 station than

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.

TN, TP, and Chl a indicate total nitrogen, total phosphate, and Chlorophyll a.

Mean (SD) of values over 4 seasons (July, October, December, and March) are provided.

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 2**). 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.

#### 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 3**. 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 4**).

The qPCR based abundance estimates of 6 major phylum/class i.e., α-Proteobacteria, β-Proteobacteria, Cyanobacteria, Actinobacteria, Firmicutes, and Bacteroidetes are shown in **Figure 5**. α-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.

**Table 2** 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 6**) 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.

## 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 3** and **Figure 7**. 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 3**). A broad phylogenetic distribution of the bacterioplankton community at the 6 stations in warm season is shown in **Figure 7A**. The analysis of 16S rRNA gene sequences from these stations revealed 7 distinct phyla (**Figure 7A**). 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

FIGURE 3 | 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.

(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 7B**). 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).

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 3**). 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

∗∗P < 0.01, ∗P < 0.05; NS, non-significant.

The values represents correlation coefficient (significant bilateral) for sample number N = 24. Total bacteria and bacterioplankton group abundance data were based on fluorescence microscopy and quantitative PCR, respectively.

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.

(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

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

W1–W6 represents the 6 stations near Qinhuangdao coastal area for which bacterioplankton clone libraries were constructed and analyzed.

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 nearshore 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

#### REFERENCES

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.

#### ACKNOWLEDGMENTS

The work was partially supported by NSFC (grant # 31670044 and 31602185) and National Marine Public Welfare Industry Special Scientific Research Project (201305022). The views expressed herein are those of the authors and do not represent the views of the funding agencies or any of its subagencies.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fmicb. 2017.01579/full#supplementary-material

**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.

Copyright © 2017 He, Sen, Zhou, Xie, Zhang, Zhang and Wang. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Long-Term Survey Is Necessary to Reveal Various Shifts of Microbial Composition in Corals

Shan-Hua Yang1† , Ching-Hung Tseng2† , Chang-Rung Huang1‡ , Chung-Pin Chen3 , Kshitij Tandon1,4,5, Sonny T. M. Lee6 , Pei-Wen Chiang1 , Jia-Ho Shiu1,7,8, Chaolun A. Chen1 and Sen-Lin Tang1 \*

1 Biodiversity Research Center, Academia Sinica, Taipei, Taiwan, 2 Germark Biotechnology Co., Ltd., Taichung, Taiwan, 3 Yourgene Bioscience, New Taipei City, Taiwan, 4 Bioinformatics Program, Institute of Information Science, Taiwan International Graduate Program, Academia Sinica, Taipei, Taiwan, 5 Institute of Bioinformatics and Structural Biology, National Tsing Hua University, Hsinchu, Taiwan, 6 Section of Gastroenterology, Hepatology and Nutrition, Department of Medicine, University of Chicago Medicine, Chicago, IL, United States, 7 Molecular and Biological Agricultural Sciences Program, Taiwan International Graduate Program, Academia Sinica, Taipei, Taiwan, 8 Graduate Institute of Biotechnology, National Chung Hsing University, Taichung, Taiwan

#### Edited by:

Hongyue Dang, Xiamen University, China

#### Reviewed by:

Luisa I. Falcon, National Autonomous University of Mexico, Mexico Fabiano Thompson, Federal University of Rio de Janeiro, Brazil

#### \*Correspondence:

Sen-Lin Tang<EMAIL_ADDRESS>†These authors have contributed equally to this work and have shared co-first authorship.

‡Deceased 17 June, 2016

#### Specialty section:

This article was submitted to Aquatic Microbiology, a section of the journal Frontiers in Microbiology

Received: 06 March 2017 Accepted: 30 May 2017 Published: 13 June 2017

#### Citation:

Yang S-H, Tseng C-H, Huang C-R, Chen C-P, Tandon K, Lee STM, Chiang P-W, Shiu J-H, Chen CA and Tang S-L (2017) Long-Term Survey Is Necessary to Reveal Various Shifts of Microbial Composition in Corals. Front. Microbiol. 8:1094. doi: 10.3389/fmicb.2017.01094 The coral holobiont is the assemblage of coral host and its microbial symbionts, which functions as a unit and is responsive to host species and environmental factors. Although monitoring surveys have been done to determine bacteria associated with coral, none have persisted for >1 year. Therefore, potential variations in minor or dominant community members that occur over extended intervals have not been characterized. In this study, 16S rRNA gene amplicon pyrosequencing was used to investigate the relationship between bacterial communities in healthy Stylophora pistillata in tropical and subtropical Taiwan over 2 years, apparently one of the longest surveys of coral-associated microbes. Dominant bacterial genera in S. pistillata had disparate changes in different geographical setups, whereas the constitution of minor bacteria fluctuated in abundance over time. We concluded that dominant bacteria (Acinetobacter, Propionibacterium, and Pseudomonas) were stable in composition, regardless of seasonal and geographical variations, whereas Endozoicomonas had a geographical preference. In addition, by combining current data with previous studies, we concluded that a minor bacteria symbiont, Ralstonia, was a keystone species in coral. Finally, we concluded that long-term surveys for coral microbial communities were necessary to detect compositional shifts, especially for minor bacterial members in corals.

Keywords: long-term survey, microbial composition, corals, Stylophora pistillata, shifts

## INTRODUCTION

Corals harbor a variety of microorganisms, including algae, fungi, bacteria, archaea, and viruses (Knowlton and Rohwer, 2003; Wegley et al., 2007; Marhaver et al., 2008). Coral-associated microbes, which live in coral mucus, tissue, and skeleton, have various interactions with their host and maintain coral holobiont function, including antimicrobial defense and nutrient acquisition (Bourne et al., 2016). With advances in high-throughput sequencing technologies, there is growing

**104**

research interest in microbial ecology in rare biospheres (Lynch and Neufeld, 2015), characterizing genetic and functional diversity supporting global ecosystem health (Yachi and Loreau, 1999; Elshahed et al., 2008).

 According to previous studies, some microbes had host specificity, regardless of geographical difference (Rohwer et al., 2001, 2002), whereas others had regional specificity rather than host types (Littman et al., 2009).

Various coral-associated microbial species reportedly shared similar functional roles in coral ecology (Li et al., 2013). However, bacterial succession in corals has not been well characterized, as the longest monitoring survey lasted only 1 year (Rohwer et al., 2002; Hong et al., 2009; Ceh et al., 2011; Carlos et al., 2013; Lema et al., 2014; Hernandez-Agreda et al., 2016). Furthermore, although some microbes were consistently detected in coral (Rohwer et al., 2002; Carlos et al., 2013), it is still unknown whether these common microbes change their abundance over longer intervals. A 2-years survey characterizing dynamics of Symbiodinium rare biosphere provided evidence for symbiont switching in reef-building corals after environmental changes (Boulotte et al., 2016), but dynamics of the bacterial rare biosphere in coral, and bacterial shifting/switching within corals, are not established.

The objective was to determine whether a long-term survey was necessary to detect compositional dynamics of coral-associated bacteria in various locations. Based on previous short-term studies, we hypothesized that coral-associated bacteria require varying time scales to manifest detectable variations in different environments. Cumulative stressors, e.g., nutrient availability and elevated sea surface temperature, have been postulated to change coral reef ecosystems; therefore, coral-associated microbes might either attenuate or exacerbate effects from stressors through their roles in coral holobiont fitness (Bourne et al., 2008, 2016). Since cumulative stressors vary from place to place, duration of detectable changes in coral-associated microbes is also likely to vary, which should be apparent in a long-term survey (Mouchka et al., 2010; Ceh et al., 2011; Chen et al., 2011; Lee et al., 2012).

To determine dynamics of coral-associated microbial communities among locations and to study temporal variations of dominant and minor bacterial groups, 16S rRNA gene amplicon pyrosequencing was used in a 2-years observational study investigating bacterial communities in coral Stylophora pistillata, which is widely distributed in Taiwan and around the globe. That Taiwan lies across the Tropic of Cancer, S. pistillata in tropical (Kenting and Lyudao) and subtropical (Yehliu) Taiwan were collected to provide a comparative basis for latitude influence on bacterial composition in corals. This 2-years survey was to assess coral-associated microbial differences among various climate zones, and determine importance of long-term observation to detect bacterial dynamics in coral hosts.

## MATERIALS AND METHODS

## Coral Samples

Coral and seawater samples were collected from three locations, Kenting, Yehliu, and Lyudao, in Taiwan (**Figure 1**) from February 2008 to January 2010. In each location and at each sampling time, three colonies (replicates) of S. pistillata at depths of 3–7 m were selected and tagged for sampling. For each colony, a coral fragment approximately 2 cm long and 2 cm in diameter was collected. Samples were washed twice with filtered seawater and stored at 4◦C until coral tissue was processed. In addition, 1 L of seawater surrounding the coral colony was simultaneously collected to analyze the microbial community in reef water. A total of 88 S. pistillata and seawater samples were collected (Supplementary Table S1). Environmental data (precipitation, daylight duration and air and water temperatures) were obtained from the Central Weather Bureau in Taiwan (Supplementary Table S1).

Coral tissues were sprayed (using an airbrush) with 2–5 mL 10× TE buffer (10 mM Tris–HCl, pH 8.5, 1 mM EDTA, pH 8.0) and collected in a sterilized bag. Coral tissues were transferred to a 15 mL sterilized tube and stored at −20◦C pending DNA extraction.

## DNA Extraction, Primers, PCR, and Sequencing

Coral tissue was recovered for genomic DNA extraction by centrifugation (10,000 × g, 15 min). After discarding the supernatant, tissue pellets were put into 600 µL of 10× TE buffer and liquid nitrogen, homogenized (mortar and pestle), and the resulting solution transferred to a clean tube. Seawater samples were filtered through membrane filters (0.2 µm pore size; Adventec, Tokyo, Japan) to collect microbial particulate. Total genomic DNA of coral tissue and seawater were extracted as outlined in the UltraClean Soil DNA Kit (MioBio, Solana Beach, CA, United States).

The PCR was performed using two bacterial universal primers: 27F (50 -AGAGTTTGATCMTGGCTCAG-30 ) and 341R (50 -CTGCTGCCTCCCGTAGG-30 ), which were designed for the V1–V2 hypervariable region of the bacterial 16S rRNA gene (Hamp et al., 2009). The reaction mixture contained 1 µL of 5 U TaKaRa Ex Taq HS (Takara Bio, Otsu, Japan), 5 µL of 10× Ex Taq buffer, 4 µL of 2.5 mM deoxynucleotide triphosphate mixture, 1 µL of each primer (10 µM), and 1–5 µL (10–20 ng) template DNA in a volume of 50 µL. The temperature program for 30 PCR cycles was the initial step of 94◦C for 3 min, 94◦C for 20 s, 52◦C for 20 s, 72◦C for 20 s, and 72◦C for 2 min (final extension after the last cycle).

The PCR amplicons of the bacterial 16S rRNA gene V1–V2 region were verified by DNA agarose gel electrophoresis with 1.5% agarose gel and 10× TE buffer. The expected DNA band (∼320 bp) was cut from the gel and DNA was recovered with a QIAEX II Agarose Gel Extraction kit (QIAGEN,

Hilden Germany) and quality verified with a Scandrop spectrophotometer (Thermo Scientific, Vantaa, Finland).

A DNA tagging PCR was used to tag each PCR product of the bacterial 16S rRNA gene V1–V2 region (Chen et al., 2011). According to Chen et al. (2011), the tag primer was designed with four overhanging nucleotides; this arrangement ensured 256 distinct tags, at the 50 end of the 27F and 341R primers for bacterial DNA. The tagging PCR conditions consisted of an initial step of 94◦C for 3 min, 5 cycles of 94◦C for 20 s, 60◦C for 15 s, 72◦C for 20 s, and a final step at 72◦C for 2 min.

Pooled 40 ng lots of each tagged bacterial 16S rRNA gene V1–V2 region DNA samples (88 samples in total) were used for massively parallel pyrosequencing using a Roche 454 Genome Sequencer FLX System (Mission Biotech, Taipei, Taiwan). Raw sequencing reads were sorted according to unique tags using an in-house script, and deposited in NCBI Sequence Read Archive database, under accession number SRP098878, after the primer was removed. Reads from replicate colonies, sampled at the same time and location, were combined into one sample.

## Data Analysis

On a per-sample basis, bacterial reads were quality-filtered using Mothur v1.36.1 (Schloss et al., 2009) to retain reads of lengths from 290 to 340 base pairs (bp), and an average quality score >27. Reads containing any ambiguous base or homopolymer >8 bp were excluded. Chimeric reads were detected and removed using USEARCH v8.1.1861 (Edgar, 2010) with reference mode (3% minimum divergence). Quality-filtered and non-chimeric reads were analyzed using UPARSE pipeline (Edgar, 2013) to generate OTUs per sample (97% identity level), which were classified with RDP classifier (v2.12) (Wang et al., 2007) with a pseudo-bootstrap threshold of 80%.

The rarefied diversity measures (e.g., Shannon, Simpson, Chao 1, ACE, and Good's coverage) were estimated using Mothur with 1000 iterations. The sample size (N) varied from 567 to 29819; therefore, 567 was chosen for rarefaction analysis.

#### Cross-Sample OTU Identification of Acinetobacter and Endozoicomonas

Bacterial OTU representative sequences of all samples affiliated with a given taxonomy were collectively clustered using USEARCH, with options –cluster\_smallmem and –id 0.97. Abundance of each centroid (or cross-sample OTU) was represented by average relative abundance of its cluster members. Phylogenetic analysis was conducted in MEGA7 (Kumar et al., 2016), with maximum-likelihood trees generated with 1000 bootstraps for aligned representative sequences.

#### Resident Types of Bacterial Lineages

The present study defined six resident types of bacteria in S. pistillata and seawater, based on relative abundance and detection frequency. Type I (prevalent abundant group), type II (moderately prevalent abundant group), and type III (non-prevalent abundant group) each had ≥1% average relative abundance and were detected in ≥80, 80–60, and <60% samples, respectively. Type IV (prevalent minor group), type V (moderately prevalent minor group), and type VI (non-prevalent minor group) had <1% average relative abundance and were detected in ≥80, 80–60, and <60% samples, respectively. Coral and seawater samples were separately considered (in terms of detection frequency).

#### Clustering Analysis of Bacterial Taxa

Relative abundance of each classified bacterial genera in individual samples was incorporated into a matrix to estimate the Bray–Curtis distance. Results were presented by non-metric multidimensional scaling (nMDS) analysis using Primer 6 (PRIMER-E, Lutton, Lvybridge, United Kingdom; Clarke and Warwick, 2001) to determine relationships of bacterial communities among times and locations. Analysis of similarity (ANOSIM) via Primer 6 was used to test bacterial diversity indices and compositions in spatiotemporal variations.

For analysis of abundant and minor bacterial members, abundant groups (type I–III) and minor groups (type IV–VI) were separately incorporated into two matrixes. These matrices were used to perform ANOSIM in a Primer 6 based Bray–Curtis distance derived from relative abundance transformed by three methods, including no transformation, square root transformation, and binary transformation (i.e., presence and absence analysis).

Bacterial compositions between sampling times in each location were compared with Pearson correlations using R package Hmisc1 . In addition, Pearson correlations were used to compare fluctuations of Endozoicomonas and Acinetobacter.

Hierarchical clustering analysis was performed using the average linkage via R package pheatmap (Kolde, 2015) and visualized on a heat map. The distance matrix used for clustering was calculated using the R package vegan (Dixon, 2003) based on a base-10 logarithmic transformation of relative abundance (in percentage), plus a pseudo-count of 1 × 10−8 (Costea et al., 2014). Genera of unknown taxonomy or <0.2% relative abundance were excluded.

## RESULTS

#### Diversity of Microbial Community

Seawater sample K1001S yielded the most operational taxonomic units (OTUs), 406.15 (in rarefaction), whereas coral sample K0808C had the fewest (58.68 OTUs; Supplementary Table S2). According to ANOSIM, bacterial diversity indices in corals and seawater samples differed (richness: R = 0.634, p = 0.001; Shannon: R = 0.298, p = 0.001) and bacterial community in seawater samples had greater richness and Shannon diversity than those in corals (Supplementary Table S2).

#### Composition of Bacterial Community

Overall, bacterial compositions in S. pistillata from Kenting, Yehliu, and Lyudao were similar at a class level,

1https://github.com/harrelfe/Hmisc

similar to seawater samples (Supplementary Figure S1). However, bacterial communities in S. pistillata and seawater differed from each other. Dominant bacterial classes in seawater were Flavobacteria, unclassified Marinimicrobia, Gammaproteobacteria, and Alphaproteobacteria. Among these, most dominant genera were unclassified, such as unclassified Marinimicrobia, unclassified Gammaproteobacteria, and unclassified Flavobacteriaceae. In addition, the genus Vibrio, which belongs to Gammaproteobacteria, was higher than those in corals.

In S. pistillata, Actinobacteria, Betaproteobacteria, and Gammaproteobacteria were dominant, and their relative abundances higher than those in seawater. In Gammaproteobacteria, dominant genera were Endozoicomonas, Acinetobacter, Pseudomonas, and Serratia (Supplementary Figure S1), among which Endozoicomonas was particularly more dominant in Kenting and Lyudao compared to Yehliu. Propionibacterium and Actinophytocola were the two dominant Actinobacteria genera, and Massilia Betaproteobacteria in coral.

Bacterial community structure in S. pistillata and seawater at various sampling times is shown (**Figure 1**). Some genera, such as Ralstonia, Cupriavidus, and Aquabacterium, were only detected in coral samples, whereas Pseudoalteromonas was common only in seawater, and Gramella was exclusively retrieved from seawater samples. In addition, among coral samples, although some bacterial taxa were prevalent at every collection, their relative abundance varied from time to time (**Figure 1**). For example, Endozoicimonas and Acinetobacter fluctuated in abundance throughout this survey.

## Variation of Bacterial Composition in Different Time