Diversity and community structure of cyanobacteria and other microbes in recycling irrigation reservoirs

Recycling irrigation reservoirs (RIRs) are emerging aquatic environments of global significance to crop production, water conservation and environmental sustainability. This study characterized the diversity and population structure of cyanobacteria and other detected microbes in water samples from eight RIRs and one adjacent runoff-free stream at three ornamental crop nurseries in eastern (VA1 and VA3) and central (VA2) Virginia after cloning and sequencing the 16S rRNA gene targeting cyanobacteria and chloroplast of eukaryotic phytoplankton. VA1 and VA2 utilize a multi-reservoir recycling irrigation system with runoff channeled to a sedimentation reservoir which then overflows into transition and retention reservoirs where water was pumped for irrigation. VA3 has a single sedimentation reservoir which was also used for irrigation. A total of 208 operational taxonomic units (OTU) were identified from clone libraries of the water samples. Among them, 53 OTUs (358 clones) were cyanobacteria comprising at least 12 genera dominated by Synechococcus species; 59 OTUs (387 clones) were eukaryotic phytoplankton including green algae and diatoms; and 96 were other bacteria (111 clones). Overall, cyanobacteria were dominant in sedimentation reservoirs, while eukaryotic phytoplankton and other bacteria were dominant in transition/retention reservoirs and the stream, respectively. These results are direct evidence demonstrating the negative impact of nutrient-rich horticultural runoff, if not contained, on natural water resources. They also help in understanding the dynamics of water quality in RIRs and have practical implications. Although both single- and multi-reservoir recycling irrigation systems reduce the environmental footprint of horticultural production, the former is expected to have more cyanobacterial blooming, and consequently water quality issues, than the latter. Thus, a multi-reservoir recycling irrigation system should be preferred where feasible.


Introduction
Agriculture uses 70% of global consumptive water [1], and generates huge quantities of runoff which is a major source of impairments to rivers, lakes and wetlands in the United States a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 (http://water.epa.gov/polwaste/nps/agriculture_facts.cfm) and elsewhere [2][3][4]. Capture and reuse of irrigation and precipitation runoff is of global significance to water conservation and environmental sustainability by increasing the efficiency of agricultural water use and reducing nonpoint source pollution. These benefits, however, do not come without cost. Such practices potentially recycle and spread plant pathogens from isolated infections to an entire farm and from a single farm infestation to other farms sharing the same water source [5,6]. Some pathogens may survive better in recycled water than in clean water [7][8][9][10][11]. Water quality may deteriorate, resulting in reduced availability of nutrient elements, underperformance of applied pesticides [12][13][14], and clogging of irrigation pipelines or emitters by microbial mats [15,16].
Photosynthetic agents are considered the driving force affecting several water quality parameters in recycling irrigation reservoirs (RIRs) [15][16][17]. However, the diversity and structure of the phytoplankton community in RIRs is largely unknown. Current research on phytoplankton has focused on non-agricultural freshwater ecosystems. Studies ranged from physical and metabolic capabilities, toxin identification, genome sequencing of cyanobacteria to phytoplankton diversity [18][19][20][21][22][23][24][25][26][27][28]. For example, a recent survey of 46 natural lakes, reservoirs and rivers in Virginia, USA identified more than 300 phytoplankton taxa including cyanobacteria and eukaryotic chlorophytes, diatoms and cryptophytes [28]. Similarly, surveys of five Swedish lakes revealed a linkage between bacterioplankton community composition and environmental variables, providing information on community ecology [20,21]. While inferring the important roles of these microorganisms in global carbon and nitrogen cycles and the quality of public water supplies, these studies have also provided valuable insights into the ecology and potential applications of these microorganisms [18][19][20][21][22][23][24][25][26][27][28].
RIRs may differ in their phytoplankton community composition from natural water bodies for several reasons. First, RIRs receive irrigation runoff containing residues of fertilizers, pesticides and soil microbes as well as plant debris. Second, RIRs may reconfigure runoff through blending with other water sources, progressive dilution in the same reservoir, or channeling to additional reservoirs which may change the composition of the microbial community.
The primary objective of this study was to determine the complexity of cyanobacteria communities in RIRs and its relation to recycling irrigation system setup and other environmental factors. Eight reservoirs and one adjacent stream at three ornamental crop nurseries were surveyed mid to late spring, the busiest growing season of the year for these production facilities.

Selection and configuration of reservoirs and recycling irrigation systems
Eight reservoirs and one adjacent spring-fed stream at three ornamental crop nurseries in eastern (VA1 and VA3) and central (VA2) Virginia were included in this study. The nursery locations are detailed in Table 1 by latitude and longitude. VA1 and VA2 utilize a multi-reservoir recycling irrigation system with a stepwise water flow while VA3 has a single-reservoir system (Fig 1). At VA1 the first reservoir, labeled VA11, receives irrigation and precipitation runoff directly from ornamental crop production areas. Water flows from VA11 to VA12 through a culvert. When VA12 is full, water can flow into VA13 by opening a sealed culvert. Both VA12 and VA13 are pumped for irrigation. VA10 is a small shaded stream that flows along the perimeter of the nursery; it does not receive any runoff water from the nursery production areas and was occasionally pumped to replenish VA11. Three reservoirs at VA2 (VA21, VA22 and VA23) have a similar water flow arrangement but the top reservoir VA21 has a large holding capacity thus water in this reservoir seldom overflowed into VA22 unless there was a severe storm with heavy rain. VA23 receives water from VA22 through a connecting culvert and is the only reservoir used for irrigation. VA31 and VA3X at nursery VA3 are not connected. VA31 directly receives irrigation and precipitation runoff from production areas at this horticultural facility and is pumped for irrigation to adjacent plots. VA3X receives runoff water from nearby agronomic crop fields outside of the property.

Water sampling and environmental data collection
Environmental variables including reservoir function, physical features, surface water quality and nutrient levels are listed in Table 1. Water quality was measured with a 6600V2-4 Multiprobe Datasonde (YSI Inc., Yellow Springs, Ohio, USA) which was deployed in a buoy anchored at the center of each reservoir for continuous data collection. Water samples were taken using 1 liter bottles 20 cm below the surface in a triangular pattern around the buoy at the same time of day in April or May of 2012 or 2013. Water quality spot logs were recorded simultaneously. The collected samples were placed in a cooler and transported back to the laboratory for analyses. Ten ml aliquots of each water sample were immediately filtered through a 0.22-μm membrane and frozen at -20˚C before shipping to Brookside Labs (New Bremen, Ohio, USA) for nutrient analysis. Three replicate measurements/samples were taken from each reservoir or stream.

DNA extraction, PCR and clone library construction
For DNA extraction, 1-liter samples were centrifuged individually at 12,230 g for 50 min in an Avanti J-26XP (Beckman Coulter, Atlanta, GA, USA). The pellets were collected and the supernatant was filtered through 0.22-μm membranes to catch any remaining microorganisms in the sample. DNA was extracted from the pellet and filters with the Power Water 1 DNA Isolation kit (MO BIO Laboratories, Inc, Carlsbad, CA, USA). An equal volume of the three w Reservoirs coded with "1" in the last digit receive runoff directly from production areas of ornamental crops, VA3X received runoff from production areas of agronomic crops. The two reservoirs in VA3 are not connected. Water flows from VA11 to VA13 and from VA21 to VA23.
x All water quality parameters were recorded at sampling. Chl = chlorophyll; DO = dissolved oxygen; EC = electrical conductivity; T = temperature; Turb = turbidity. y NO 2 was not included which was <0.1 mg/L in all reservoirs except VA3X (0.3 mg/L). z Reservoirs from which water was pumped for irrigation. https://doi.org/10.1371/journal.pone.0173903.t001 replicate DNA samples from the same reservoir was pooled and quantified using a DU 800 Spectrophotometer (Beckman Coulter, Atlanta, GA, USA). 20 ng of each pooled DNA sample was used as template in PCRs with cyanobacteria primers CYA359F and CYA781R [29] to amplify the 16S rDNA genes in cyanobacteria and in chloroplasts of eukaryotic phytoplankton. In addition to the template DNA, the 25-μl reaction contained PCR buffer, 0.5 μl of 10 mM each primer, 1 μl of 25 mM MgCl 2 , 1μl of 2 mM dNTPs, 0.25 μl 100×BSA, 0.5 μl Taq polymerase (TaKaRa Bio Inc, Japan), and dH2O. The reaction was run in an Eppendorf Master Cycler with an initial 2 min denaturation at 94˚C, followed by 30 cycles of 30 sec at 94˚C, 30 sec at 55˚C, and 30 sec at 72˚C. PCR products were checked for quantity and size by electrophoresis and ethidium bromide staining. Clone libraries of the PCR amplicons from each reservoir were generated with the TOPO 1 TA Cloning Kit (Life Technology-Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instruction. Over 100 clones per library were randomly selected for sequencing after the 16SrDNA fragment was amplified in the vector with T3 and T7 primers using the same PCR protocol described above. Sequencing was conducted with T3 primer at University of Kentucky's AGTC sequencing facility.

Sequence quality control and editing
Poor quality sequences were discarded according to the quality control report from the sequencing facility. The remaining sequences were aligned, reverse complemented and chopped of cloning vector sequences with Mega6 [30]. Chimera sequences were removed with DECIPHER [31] or Mothur software package V.1.33.3. Other sequence processing steps included alignment, screening, filtering, and preclustering with Mothur. Name and group files were created or edited with BioEdit V7.2.5 (http://www.mbio.ncsu.edu/BioEdit/bioedit.html) [32] and Microsoft Excel.

Construction of reference library for sequence alignment and classification
To construct a reference library for use with Mothur, the taxonomy database (nucleotide) at NCBI (http://www.ncbi.nlm.nih.gov) was searched for "16S rDNA cyanobacteria". The resulting sequences including accession numbers were retrieved with Batch Entrez and then searched on Silva databases (http://www.arb-silva.de/aligner) and classified using SINA [33]. The taxa were then further screened and filtered with Mothur which resulted in 1161 sequences. These sequences were used to construct template files for taxonomic and phylogenetic analysis of the microbial communities in the reservoirs.

OTU construction, classification, phylogenetic analysis and taxonomy
Operational taxonomic units (OTUs) were built with Mothur based on a distance matrix of the sequences from all water samples at a pairwise distance cutoff of 0.15 (a percent similarity threshold of 85%). To analyze phylogeny of OTUs, sequences clustered into the OTUs at cutoffs of 0.03 were used as representative OTUs to make a neighbor joining tree.
OTU sequences were classified through Mothur with the constructed cyanobacteria reference library and Arb-Silva-SINA against Silva and Greengenes 16S gene libraries. When identities of these sequences were larger than the default of 70%, their taxa were accepted. Taxa of "unclassified" were predicted based on classification of nearby known taxa in the phylogenetic tree (S1 Fig) and assigned into one of the three groups identified in this study (Table 2).

Community microbial diversity analysis
Alpha-diversity, an indicator of complexity or richness of the samples, was measured by the rarefaction curve and inverse Simpson index through Mothur. The results were graphed with Microsoft Excel. Beta-diversity (the membership) and microbial community composition as well as their relation to environmental variables were measured with OTU-and taxon-based approaches. For OTU-based analysis, a heat map was generated with the shared file from OTUs of individual reservoirs at cutoffs of 0.03 using Mothur [34] and R version 3.1.0 (http:// www.R-project.org/) [35]. The map included the relative abundance of dominant OTUs that contained five or more clones in each reservoir and a dendrogram showing the relationship between the reservoirs. Variations of taxon composition in the reservoirs of single and multiple reservoir systems were investigated with non-metric multidimensional scaling (NMDS) through the vegan "metaMDS" function in R [36]. To analyze correlation between community composition and environmental variables including reservoir features, water quality and nutrient levels (Table 1), the vegan "envfit" function was used to fit the variables into the NMDS ordination.

Results and discussion
Taxon groups identified from the reservoirs and stream A total of 856 quality sequences were obtained from the clone libraries of the eight RIRs and one stream. VA10 had the least (82) and VA31 had the most (107) clone sequences. When classified through Silva 16S rDNA databases and Mothur with the cyanobacteria reference constructed in this study, 41.2% of these clone sequences belonged to cyanobacteria (Table 2) with the remaining classified as eukaryotic phytoplankton (45.1%) and other bacteria (13.3%). Many algae contain chloroplasts that evolved from endosymbiotic cyanobacteria [37,38] and 16S rDNA primers for cyanobacteria have successfully detected green algae by targeting the chloroplast [39]. All of the identified eukaryotic phytoplankton taxa are listed under the division of cyanobacteria in the current microbial 16S rDNA database. It must be noted that the PCR annealing temperature used in this study was 55˚C which is lower than the recommended 60˚C for cyanobacteria in the original protocol [29]. These results suggest that these cyanobacterium primers may also detect some bacteria under less stringent PCR conditions. Detection of eukaryotic phytoplankton and other bacteria widens the scope of this study to a limited extent. These detections could be a good starting point in understanding microbial community diversity in RIRs. However, if a study's aim is to investigate the entire phytoplankton community, these cyanobacteria primers alone would be inadequate. For such studies, diversitydriven genome sequencing or metagenomic methods with primers or barcodes for other groups of phytoplankton would be better utilized.

Membership of identified groups
Because of the limitation of taxon-based classification, operational taxonomy units (OTUs) were used to classify and count the members in the three microbial groups identified. A total of 208 OTUs were identified from 856 clones. Fifty-seven OTUs contained at least two clones (S1 Table). Twenty-two of these OTUs with 5 to 148 clones each hereafter are referred to as dominant OTUs. The remaining 151 OTUs had a single clone. Among them, 36 belonged to cyanobacteria, 81 belonged to other bacteria and 34 belonged to eukaryotic phytoplankton. Phylogenetic relationships among these OTUs are shown in Fig 2 or S1 Fig  Similar numbers of OTUs were assigned to cyanobacteria and eukaryotic phytoplankton, accounting for 25.5 and 28.4% of total clone sequences, respectively ( Table 2). There were 53 OTUs in the cyanobacterial group with Synechococcus as the most dominant, accounting for 28.3% (79.3% of clones) ( Table 2). The 2 nd to 4 th dominant genera were Pseudanabaena (11.3% of OTUs and 7.3% of clones), Leptolyngbya (11.3% of OTUs and 2.5% clones) and Phormidium (7.5% of OTUs and 1% of clones), respectively. Other genera detected included Arthronema, Chamaesiphon, Cyanothece, Dolichospermum and Microcystis with much fewer OTUs. However, 28.3% of OTUs in the cyanobacteria group remain unclassified with only 81.6% to 93.3% sequence similarities to those in the Silva 16S rRNA database (S1 Table).
Fifty-nine OTUs were identified in the eukaryotic phytoplankton group. Among them, Stramenopiles including diatom, golden algae, brown algae, yellow-green algae and oomycetes were predominant, accounting for 28.8% OTUs and 25.3% clones in the group (Table 2). Chlorophyta or green algae including Trebouxiophyceae were the second most dominant with 20.3% OTUs and 39.3% clones. Cryptophyta, single-celled small algae with two flagella, were the third most dominant taxa with 18.6% OTUs and 15.2% clones. Like cyanobacteria, 26% OTUs of eukaryotic phytoplankton were unclassified.
Other bacteria accounting for 46.2% of total clone sequences included 96 OTUs, almost twice as many as the other two groups. In addition, the percentage of unclassified OTUs was also significantly higher. 99% of these taxa including Candidate (97.3% clones) were not classified (Table 2), however some of them such as MLE1-12 or Melainabacteria have recently been assigned to cyanobacteria [40]. This suggests that other bacteria may be important sources of microbial biodiversity in RIRs, and the high rates of unknown taxa in all three groups may reflect uniqueness of microbial communities in RIRs.

Complexity and structure of microbial communities
Community variations were observed among RIRs. Rarefaction curves based on observed OTUs and detected clone sequences in the communities showed relatively higher diversity in VA10 and VA11 than other reservoirs (Fig 3A). Similar results were obtained when the diversity was measured with inverse Simpson index (Fig 3B). VA10 had the highest index, VA11 had the next highest index, followed by VA23, VA22, VA12, VA31, VA21, VA 3X and VA13. VA10 is a stream that never receives runoff from the production nursery and was occasionally pumped to replenish VA11 during the growing season. Thus, this study provides direct evidence on the negative impact of nutrient-rich agricultural runoff on microbial diversity in these aquatic ecosystems.
There was a general association between microbial community composition and reservoir function (Table 1, Fig 4). First, sedimentation reservoirs (VA11, VA21, VA31 and VA3X) had the highest percentage of cyanobacteria clones ( Fig 4A) and OTUs (Fig 4B) among the three groups. Second, transition/retention reservoirs (VA12, VA13, VA22 and VA23) were  Microbial complexity in recycling irrigation reservoirs predominantly eukaryotic phytoplankton. Third, VA10, the only runoff-free water was predominantly other bacteria. Multi-reservoir water recycling systems and blending with water from streams, rivers or wells may be effective options to avoid irrigation water quality deterioration resulting from cyanobacterial blooms.
There were also site-specificities within the dominant microbial groups. For cyanobacterium-rich sedimentation reservoirs, Synechococcus (OTU1, OTU3 or OTU4) were detected in VA21, VA31 and VA3X but not in VA11 where Pseudanabaena (OTU10) was predominant (Fig 5, S1 Table). Similarly, for eukaryotic phytoplankton-rich transition/irrigation reservoirs, few OTUs were present in both VA1 and VA2. VA13 and V12 of VA1 were predominantly Stramenopiles (OTU2), while inVA22 and 23 of VA2 several OTUs of Cryptophyta were predominant. These site-specificities are likely due to the differences in local soil type, weather conditions and horticultural practices including application of fertilizers, pesticides and plant growth regulators.
Some taxa seemed to have moved along with runoff movement from sedimentation to other reservoirs at nurseries VA1 and VA2. For example, at VA1, Stramenopiles (OTU2, 12, and 17) and Chlorophyta (OTU 6) were detected in both sedimentation reservoir VA11 and transition reservoir VA12; Stramenopiles carried over into the irrigation reservoir VA13 ( Fig  5). Similarly, at VA2, the predominant Synechococcales (OTU 3) in VA21 was also detected in VA22 and VA23. Whether these newcomers in transition/irrigation reservoirs become dominant over time is not known at this time.

Ordination of microbial diversity with environmental variables
Reservoir features, water quality and nutrient levels ( Table 1) were fitted into the community NMDS ordination to determine whether and how environmental variables may affect microbial diversity. Among the 11 investigated environmental variables, nitrate load (NO 3 − ) was the only one substantially influencing the ordination (P < 0.055) ( Table 3), notably in the presence of bacterial Candidate OD1, Parcubacteria, a group favoring anoxic environments (Fig 6). Other variables had relatively lower impacts, and chlorophyll level was least significant. Among these low impact variables, pH, temperature and depth positively associated with the dominance of Synechococcus and Dolichospermum and EC negatively associated with the dominance of Microcystis and Pseudanabaena (Fig 6). These findings are consistent with the previous observations that cyanobacteria are tolerant of high pH and temperature [41] and that Synechococcus grows well under low light [42]. The majority of Synechococcus were found in VA21 and VA3X which had relatively deeper water columns (Table 1). This also explains the lack of Synechococcus in the VA1 communities where reservoirs were shallow at two meters or less.
Detection of eukaryotic phytoplankton and other bacterial species may be associated with some other environmental variables. For example, Chlorophyta, Stramenopile and Trebouxiophyceae dominated VA13, the largest reservoir surveyed in this study. Cryptophyta dominance was positively associated with DO level, as was the dominance of Streptophyta with turbidity, ammonium and orthophosphate. These associations might be indicative of the impact of environmental variables on RIR microbial diversity. Further investigations into whether and how these associations may become significant with time are warranted.

Conclusions
This study characterized the diversity and composition of cyanobacteria and other detected microbes in RIRs at three ornamental crop nurseries and their relation to reservoir function and environmental variables. The diversity of cyanobacteria detected in RIRs is quite different   from that of natural lakes and rivers [28]. This difference highlights the importance of investigation into microbial diversity and population dynamics in RIRs, emerging aquatic ecosystems of global significance. Within the recycling irrigation systems investigated in this study, cyanobacteria was found dominant in the sedimentation reservoirs, whereas eukaryotic phytoplankton were dominant in the transition and retention reservoirs, and other bacteria were dominant in the runoff-free stream. These results are direct evidence that nutrient-rich horticultural runoff, if not captured and retained on the production properties, will negatively influence microbial diversity in natural water resources including lakes, streams and rivers. Practically, both single-and multi-reservoir recycling irrigation systems reduce the environmental footprint of horticultural production; the former is expected to have more cyanobacterial blooming, consequently water quality issues, in RIRs than the latter. Thus, a multi-reservoir recycling irrigation system is preferred where feasible. Further investigations into how microbial community diversity and structure in RIRs may change seasonally are warranted.