Prorocentrum donghaiense was kindly provided by the Collection Center of Marine Algae, Xiamen University, China. All the cultures were aseptically handled and grown in autoclaved and 0.22 μm filtrated seawater with f/2 -silicate medium (Guillard, 1975). Seawater used in this study was collected from the South China Sea and its background DIP concentration was measured (see below for the measurement method) and found to be low (0.73 μM). Cultures were maintained at 20°C under a light: dark cycle of 14:10 h with a photon flux of 100 μmol·m⁻²·s⁻¹.

Stock solution of pure glyphosate (99.9%, non-derived compound; SIGMA -ALDRICH) was prepared at a final concentration of 36 mM, sterilized by filtration through 0.22-μm membranes, and stored at 4°C. P. donghaiense cultures were set up in five different P nutrient conditions (f/2, f/2^−P, f/2^−P+Gly, f/2^LG, and f/2^HG). Due to the lack of an accurate method for glyphosate measurement in coastal waters, it is difficult to know the glyphosate distribution in the ocean. Thus, in this paper, 36 μM glyphosate (group f/2^−P+Gly) was used, a concentration equivalent to the phosphate concentration in the f/2 medium, as the sole P-source for comparison with the DIP-depleted cultures, f/2^−P (DIP concentration was <1 μM). Furthermore, 36 μM (f/2^LG) and 360 μM (f/2^HG) glyphosate were added to the normal (f/2) media to test for toxic effects. A high concentration of glyphosate was used because: (1) we wanted to maximize the chance to elicit a response of the bacterial community in the P. donghaiense culture and (2) these concentrations were used in our previous study in which 360 μM glyphosate was found to cause severe inhibition for most phytoplankton species (Wang et al., 2016), except P. donghaiense. Each treatment was set up in triplicate and all experiments were run in a volume of 400 mL in 600 mL Cell and Tissue Culture Flasks (Jet Biofil). Because dinoflagellates do not grow well under mechanical disturbance, the cultures were not agitated but were mixed gently whenever a sample was taken. At the end of the experiment (day 10), a slight increase in cell concentration was observed in the last two media (LG and HG). To further investigate the effects of glyphosate on these cultures, f/2^LG and f/2^HG were extended to 122 days.

Experiments with glyphosate began when the stock cultures entered the exponential growth stage. For all five P-nutrient conditions described above, the initial cell concentration of P. donghaiense was set to ~10% of the highest concentration, or 8,000 cells/mL, to provide room for growth. To monitor growth, 1 mL culture was removed to perform a cell count every 2–3 days using a Sedgwick-Rafter counting chamber under the microscope (LeGresley and McDermott, 2010). Meanwhile, DIP was measured using the molybdate-ascorbic method by filtrating 25 mL of culture (Parsons, 2013) and the detection limit of the machine was 0.01 μM. These measurements were carried out for a period of 10 days, until growth of all treatments reached a plateau except for the glyphosate-treated cultures, which showed no significant growth throughout. After day 10, the cultures were checked approximately every 2 weeks. On day 122, remarkable growth was observed in the cultures with glyphosate and concentrations of DIP or DIN (dissolved inorganic nitrogen) were again determined.

On day 122, about 25 mL of culture was collected from each treatment and was filtrated by using a sterilized 0.22-μm nitrocellulose membrane. The membrane was then transferred into a 1.5 mL microcentrifuge tube with 0.8 ml of DNA extraction buffer (containing 0.1 M EDTA, 1% sodium dodecyl sulfate, 8 μg lysozyme, and 18 μg proteinase K). The samples were incubated at 55°C for 3 days, with the prolonged incubation intended to maximize cell breakage. DNA extraction was then carried out using a CTAB protocol combined with the Zymo DNA Clean & Concentrator kit (Zymo Research Corp., Orange, CA), as previously reported (Zhang et al., 2005).

From each sample, 1 μL (~80 ng) of the extracted DNA was used as the PCR template to amplify the full-length 16S rRNA gene. Amplification was run in the following program; initial denaturation at 95°C for 3 min, followed by 35 cycles of denaturation at 94°C (30 s), annealing at 56°C (30 s) and extension at 72°C (1 min), with a final extension cycle at 72°C for 10 min. Degenerate primers 27F [5′-AGAGTTTGATCMTGGCTCAG-3′ (M is A or C)] and sgR (5′-TAGGGTTACCTTGTTACGACTT-3′) were used as they were previously found to amplify approximately a 1.5-Kb fragment of 16S rRNA in a wide range of bacteria (Webster et al., 2003; Xu et al., 2011; Zhou et al., 2014). The amplicons were gel purified and cloned. The resulting clones were randomly picked from each treatment for Sanger sequencing.

Although Sanger sequencing provides long reads (~800 bp) and hence higher taxonomic resolving power, its data output from a feasible sequencing scale is limited and as such, limits the depth coverage of existing taxonomic diversity. To address this issue, we also conducted Illumina high-throughput sequencing. V4-V5 specific primers, 515F (forward primer, 5′-GTGCCAGCMGCCGCGG-3′) and 907R (reverse primer, 5′-CCGTCAATTCMTTTRAGTTT-3′), were used to amplify a 390-bp 16S rRNA gene fragment from the DNA template. PCR was carried out as the following: 95°C for 2 min, followed by 30 cycles of incubation at 95°C (30 s), 55°C (30 s) and 72°C (30 s), with a final extension cycle of 5 min at 72°C. The amplicons obtained through three independent PCRs for each of the DNA template preparations were purified and subjected to the Illumina Miseq (Majorbio Co. Ltd., Shanghai, China).

Nucleotide sequences acquired through Sanger sequencing were subjected to a chimera check using the Bellerophon v.3 program (Huber et al., 2004; http://greengenes.lbl.gov) with default settings. Sequences flagged as potential chimeras were examined further by comparing the BLAST results of both the sequences flanking the breakpoint with suspected parental sequences. If both the parental sequences were from the same PCR reaction of our own sample, we ranked the sequence as a potential chimera. Next, the end sequences of the same potential chimeric sequence were BLASTed against GenBank database separately. If the two fragments were highly identical (>97%) to different species of bacteria, the sequence was ranked as a chimera. After the removal of the chimera sequences, the remaining sequences were aligned using MUSCLE on MEGA (Hall, 2013) and the criteria to define an operational taxonomic units (OTUs) was set as 97% similarity.

All the sequences obtained from MiSeq are available in the Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra). We used the Trimmomatic (Bolger et al., 2014) and FLASH (Magoč and Salzberg, 2011) software packages to conduct quality filtering of raw reads. Reads with trailing bases of below quality 20 (Phred Quality Score) were trimmed and a sliding window of 50 bp was used to trim tail reads that showed an average quality within the window < 20. Reads shorter than 50 bp were discarded. Setting the minimum overlap length as 10 bp and maximum mismatch ratio as 0.2, pair-end reads were assembled into contigs. These contigs were then assigned to treatment groups according to barcode with maximum allowed mismatch of 2 bases.

The contigs were clustered at 97% identity cutoff (generally considered to define microbial species) into OTUs to calculate rarefaction and non-parametric estimators (such Chao1 and ACE) using UPARSE (version 7.1) (Edgar, 2013). Chimeras were detected and removed using the UCHIME pipeline (Edgar et al., 2011). Taxonomic assignments of the OTUs were performed using the RDP classifier (version 2.2) (Wang et al., 2007). The SILVA database (release 115) (Quast et al., 2012) was used for the retraining of RDP classifier and the minimum confidence threshold was set to 0.70. Only five sequences in the V4-V5 library remained unclassified below the domain Bacteria (Table 1A).

For the prediction of functional profiles of the bacterial community based on our 16S rRNA gene sequences, data from each treatment and the control were run through the newly developed open-source R package Tax4Fun (Aßhauer et al., 2015) after being processed with Qiime (based on SILVA).

Copy numbers of two phylotypes (16S rRNA gene) and phnJ were determined using qPCR (Quantitative PCR), which was conducted on the CFX96 Real-Time System (Bio-Rad, USA). The two 16S rRNA gene phylotypes belonged to Actinobacteria and Gammaproteobacteria, which were not detected in the control but had high copy numbers in the LG or HG groups. Specific primers for the 16S rRNA genes of these two phylotypes were designed based on the sequences we obtained. For phnJ, we first used degenerate primers (PhnJF1 and PhnJR1; Table 2) that were based on the alignment of thirteen phnJ sequences from GenBank's Gammaproteobacteria to identify conserved regions of this gene (Figure S2). Using these primers in PCR, we obtained a gene fragment of about 340 bp and sequenced 10 clones which revealed that all the eight sequences (Figure S3) we obtained from the colonies hit the phnJ gene in the database (the other 2 were not a resolved single sequence and were discarded). Based on these sequences, we designed several sets of primers for qPCR. The specificity of these primers was tested using qPCR melt curve analysis. One of the primer pairs was found to amplify one single product (melt curve showing single peak; Figure S4) and was selected for use in qPCR (PhnJF2 and PhnJR2; Table 2). The PCR products were purified and prepared in successive 10-fold dilutions to be used as standards. The same molar quantity of genomic DNA (5 ng) was used as the template to amplify the target genes from each sample group (C, LG, and HG). Their copy numbers were calculated based on the standards, which were amplified on the same PCR runs as the samples. Each PCR reaction was carried out in a total volume of 12 μL containing 6 μL of 2 × iQSYBR Green supermix (Bio-Rad, USA), 375 nM of each primer, and 6 μL of 5 ng/μL DNA. The PCR program was composed of a denaturation step of 3 min at 95°C, followed by 40 cycles of 95°C for 10 s and 60°C for 32 s. In each run, negative controls were set up without DNA templates. Each reaction had three technical replicates. At the end, to confirm primer specificity, all the PCR products were subjected to melting curve analysis.

For all variables, differences among treatments were examined using repeated measure analysis of variance (RM ANOVA; Winer et al., 1971), with three glyphosate treatments (f/2, f/2^−P, f/2^−P+Gly and f/2, f/2^LG, f/2^HG, separately) and sampling times as treatment matrices. RM ANOVA analyses were followed by all pairwise multiple comparisons (post hoc testing), using the Holm–Sidak method (P < 0.05).

Starting with equal cell concentrations, the f/2-group grew exponentially and peaked on day 8, while the f/2^−P and f/2^−P+Gly treatments entered the plateau phase on day 4 (Figure 1A). Statistical analysis showed that there was a significant difference in growth between groups f/2 and f/2^−P (p < 0.05, RM ANOVA) from day 1 to day 10, while no significant difference was observed between f/2^−P+Gly and f/2^−P (p > 0.05, RM ANOVA). Meanwhile, DIP concentrations in the f/2 cultures decreased steadily whilst those in the f/2^−P and f/2^−P+Gly treatments were low to begin with (limited to the carry-over from the culture inoculum or seawater) and further dropped to below the detection limit from day 4 (Figure 1B). These results indicated that P. donghaiense could not hydrolyze glyphosate to obtain phosphate. On day 64, the f/2 and f/2^−P almost died due to depletion of available P-source whereas the f/2^−P+Gly group could still maintain slight growth (Figure 1A), showing a positive albeit delayed effect caused by glyphosate. Furthermore, on day 122, significant increasing growth was observed in f/2^−P+Gly compared to f/2 and f/2^−P. However, DIP in f/2^−P and f/2^−P+Gly cultures was still below the detection limit on this day. The late-stage increased growth in f/2^−P+Gly might be because glyphosate induced bacteria in the culture to degrade glyphosate and release phosphate, which was taken up immediately by P. donghaiense. We also noticed that as early as day 10, the cell concentrations of both the f/2^LG and f/2^HG groups showed higher (1.4- and 1.7-folds, respectively) than that in the f/2 group (Figure 1C). On day 64 and 122, when the f/2 cultures declined dramatically due to nutrient depletion, groups f/2^LG and f/2^HG exhibited a bloom equivalent condition, with a higher biomass in f/2^HG, and this trend lasted for at least 2 months (Figure 1C). Correspondingly, the color of cultures became increasingly darker with increasing concentrations of glyphosate among these three groups (pictures not shown). And on day 122, DIP concentrations in all three groups were < 4 μM, indicating that DIP was nearly depleted. DIN in f/2^LG and f/2^HG was also depleted (4.5 and 4.7 μM, respectively) compared to the abundant presence in f/2 (203 μM; Figure 1D), indicating further nitrogen consumption and growth in the glyphosate-based cultures, beyond what f/2 cultures could achieve.

Bacterial communities in the C, LG, and HG cultures were compared on day 122 based on 16S rRNA gene. We randomly picked and sequenced 150 clones derived from 16S rRNA gene amplicons. After excluding 27 sequences because of evidence of sequencing error or likely being chimeric DNA, the remaining good-quality sequences were classified into different phylotypes to characterize changes in major functional groups among C, LG, and HG (Figure 2). It is clear that the microbial community was composed of Proteobacteria (in declining richness order from Alphaproteobacteria, Gammaproteobacteria to Deltaproteobacteria), Bacteroidetes (Class Sphingobacteria), and Planctomycetes.

Results showed that there were 17 phylotypes identified in the control, while 9 and 22 phylotypes in LG and HG respectively. A similar trend was found for Shannon–Wiener diversity index, a commonly used diversity index in ecological literature (Spellerberg, 2008), which was 2.28, 1.80, and 2.71 for C, LG, and HG, respectively. This result indicates that 36 μM glyphosate caused a decrease in the bacterial diversity but 360 μM glyphosate caused an increase. Of the bacteria detected, Alphaproteobacteria were dominant and their richness and abundance went up along with the increase in glyphosate concentration. As for Gammaproteobacteria, a higher richness was observed in C in comparison to LG and HG, displaying a depressed effect of glyphosate on this group. By contrast, we identified Deltaproteobacteria only in C and HG. No significant difference was observed in richness and abundance among group C, LG, and HG for Sphingobacteria, showing that this might be a core group of bacteria constitutively associated with P. donghaiense. Only one phylotype was identified in Planctomycetes and it only occurred in HG.

More than 53,000 16S rRNA gene V4-V5 sequences (Table 1A) were obtained using Illumina high-throughput sequencing. Most of the curated sequences (exceeding 99.9%) were classified to a particular bacterial phylum and class (or lower taxon). Sequences that were unclassified below the level of Bacteria probably represented a mixture of true biodiversity of novel taxa and complex chimeras that were not identified by UCHIME pipeline (see section Materials and Methods). In this dataset, different read numbers were obtained in groups C, LG, and HG (Table 1B), but the rarefaction curves were nearly asymptotic (Schopf, 1975) for C and LG, while our sampling for HG was less but still nearly complete (Figure S1).

Results from Sanger sequencing of the cloned libraries indicated a significant shift in the bacteria community present in the P. donghaiense cultures. In the glyphosate-treated cultures, Alphaproteobacteria and Planctomycetes were promoted markedly, because the abundance of the former increased significantly in LG and HG compared with the control and Planctomycetes were only detected in HG. In contrast, Gammaproteobacteria showed reduced abundance with the presence of glyphosate. These results indicate that Alphaproteobacteria and Planctomycetes responded positively while Gammaproteobacteria responded negatively to glyphosate. The latter is at odds with previous studies that showed Gammaproteobacteria were able to utilize glyphosate as sole P-source (Moore et al., 1983; Quinn et al., 1988). We further sequenced the V4-V5 region in 16S rRNA gene using Illumina high-throughput sequencing for a more complete picture of the microbial community in the P. donghaiense cultures.

The high-throughput sequencing result showed that the most abundant phylotypes shared among these three groups belonged to the family Rhodobacteraceae (Figure 3), which are aquatic bacteria that frequently thrive in the marine environments (Pujalte et al., 2014). It is reported that Rhodobacteraceae is one of the three most dominant algal bloom-associated bacterial classes and several laboratory studies have identified specific associations between phytoplankton and certain species of Rhodobacteraceae, making it a major model for the study of microorganism-phytoplankton interactions (review in Buchan et al., 2014). The presence of this group of bacteria in all treatment groups suggests that Rhodobacteraceae is also the major bacteria co-occurring with P. donghaiense.

Among the five clusters we detected, clusters C and E responded to glyphosate positively (Figure 4) indicating that they are better adapted to glyphosate and may be able to metabolize glyphosate. Most of the bacteria in these two clusters came from Alphaproteobacteria and Gammaproteobacteria. Other phylotypes identified in these two clusters included Firmicutes, Flavobacteria, Actinobacteria, and Planctomycetes. Of these, Flavobacteria contains species that are usually most abundant in coastal waters and in phytoplankton blooms (reviews in Buchan et al., 2014), suggesting these marine bacteria, if proven to be able to utilize glyphosate, are likely to play an important role in glyphosate metabolism.

By contrast, cluster A showed hormetic response to glyphosate, a dose response phenomenon characterized by a low dose stimulation and high dose inhibition (Axelrod et al., 2004). For instance, Methylotenera spp. in cluster A was more abundant in LG than in the control but its abundance decreased by 70-fold in HG, suggesting that Methylotenera spp. is resistant to glyphosate at micromolar concentrations but sensitive at millimolar levels. This is common in bacteria: the growth of many cyanobacteria strains were only inhibited by glyphosate at concentrations exceeding 0.3 mM (Forlani et al., 2008). Pseudomonas and many other bacteria isolated from sprayed sites all showed high glyphosate resistance even under the concentration of 10 mM (Quinn et al., 1988; Staub et al., 2012). Furthermore, cluster B showed another dose-dependent response to glyphosate. Almost all of these phylotypes were from the aquatic environment and are usually abundant in coastal regions and phytoplankton blooms (Bratbak et al., 1996; Fandino et al., 2001; Pinhassi et al., 2004; Alonso et al., 2007; Buchan et al., 2014). They all showed dose-dependent increases in mortality caused by glyphosate and this dose-dependent induction effect is commonly found in some cyanobacteria and eukaryotic algae (Forlani et al., 2008; Lipok et al., 2010) as well as fish and other amphibious (Annett et al., 2014). Additionally, glyphosate could cause reduction (or disappearance) of some lineages of bacteria such as Pseudomonas spp., Legionella spp., Aquarelle spp., Xanthobacter spp., and Marinobacter spp., suggesting that they are quite sensitive to glyphosate. Majority of them belong to Gammaproteobacteria, which has been described above as being able to degrade glyphosate. This indicates that glyphosate could change the bacterial composition significantly. Moreover, shifts in microbial community composition were also observed in soil samples exposed to high concentrations of glyphosate (>10 mM; Araújo et al., 2003; Ratcliff et al., 2006; Lancaster et al., 2010). Our study shows that glyphosate can also exact influence on the structure of marine bacteria associated with phytoplankton. All the results in concert support our first hypothesis that glyphosate can induce a shift in bacterial community associated with P. donghaiense.

If an alga could utilize glyphosate, an immediate growth effect would be observed on the alga, as previously reported (Wang et al., 2016). Therefore, our observation that P. donghaiense did not grow in the first 60 days clearly indicates that this dinoflagellate is unable to utilize glyphosate directly. This result is consistent with the previous finding that most eukaryotic phytoplankton examined especially dinoflagellates could not utilize glyphosate or other phosphonates (Cui et al., 2015; Wang et al., 2016). So far many cyanobacteria and heterotrophic bacteria have been shown to be able to hydrolyze glyphosate to scavenge phosphate or use the organic carbon (Hove-Jensen et al., 2014), while only few eukaryotic species have been shown to be able to grow with glyphosate provided as sole P-source (Wang et al., 2016). Our observation that P. donghaiense showed remarkable growth after 120 days suggests that after the DIP was used up, some specific lineages of bacteria were induced by glyphosate and grew and degraded glyphosate to release sufficient phosphate to support the growth of P. donghaiense. This suggests that the bacteria might have originally been in the culture of P. donghaiense in very low abundance. Although it is also possible that these bacteria might have been introduced from outside the culture subsequently, this is not so likely because the cultures have not been disturbed since day 10.

Both our physiological and molecular data support our second hypothesis that the bacterial community shift induced by glyphosate involves the increase in glyphosate-degrading species. First, the most remarkably induced groups of bacteria (Clusters A and C) in the presence of glyphosate were from Alphaproteobacteria and Gammaproteobacteria, which have already been reported to utilize glyphosate as a sole P-source (Sviridov et al., 2014). Other phylotypes identified in these two clusters included Firmicutes, Flavobacteria, Actinobacteria, and Planctomycetes. Several species in Firmicutes have been shown to contain the phn cluster, a group of genes responsible for phosphonate degradation (Huang et al., 2005). Furthermore, the phosphonatase pathway genes responsible for phosphonate degradation were also detected in many species in Alphaproteobacteria, Gammaproteobacteria, Actinobacteria, and Planctomycetes. Planctomycetes have also been shown to utilize 2-AEP and methylphosphonate as a sole P-source to support growth (Martinez et al., 2010).

Second, a search for phn genes (the C-P lyase pathway) against available DNA and protein sequence databases also found significant hits to phnJ (the core gene in C-P lyase pathway) in Alphaproteobacteria (Villarreal-Chiu et al., 2012). Furthermore, our qPCR results indicated increase in the phnJ gene copy number in the glyphosate treated groups. The phnJ copy number from the samples might have been underestimated, however, because the degenerate primers we used were designed from only sequences of Gammaproteobacteria (Table S2), leading to only a dominant phnJ sequence for qPCR primer design. However, the Gammaproteobacteria were a dominant group of the bacteria in our cultures and the increasing trend of phnJ copy number in the glyphosate-treated cultures was likely representative of the general trend of phnJ gene abundance under glyphosate treatment. In addition, our Tax4Fun analysis predicted that the abundance of genes in the bacterial assemblage involved in phosphonates uptake and metabolism consistently increased as glyphosate concentration went up, with only one exception (Table S1). This too might have underestimated the ability of the bacterial assemblage to degrade glyphosate because Tax4Fun can only cover taxa that are represented in the genome database and yet not all of our retrieved 16S rRNA sequence data could be mapped to the existing genome database. This is evident as phnJ, which was detected in our cloning and qPCR analyses, was not in the predicted genes (Table S1). Nevertheless, our observation of significantly higher abundances of the predicted glyphosate-metabolizing genes (except one) in glyphosate-amended cultures when compared to the control is indicative of an increase in glyphosate-degrading bacteria; thus, supporting our second hypothesis. In sum, presence of glyphosate in the environment can impact the acquisition of P by the co-occurring phytoplankton species on the one hand and affect the microbial processes in the marine ecosystem on the other.

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.