Abstract
Standard aquaculture generates large-scale pollution and strains water resources. In aquaculture using zero discharge systems (ZDS), highly efficient fish growth and water recycling are combined. The wastewater stream is directed through compartments in which beneficial microbial activities induced by creating suitable environmental conditions remove biological and chemical pollutants, alleviating both problems. Bacterial predators, preying on bacterial populations in the ZDS, may affect their diversity, composition and functional redundancy, yet in-depth understanding of this phenomenon is lacking. The dynamics of populations belonging to the obligate predators Bdellovibrio and like organisms (BALOs) were analyzed in freshwater and saline ZDS over a 7-month period using QPCR targeting the Bdellovibrionaceae, and the Bacteriovorax and Bacteriolyticum genera in the Bacteriovoracaeae. Both families co-existed in ZDS compartments, constituting 0.13–1.4% of total Bacteria. Relative predator abundance varied according to the environmental conditions prevailing in different compartments, most notably salinity. Strikingly, the Bdellovibrionaceae, hitherto only retrieved from freshwater and soil, also populated the saline system. In addition to the detected BALOs, other potential predators were highly abundant, especially from the Myxococcales. Among the general bacterial population, Flavobacteria, Bacteroidetes, Fusobacteriaceae and unclassified Bacteria dominated a well mixed but seasonally fluctuating diverse community of up to 238 operational taxonomic units, as revealed by 16S rRNA gene sequencing.
Introduction
Aquaculture has become a major source of food, generating more than 45% of global fish production (FAO, ). This success comes with a heavy environmental price, as aquaculture requires huge quantities of water for operational uses, discharges heavy loads of organic and inorganic contaminants to water bodies, and consumes natural coastlines and their ecosystems (Subasinghe et al., ). Inland water recycling systems help reduce this environmental burden, shifting the operations away from sensitive areas, and diminishing the demand for water. Yet, pollutants are still discharged into the environment, as organic and inorganic nitrogen compounds are mostly untreated (van Rijn, ). In zero discharge systems (ZDS) where fish basins (FBs) and treatment compartments are connected to form continuous water recycling units, nitrogen, sulfate and organic matter are removed by microbial processes (Neori et al., ). ZDS are thus based on structured but connected environments that form mesocosms reflecting rather varied ecological conditions. This structure brings about continuous and efficient water purification leading to adequate water quality conditions for aquacultural needs over time, with water requirements limited to compensation for evaporation losses. Significantly, high fish yields are obtained under relatively disease-free conditions (Shnel et al., ; Gelfand et al., ).
The major processes responsible for water purification occurring in ZDS are nitrification, denitrification, sulfur oxido/reduction and degradation of organic matter. They are performed by numerous groups of bacteria, some of which have been identified and partially characterized (Cytryn et al., , , b). Other ecological processes, such as microbial predatory interactions may affect the functioning of the microbial communities sustaining water purification, for example predation by viruses and protists in microbial food webs cause large-scale bacterial mortality (Fuhrman & Noble, ; Pernthaler, ). In aquatic environments, microbial biomass can be consumed at approximately the same rate as it is produced (Fuhrman & Noble, ; Pernthaler, ) and as such, microbial predation has a significant ecological effect. It has also been shown that microbial predation affects the evolution and the adaptation of both predators and prey, and modifies the structure and function of microbial communities (Gallet et al., , ).
The understanding of the ecological relevance of microbial predatory interactions reposes with research with protists and with phages. Much less is known on the effects of bacterial predators of bacteria in natural or man-made systems, although these organisms are commonly found in the environment and are widely distributed (Jurkevitch & Davidov, ). Recent studies have shown that Bdellovibrio and like organisms (BALO), a major group of predatory bacteria, rapidly respond to changes in the structure of bacterial communities (Chauhan et al., ; Chen et al., ).
BALOs mainly belong to the Deltaproteobacteria, where they form two distinct and very diverse families. The Bdellovibrionaceae, with two defined species, Bdellovibrio bacteriovorus and Bdellovibrio exovorus, are in freshwater (including wastewater) and soils. The Bacteriovoracaceae comprise the halophilic species Bacteriovorax marinus sp. and Bacteriovorax litoralis sp. (Baer et al., ; Davidov & Jurkevitch, ), the freshwater/soil Bacteriolyticum stolpii (Pineiro et al., ), and the Peredibacter genus (Davidov & Jurkevitch, ). An additional BALO, the obligate predator Micavibrio sp., forms a distinct group in the Alphaproteobacteria (Davidov et al., ). As to their effect on microbial communities, BALOs differ in prey range, from generalists to being rather restricted in range (Jurkevitch et al., ; Shemesh et al., ; Chen et al., ). BALOs may exert significant yet unexplored effects on the marine ‘microbial loop’ (Pomeroy et al., ), altering the structure of bacterial communities and affecting population dynamics. More specifically, BALOs have the potential to alter particular ecological functions by affecting the abundance of prey populations expressing them. As an example, BALOs isolated from fish ponds were shown to reduce disease incidence caused by the fish pathogens Aeromonas hydrophila and Vibrio alginolyticus (Chu & Zhu, ; Wen et al., ) and contributed to the improved growth performance of fish, shrimp, crab and sea cucumber (Qi et al., ). Thus, BALOs may affect the functioning of water recycling systems, and in aquaculture may help keep fish stocks healthy.
To assess the impact of BALOs on microbial processes it is imperative to identify the predators and to quantify their absolute and relative abundance. In this study, we explored the dynamics of BALO populations, including the Bdellovibrionaceae and the Bacteriovorax and Bacteriolyticum genera of the Bacteriovoracaceae, in two ZDS with different salinities over a period of 7 months. In the freshwater system, we examined the community composition of the total bacterial community in detail and provided evidence for the hypothesis that environmental parameters may affect BALO predatory dynamics and population structure. Although the predators do not dominate numerically, they form fairly abundant populations that fluctuate over time and differ according to salinity levels. Remarkably, freshwater Bdellovibrionaceae were actually found to populate the saline ZDS also.
Materials and methods
The zero discharge systems
A pilot plant comprising two parallel ZDS operates at the Rehovot campus (Supporting Information, Fig. S1). One ZDS operates with freshwater, the other with water maintained at a salinity level of 20 p.p.t. Each system consists of a 5-m3 fish culture basin (FB) from which water is withdrawn through two parallel biofiltration loops: (1) Water from the upper layers of the FB is treated in a 7-m3 trickling filter (TF), filled with PVC substratum, for oxidation of ammonia to nitrate and for degassing of CO2. The retention time of the water flowing through this loop is 30 min. (2) Water rich in particulate organic carbon from the bottom of the FB is diverted into a 2-m3 digestion/sedimentation basin (DB) where oxidation of organic material is mediated by a variety of electron acceptors including oxygen, nitrate and sulfate. Effluent from the DB is then pumped into a 13-L fluidized bed reactor (FBR) for additional nitrate and sulfide removal. The flow rate through this loop is 0.8 m3 h−1. Following treatment, purified water from both loops is returned to the FB. Fish were stocked on 31 October 2010. The fresh water basin was stocked with 650 tilapia hybrids (Oreochromis niloticus× Oreochromis aureus) with an average weight of around 50 g. The salt water FB was stocked on 31 October 2010 with 738 gilt-head sea bream (Sparus aurata) with an average weight of 1.5 g. The fish were fed a high protein and lipid diet (45% and 19%, respectively). The diet is assimilated to up to 60%, but about 80% of the nitrogen is excreted in the form of ammonium. The tilapia and gilt-head sea bream reached an average of about 250 and 240 g at harvest, respectively. Harvesting was performed in July 2011.
Sampling
Water samples were collected in triplicate from each of the compartments of the fresh and salt water systems on a monthly basis, totalling 21 samples per time point, from December 2010 to June 2011, except for February 2011. Samples of 500 mL from the fish and DBs were sterilely collected from the upper 25 cm of the water phase using a self-constructed sampler. For the TF, water trickling down at the bottom of the filter was collected, whereas for the FBR water and some sand were collected from the base of the reactor. Collected samples were immediately stored at 4 °C before processing. Biofilm samples were collected by inserting a long plastic strip in the TF and allowing a biofilm to develop on it. To relate biofilm measurements with planktonic phase measurements, the biofilm was assumed to be 20 μm thick. Water quality parameters in the FBs were measured daily (temperature and dissolved oxygen) or twice monthly (pH, nitrite, nitrate, total ammonia, orthophosphate and sulfide levels).
DNA extraction
A 50-mL sample of those collected was filtered through a 0.2-μm filter (Schleicher and Schuell) to capture the bacteria. Half of the filter was cut into pieces and used for DNA extraction using a Powersoil™ DNA Isolation Kit (Mo Bio Laboratories, Inc.) according to the manufacturer's protocol. DNA was isolated from biofilms by cutting 10-cm2 strips into sterile water followed by vigorous vortexing until the attached material was cleaned off the surface, and then proceeding as above. Finally, the DNA bound to the silica gel of the extraction column was eluted with 50 μL PCR grade double-distilled water (DDW). DNA concentration was measured with a Nanodrop 2000 spectrophotometer (Olympus) and stored at 4 °C until use.
BALO-specific PCR
PCR was performed with the DNA samples extracted from the ZDS using 16S rRNA gene primers designed by Davidov et al. (). These primers specifically target the Bdellovibrionaceae (Bd) family and the following clades within the Bacteriovoracaceae (Bc) family: Bacteriovorax (the halophilic BALOs) and Bacteriolyticum (freshwater and soil BALOs) (Table S2) but not Peredibacter, as confirmed by PCR testing using two strains from this genus (data not shown). The PCR reaction included 12.5 μL of PCR mastermix (Cat no. D123P; Lambda Biotech), 1 μL of each primer (10 μm), 2 μL of DNA and 8.5 μL of PCR grade DDW. Positive (BALO DNA) and negative (no DNA) controls were always included. The PCR was performed with a Mastercycler Gradient (Eppendorf AG, Germany). Thermal profiles for Bdellovibrionaceae were: denaturation and enzyme activation at 95 °C, 2 min, followed by 36 cycles of 95 °C for 30 s, annealing at 51 °C for 30 s, extension at 72 °C for 30 s, and final extension at 72 °C for 5 min. For the Bacteriovoracaceae the thermal profiles were similar except that annealing was at 55 °C for 30 s followed by an extension step at 72 °C for 45 s. PCR products (6 μL) were verified by gel electrophoresis [1% agarose gel in 0.5% Tris acetate-EDTA (TAE) buffer stained with 1 μL of ethidium bromide; 100 V for 30 min] and the gel was digitized with a B.I.S Bioimaging System (Pharmacia Biotech, Sweden).
Standard curves for BALO-specific QPCR
The predatory strains used as reference strains in the PCR and QPCR assays (Table S1) were grown as previously described (Jurkevitch, ). Bd and Bc abundances were calculated using standard curves generated by serially diluting plasmid preparations containing a fragment of the 16S rRNA gene from B. bacteriovorus 109J or from B. stolpii UKi2, respectively. The plasmid containing the B. bacteriovorus gene fragment was constructed by inserting a 492-bp fragment of the gene, amplified with primers BbsF216 and BbsR707 (Table S2), and using the following thermal cycling profile (modified from Van Essche et al., ): 95 °C, 2 min, 29 cycles of 95 °C for 30 s, 56 °C for 30 s, 72 °C for 30 s, with a final extension at 72 °C for 5 min. The plasmid containing the B. stolpii gene fragment was constructed using a 981-bp fragment from the 16S rRNA gene amplified with the primers Bac69F and Bac1049R (Table S2) designed using the NCBI primer-blast program (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). The PCR reaction parameters used to amplify this fragment were as above, except that annealing was performed at 56 °C and extension was run for 60 s.
BALO-specific QPCR
Bdellovibrionaceae
QPCR for the Bd was performed based on the TaqMan technology, as in Van Essche et al. (), using primers Bd347F and Bd549R and probe Bd396P, labeled with FAM (6-carboxyfluorescein) at the 5′ end, and TAMRA (6-carboxytetramethylrhodamine) at the 3′ end (Sigma Aldrich; Table S2). Reaction tubes contained 12.5 μL of Jumpstart Taq Readymix for High Throughput Quantitative PCR (D6442; Sigma-Aldrich), 2 μL of each primer (10 μM), 1.25 μL of the Taqman probe (1 μM), 5 μL of DNA and 2.25 μL PCR grade DDW in 25 μL tubes. The thermal profile was as in Van Essche et al. (). The data were collected during the annealing and extension phase of the reaction. The reactions were performed in a MicroAMP 8-tube stripe (0.2 mL) with MicroAMP Optical 8-Cap Strip (Applied Biosystems). Reactions were run in a 7300 Real-Time PCR System and the data were collected and analyzed by the sds v1.4 software (Applied Biosystems).
Bacteriovoracaceae (Bacteriovorax and Bacteriolyticum)
This series of specific QPCR was performed based on SYBR green I chemistry. The reaction was carried out according to Zheng et al. () using primers Bac 519F and Bac 677R, with some adjustments. Briefly, each 25-μL reaction consisted of 12.5 μL ABsolute Blue QPCR SYBR Green ROX Mix, 0.5 μL of each primer (10 μM), 1 μL of DNA and 10.5 μL of PCR grade DDW. Thermal profiles were 94 °C for 15 min (enzyme activation), and 45 cycles of 94 °C for 30 s, 62 °C for 30 s and 72 °C for 30 s followed by a dissociation stage. The data were collected during the extension phase of the reaction. The reactions were carried out in a 96-well plate (ABGene; Thermoscientific) with an adhesive seal applicator (AB-1391; Thermoscientific) on a Stratagene-MX 3000P Real-time PCR System (Stratagene). For both Bd and Bc reactions, seven standard plasmid samples (9 to 9 × 107 copies per reaction) were included in triplicates in every run, along with three controls without template (NTC).
QPCR for Bacteria
This reaction was performed to measure the total abundance of bacteria in the sample, enabling the calculation of the relative abundance of each BALO group. Standards were prepared by inserting a 1467-bp fragment of the E. coli ML-35 16S rRNA gene amplified with primers 27F and 1494R into a Topo Cloning Vector (Invitrogen). The cloning, transformation and plasmid extraction and quantification were done as above. DNA samples were diluted five times before use in the reactions, which were performed with primers 519F and 907R (Muyzer et al., ; Table S2). Each 25-μL reaction consisted of 12.5 μL of ABsolute Blue QPCR SYBR Green ROX Mix (Thermoscientific) 1 μL of each primer, 1 μL of DNA and 9.5 μL of PCR grade DDW. Reactions were ran at 95 °C for 15 min followed by 40 cycles of 95 °C for 30 s, 58 °C for 30 s and 72 °C for 30 s. Data were collected during the extension phase of the reaction and included 24 standards (three replicates for building a standard curve with copy number of 102–109 per reaction) and six wells without template. The reactions were carried out in a 96-well plate (ABGene; Thermoscientific) with an adhesive seal applicator (AB-1391; Thermoscientific). QPCR was performed on a Stratagene-MX 3000P Real-time PCR System (Stratagene). Total QPCR counts (16S rRNA bacterial genes mL−1) were estimated based on a standard curve with 100% efficiency factoring in the cycle threshold value of the particular sample obtained, taking into account that 25 mL of sample was used for DNA extraction, DNA was eluted in 50 μL of DDW and 1 μL of this DNA was used in the QPCR reaction.
High throughput sequencing
Samples were collected and DNA extracted as above from the water phase of the TF, the DB, and the FBR, in duplicates, in December 2010 and April 2011. The April 2011 sampling set included one sample from the water phase and one from the biofilm phase of the TF. High throughput pyrosequencing (HTS) was performed at the Research and Testing Laboratory (Lubock, TX) on a 454 FLX Genome Sequencer System (Roche Applied Science, IN). All reads were cleaned using the mothur software (Schloss et al., ). First, Fasta and quality data were extracted from the raw SFF file. Sequences were grouped according to barcode and primer, allowing one mismatch to the barcode and two mismatches to the primer. De-noising was achieved with the AmpliconNoise algorithm (Quince et al., ), removing both 454 sequencing errors and PCR single base errors. Next, sequences were trimmed to remove barcode and primer sequences, all sequences with homopolymers (i.e. AAAA) longer than 8 bp, and all sequences < 100 bp long. All sequences were aligned to the silva reference alignment database (Schloss et al., ) and filtered so that they all overlapped perfectly (with no overhang or no-data base pairs). To further reduce sequencing errors, pre-clustering of the sequences was performed, based on the algorithm of Huse et al. (). Finally, chimeric reads were removed with the UChime method (Edgar et al., ), and all chloroplast, mitochondria and ‘unknown’ (i.e. unclassified at the Domain level) reads were deleted.
Sequences have been deposited at the MG-RAST metagenomic database, project #4588 (http://metagenomics.anl.gov/linkin.cgi?project=4588).
Statistical analysis
All 454 reads were clustered into genera based on current RDP-II taxonomy (Cole et al., ). The resulting genera were arranged in a data matrix where each row is a single sample and each column a genus; each data point in the matrix represents the abundance of the particular genus in the particular sample, relativized to the sampling effort (i.e. the number of 454 reads obtained from that sample). Multivariate analysis was performed in pc-ord v6.12 (mjm Software) with Sorensen distances. Ordinations were performed with non-metric multidimensional scaling (NMDS; Mather, ) at 500 iterations, and cluster analyses were performed with flexible beta linkages (β = −0.25). Differences between sample groups were calculated with multi-response permutation procedure (MRPP; Mielke, ), a test based on the assumption that if two groups are different from each other, the average within-group difference is smaller than the average between-group distance. The size of the difference between groups is represented by the A-statistic of the MRPP test, and its significance is identified by the MRPP P-value.
To detect which genera were mainly responsible for differences between groups, we use the method of Dufrêne & Legendre (). The basis for this procedure is the computation of Indicator Values (IVs), which are a combination of the frequency of occurrence and of the abundance of each genus in each group; IV spans between 0 and 100, and is larger if a genus is more frequent and/or more abundant in a given group than in the other group.
Results
QPCR sensitivity
Standard QPCR curves for Bd (Fig. S2a) and Bc (Fig. S2b) were generated by serial dilution of plasmids containing cloned fragments of the 16S rRNA gene amplified from B. bacteriovorus and B. stolpii reference strains for the Bd and Bc reactions, respectively. Reaction efficiency varied between 96% and 100% with correlation coefficients always > 0.99. The minimal detection limit for both Bd and Bc was 10 copies per 5-μL reaction, corresponding to 2 predatory cells mL−1 in the sample, considering that 25 mL of sample was used, DNA was eluted in 50 μL of DDW, 5 μL of this DNA was used per reaction, and there are two 16S rRNA genes in the Bdellovibrio and Bacteriovorax genomes (Rendulic et al., ; Crossman et al., ). However, highly reproducible results were obtained from 90 copies per reaction. Hence, the quantification limit was set at 100 copies per reaction, or 20 predatory cells mL−1 in the sample. Counts of B. bacteriovorus 109J cells growing on an E. coli prey were about 2.5 times higher when measured by QPCR than by counting plaques in a double-layered agar (Fig. S3). This difference may be due to decreased predation in soft agar, to the ability of QPCR to detect predators not forming plaques or to limited detection of small plaques by visual inspection (Van Essche et al., ). The correlation coefficients between copy numbers and cycles of standard QPCR curves for Bacteria were very high (r2 > 0.99).
Quantification of the bacterial community using QPCR
Total bacterial abundance was measured at three time points in the FB, TF and DB in both the FW and SW ZDS. In the former, 7 × 106 to 1 × 108 16S rRNA gene copies mL−1 and in the latter, 4 × 106 to 1.2 × 108 16S rRNA gene copy mL−1 were measured (Table 1). The 16S rRNA gene copy number per genome varies greatly between bacteria, from 1 to 15 (Lee et al., ). Nonetheless, it has been estimated that cells in wastewater plants bear an average of 3.6 copies of the 16S rRNA gene per cell (Dionisi et al., ). Accordingly, cell concentrations in the compartments varied between 1.9 × 106 to 2.8 × 107 cell mL−1, and 1.1 × 106 to 3.3 × 107 cell mL−1 in the FW ZDS and in the SW ZDS, respectively.
Bacterial 16S rRNA gene mL−1 counts in ZDS samples over time
| Compartment | December | April | June |
| FB/FW | 8 × 107 ± 3.1 × 106A1 | 7.3 × 107 ± 1.7 × 106A1 | 1.3 × 107 ± 6.3 × 106B1 |
| TF/FW | 1 × 108 ± 3.1 × 106C2 | 5.7 × 107 ± 7.1 × 106CD1 | 6.7 × 106 ± 9 × 105DE1 |
| DB/FW | 2.2 × 107 ± 5.3 × 105F3 | 3.6 × 107 ± 2.8 × 106F1 | 7.3 × 106 ± 4.9 × 105G1 |
| FBR/FW | 2.6 × 107 ± 5.4 × 1053 | NA | NA |
| FB/SW | 4 × 106 ± 1.4 × 105H4 | 4.2 × 107 ± 2.4 × 106I1 | NA |
| TF/SW | 6 × 106 ± 4.2 × 104J4 | 3.3 × 107 ± 3.4 × 106K1 | 5.5 × 107 ± 4.3 × 106K2 |
| DB/SW | 4.8 × 106 ± 7.2 × 104L4 | 1.2 × 108 ± 1.7 × 107L1 | 3 × 107 ± 1.5 × 106L3 |
| FBR/SW | NA | 7.8 × 107 ± 6.6 × 1061 | NA |
| Compartment | December | April | June |
| FB/FW | 8 × 107 ± 3.1 × 106A1 | 7.3 × 107 ± 1.7 × 106A1 | 1.3 × 107 ± 6.3 × 106B1 |
| TF/FW | 1 × 108 ± 3.1 × 106C2 | 5.7 × 107 ± 7.1 × 106CD1 | 6.7 × 106 ± 9 × 105DE1 |
| DB/FW | 2.2 × 107 ± 5.3 × 105F3 | 3.6 × 107 ± 2.8 × 106F1 | 7.3 × 106 ± 4.9 × 105G1 |
| FBR/FW | 2.6 × 107 ± 5.4 × 1053 | NA | NA |
| FB/SW | 4 × 106 ± 1.4 × 105H4 | 4.2 × 107 ± 2.4 × 106I1 | NA |
| TF/SW | 6 × 106 ± 4.2 × 104J4 | 3.3 × 107 ± 3.4 × 106K1 | 5.5 × 107 ± 4.3 × 106K2 |
| DB/SW | 4.8 × 106 ± 7.2 × 104L4 | 1.2 × 108 ± 1.7 × 107L1 | 3 × 107 ± 1.5 × 106L3 |
| FBR/SW | NA | 7.8 × 107 ± 6.6 × 1061 | NA |
Values shown are means and standard deviations of three biological replicates. Cells within a row connected with different capital letters in superscript are significantly different. Within a column, cells connected with different numbers in subscript are significantly different (P= 0.05).
FW, fresh water; NA, not available for sampling; SW, salt water.
Bacterial 16S rRNA gene mL−1 counts in ZDS samples over time
| Compartment | December | April | June |
| FB/FW | 8 × 107 ± 3.1 × 106A1 | 7.3 × 107 ± 1.7 × 106A1 | 1.3 × 107 ± 6.3 × 106B1 |
| TF/FW | 1 × 108 ± 3.1 × 106C2 | 5.7 × 107 ± 7.1 × 106CD1 | 6.7 × 106 ± 9 × 105DE1 |
| DB/FW | 2.2 × 107 ± 5.3 × 105F3 | 3.6 × 107 ± 2.8 × 106F1 | 7.3 × 106 ± 4.9 × 105G1 |
| FBR/FW | 2.6 × 107 ± 5.4 × 1053 | NA | NA |
| FB/SW | 4 × 106 ± 1.4 × 105H4 | 4.2 × 107 ± 2.4 × 106I1 | NA |
| TF/SW | 6 × 106 ± 4.2 × 104J4 | 3.3 × 107 ± 3.4 × 106K1 | 5.5 × 107 ± 4.3 × 106K2 |
| DB/SW | 4.8 × 106 ± 7.2 × 104L4 | 1.2 × 108 ± 1.7 × 107L1 | 3 × 107 ± 1.5 × 106L3 |
| FBR/SW | NA | 7.8 × 107 ± 6.6 × 1061 | NA |
| Compartment | December | April | June |
| FB/FW | 8 × 107 ± 3.1 × 106A1 | 7.3 × 107 ± 1.7 × 106A1 | 1.3 × 107 ± 6.3 × 106B1 |
| TF/FW | 1 × 108 ± 3.1 × 106C2 | 5.7 × 107 ± 7.1 × 106CD1 | 6.7 × 106 ± 9 × 105DE1 |
| DB/FW | 2.2 × 107 ± 5.3 × 105F3 | 3.6 × 107 ± 2.8 × 106F1 | 7.3 × 106 ± 4.9 × 105G1 |
| FBR/FW | 2.6 × 107 ± 5.4 × 1053 | NA | NA |
| FB/SW | 4 × 106 ± 1.4 × 105H4 | 4.2 × 107 ± 2.4 × 106I1 | NA |
| TF/SW | 6 × 106 ± 4.2 × 104J4 | 3.3 × 107 ± 3.4 × 106K1 | 5.5 × 107 ± 4.3 × 106K2 |
| DB/SW | 4.8 × 106 ± 7.2 × 104L4 | 1.2 × 108 ± 1.7 × 107L1 | 3 × 107 ± 1.5 × 106L3 |
| FBR/SW | NA | 7.8 × 107 ± 6.6 × 1061 | NA |
Values shown are means and standard deviations of three biological replicates. Cells within a row connected with different capital letters in superscript are significantly different. Within a column, cells connected with different numbers in subscript are significantly different (P= 0.05).
FW, fresh water; NA, not available for sampling; SW, salt water.
Abundance and dynamics of specific BALOs in the ZDS
Bd populations were quantified by QPCR in both the FW and the SW systems, over 7 months. Bd was found in all the compartments of both systems but the predator was more abundant in the FW ZDS than in the SW ZDS at most sampling times, ranging from 1.5 × 103 to 4.4 × 105 cell mL−1 and from 1.5 × 102 to 4.9 × 103 cell mL−1, respectively (Table S3). The presence of Bd in the SW ZDS is remarkable, as to date this genus was thought to be confined to non-saline environments (Baer et al., ; Davidov & Jurkevitch, ). Bd abundance within compartments in each system fluctuated significantly with time. Noteworthy, Bd seasonal dynamics in the FB and TF were similar, differing from those of the DB (Fig. 1a, b and d). Additionally, whereas FB and TF populations in the FW ZDS fluctuated by less than six-fold during the sampling period, in the SW ZDS they increased with time by more than 30-fold (Fig. 1a and b; Table S3). A biofilm that developed on PVC clips introduced in the FW TF supported high Bd populations: Bd abundance calculated assuming a depth of 20 μm deep was 1.3 × 107 cell mL−1, but reached 1.2 × 104 cell mL−1 (difference is significant at P= 0.001) in samples from the planktonic phase of the same compartment taken at the same time.
Bdellovibrio (a, b) and Bacteriovorax (c) abundance (cell mL−1) in the fresh water (FW) and salt water (SW) ZDS, as measured by clade-specific QPCR targeting the 16S rRNA gene, over time. (d) Cluster analysis relating the Bray–Curtis similarity between samples, using unweighted pair-group average (upgma) clustering. SE, standard error of mean.
Populations of Bc predators were more alike and stable than the Bd populations in both the FW and SW systems, and between the compartments in each system, spanning from 1.2 × 104 to 1.8 × 105 cell mL−1 (Fig. 1c, Table S4). Overall, and in both systems, Bc populations were almost always more abundant than Bd populations. The restricted number of FBR samples precluded comparative analysis but suggested that Bd population levels were similar to those of the other basins, whereas those of the Bc were lower (Fig. 1).
Correlation between specific BALO abundance, biotic and abiotic parameter
The relative abundance of the detected BALOs (Bd + Bc) in both ZDS was separately calculated as the 16S rRNA gene count ratio of Bd and of Bc over total Bacteria. The total of the targeted BALOs accounted for 0.07–0.78% of the total bacterial 16S rRNA gene count, or 0.13–1.4% of the total cell counts, using 2 and 3.6 gene copies genome−1 for BALOs and for total bacteria, respectively. A significant correlation between BALOs and the total bacterial counts was found in the FW ZDS (r2 = 0.78, P = 0.012) but not in the SW ZDS (r2 = 0.51, P= 0.19), hinting at direct coupling between predator and prey abundance in the former system. The relative abundance of BALOs (Bd + Bc) was highest in most of the DB samples in both the FW and SW systems, reflecting the high population levels of predators reached in this compartment where heterotrophic degradation of organic matter occurs (Figs 1 and 2). A significant correlation was found between Bd abundance and reactive phosphate in the SW system (r2 = 0.92, P < 0.0024). Since no information on phosphate thresholds in Bd exists, the reason for this correlation is unclear; it is noteworthy that sediment-trapped phosphate was released into the water starting about March. No other correlation between the targeted BALO abundance and temperature, oxygen, pH, nitrates, nitrites, ammonia, or sulfides was detected.
Percent Bd (a) and Bc (b) and total BALO (Bd + Bc) (c) from the total Bacterial community as determined by clade-specific (Bdellovibrio-Bd, Bacteriovorax + Bacteriolyticum-Bc) and general Bacteria QPCR targeting the 16S rRNA gene, over time. FW, fresh water; SE, standard error of mean; SW, salt water.
Sequence-based community analysis
The bacterial composition of the FW ZDS, which includes the populations preyed upon by BALOs, was explored by HTS of the 16S rRNA gene. In total, 74 003 curated sequences were obtained at two time points, and from 10 samples: c. 6000–9600, 7500–9300, 5500–7100, and 5400 sequences from four, three, two, and one samples from the DB, the TF, the FBR, and the TF biofilm, respectively. The diversity was essentially covered, reaching 95–98% (Table S5), and totalling 238 OTUs. The DB was richer (125 OTUs) and the TF had fewer of OTUs (89). The structure of the bacterial communities did not differ between the different compartments, as revealed by NMDS and MRPP tests (P > 0.1; Fig. 3). However, a seasonal effect was clearly seen, with April and December samples clustering separately (MRPP A = 0.22, P= 0.002; Fig. 3), a separation contributed by 7 and 8 OTUs, respectively (Table S6).
NMDS ordination of 16S rRNA gene sequences from December (▲) and April (△) samples. Convex hulls connect samples from each sampling time. Ordination axes do not possess biological meaning; stress = 8.4%, confirming a good correlation with the data.
The most abundant four OTUs constituted over 60% of the total sequences and belonged to: Flavobacteria (29.7 ± 12.9%), which strongly dominated all samples; unclassified Bacteria (12.8 ± 7.1%), which were mainly found in the BF, the DB and the FBR, and significantly less in the TF; Bacteroidetes (10.4 ± 3.2%), which were more abundant in the BF, DB, and FBR, and appeared to decrease in the TF; and Fusobacteriaceae (8.6 ± 5.0%), which were less common in the BF than in the other samples (Table S7). The 20 most commonly retrieved OTUs, representing 85.4 ± 4.7% of all sequences, included taxa commonly found in freshwater [Flavobacteriaceae, Cytophagaceae, Polynucleobacter (Burkholderiaceae)], Lacibacter (freshwater fish) (Cetobacterium – Fusobacteriaceae; Tsuchiya et al., ), Bacteroidetes and Flavobacteria (Kim et al., ; Sullam et al., ); and wastewater and sludge [Verrucomicrobia, Flavobacteriaceae, Novosphingobium (Sphingomonadaceae), Burkholderiales, Rhizobiales; Zhang et al., ]. Thus, the composition of the ZDS microbiota appears to reflect its dual purpose. Confirming this assumption, a phylogenetic analysis of Flavobacteria, the system's dominant phylum, revealed that the large majority of these sequences (79%) are more closely related to those of FW fish Flavobacteria than to those from other habitats (Fig. S4).
To further characterize the ZDS, potentially important ecological functions, that is nitrification, sulfate reduction, and predation were assigned, based on the identification of the populations in the different compartments (Table S8). The ammonium oxidizers Nitrosomonadaceae were barely detectable in the FBR, and only sporadically in the DB and TF, whereas the nitrite oxidizers (Nitrospira) formed dominant populations (c.1–2% of total) in the FBR, and were also abundant in the BF. Sulfate reducers were only detected in the DB and at a lower concentration in the FBR.
The estimated fraction of BALO sequences detected by HTS was between two-thirds to a fifth of that found by QPCR, the relatively low number of sequences inducing a potentially large error (Fig. 2, Table S8). BALOs constituted 20–40% of the potential Deltaproteobacteria predators (Myxococcales); other potential predators included Herpetosiphon, which were present in the BF and DB, and Ensifer (Jurkevitch & Davidov, ), present in part of the DB and FBR samples. Some Flavobacteria strains have been shown to be predatory (Banning et al., ) but no 16S rRNA gene sequences closely related to these strains were found in our dataset.
A phylogenetic analysis of all the BALO sequences was performed (Fig. 4). The restricted length of the sequenced fragments precluded a precise analysis at the species level (Aharon et al., ) in both the Bdellovibrionaceae and the Bacteriovoracaceae. Yet, it could still be discerned that all Bdellovibrio sequences were similar to strain B. bacteriovorus HD100 but that no sequences related to the epibiotic predators B. exovorus or Micavibrio were found.
Phylogenetic tree of the BALO in the ZDS. All 16S rRNA gene sequences of BALO in the ZDS were compared with the BALO type-strain sequences retrieved from GenBank. The tree was constructed by maximum likelihood based on the Tamura–Nei model. Bootstrap consensus was inferred from 100 replicates (bootstrap values are shown next to the branches). Branches corresponding to partitions reproduced in < 50% bootstrap replicates were collapsed. BF, biofilm; FBR, fluidized bed reactor.
Discussion
Predation is a major ecological process, affecting many parameters of ecosystem structure and function (Huntley & Kowalewski, ; Rodriguez-Valera et al., ; Vasseur & Fox, ; Vos et al., ; Estes et al., ). At the microbial level, bacteriophages and protists have been shown to make major contributions to bacterial mortality and turnover in aquatic environments (Fuhrman, ; Samuel et al., ); yet, our knowledge of bacterial predation by bacterial predators is still very fragmentary. Here, we analyzed the dynamics of defined BALO populations and the bacterial communities within which the predators are embedded, providing the first quantitative data for the BALO community, of its constituent populations and of their relative abundance as affected by salinity and other environmental parameters.
The relative abundance of BALOs as measured by QPCR ranged between 0.1% and 1.4% of the total bacterial population, suggesting that BALOs may be considered a relatively ‘abundant species’. The fraction of sequences obtained in the HTS analysis was at the lower end of this distribution, reaching 0.059% after curation, or about 0.1% as expressed on a cell basis (based on two copies of the 16S rRNA gene in the BALOs and 3.6 in the general population, see ), a bias possibly due to the relatively small number of sequences obtained (42 in total).
A major difference between Bdellovibrionaceae and Bacteriovoracaceae is the environments they colonize. The former has so far been found in freshwater (including wastewater), soil and sometimes in animals (Baer et al., ; Davidov & Jurkevitch, ) and it cannot withstand salinity above 5 p.p.t. (Varon & Shilo, ; Baer et al., ); the latter includes the halophilic Bacteriovorax, which populate brackish and saline to hypersaline water, and the freshwater- and soil-dwelling Peredicter and Bacteriolyticum (Baer et al., ; Pineiro et al., ). Remarkably, Bdellovibrionaceae (Bd) were detected in the SW ZDS (salinity 20 p.p.t., or about two-thirds seawater), indicating that the ecological range of this taxon may be more extended than previously thought. This is further supported by a recent screen of published 16S rRNA amplicon datasets of environmental HTS studies that showed the presence of B. bacteriovorus-related sequences in marine habitats (Dolinšek et al., ) but a full demonstration that Bd indeed populates saline water requires isolation. Bdellovibrionaceae and Bacteriovoracaceae may thus co-exist, sharing many environments, including saline ones (Davidov et al., ). Differential utilization of prey within and between clades may enable this coexistence: Chen et al. (, ) elegantly showed that when particular prey are spiked into natural samples, the abundance of specific Bacteriovorax clusters is consequently affected. Strikingly, specific Bacteriovorax clades were differently amplified by spiking with prey representing freshwater vs. prey representing saline water habitats. Altogether, most Gram-negative bacteria in the environment may be susceptible to local BALOs. Rice et al. () found that about 80% of Gram-negative isolates from an estuarine environment were susceptible to BALOs, including some Flavobacteria, the dominant taxon in the FW ZDS.
BALOs seem to differ widely in the breadth of their prey ranges. Some appear to be very restricted, preying on few different bacterial types, whereas others appear to be rather generalist predators (Chen et al., ; Koval et al., ; Chanyi et al., ). Since these BALOs can share the same habitats (Jurkevitch et al., ; Chen et al., ), the dynamics of the populations may not follow a ‘kill the winner’ pattern as observed with phages in wastewater (Shapiro et al., ). It has been shown that closely related (within families and even within species) phylogenetic clusters can form ecologically distinct populations with seasonally or spatially segregated populations. As an example, specific marine Vibrio ecotypes are free-living or are associated with zooplankton-free particles, whereas others are found with zooplankton and large particles (Hunt et al., ). Co-existence of related populations within the same ecosystem may be mediated by resource partitioning, a phenomenon observed not only in the water column (Hunt et al., ) but also in soil and in the phyllosphere (Wilson & Lindow, ; Oger et al., ). In the case of BALOs, the main partitioning resource is certainly prey. Different prey ranges may lead to sympatric separation of the predators' populations that may be reinforced as the prey themselves distribute differentially in the environment. In soils, the microbiota of the plant rhizosphere and of the bulk soil differs, selecting for BALO subpopulations exhibiting different ribotypes; these populations may further differ in prey preference (Jurkevitch et al., ). Thus, the physiological capacities to use particular resources (e.g. plant exudates or suspended organic matter) may result in resource partitioning of consumer populations, leading to spatially segregated (potential) prey populations, which in turn affect BALO distribution and dynamics in the environment.
Correlative data between BALO community structure and environmental parameters have so far yielded few additional insights beyond the effects of salinity and temperature (Kelley et al., ; Pineiro et al., ). Our data suggest that water temperature (20 ± 1 °C in winter to 24 ± 1 °C in spring) had no direct effect on ZDS BALO populations as it does in natural ecosystems subjected to larger fluctuations (Kelley et al., ) and, for reasons that are not clear, only phosphate was found to correlate with Bd in the SW ZDS. In contrast, the general bacterial community was affected by season (Fig. 3), possibly due to temperature variations or to changes in stock density (fish biomass increased with time), a variable known to affect the microbiota (Zhou et al., ). The environmental conditions observed in the ZDS compartments are brought about by engineering, physically designing the specific habitats that enable different microbial activities (Cytryn et al., , ; Gelfand et al., ) while connecting them with a constant water flow. In ZDS, mixing may then be sufficient to blur local ecological selection. Indeed, the overall microbiota community structure of the FW ZDS remained similar in all basins (Fig. 3); however, the dominant taxa reflected the dual use of the recirculating system as a fish culture and wastewater treatment plant (Table S7; and see ). No differential distribution of specific BALOs between pools was found based on the short pyrosequencing reads obtained. Deeper sequencing with longer reads, functional analysis of BALO activities in the different components of the system, and a detailed examination of spatial distribution (e.g. biofilm and flocs vs. planktonic) may yet reveal more complex patterns of abundance. A striking example of such specific predator –prey interactions in confined niches was discovered by Dolinšek et al. (): Micavibrio-like cells that were barely detected in the planktonic phase of a nitrifying bioreactor were specifically associated with sub-lineage I Nitrospira clusters but with no other bacteria. Although only about 4% these Nitrospira aggregates were surrounded by Micavibrio-like bacteria, it was posited that the predators may prevent any particular nitrifier strain from becoming dominant, thus driving diversification and functional redundancy, and increasing resistance to perturbations. In the FW ZDS, Nitrospira sequences were quite abundant in all samples (Table S7) but Micavibrio-related sequences were not found, and it is not known whether other BALOs prey upon Nitrospira. Therefore, while the impact of BALOs on their potential prey populations as a whole may be relatively minor (Data S1), the effects of BALO predation may be significant when spatial structuring is taken into account. The impact of bacterial predation may not be restricted to BALOs. HTS data suggest that a potential and diverse bacterial predatory community may be at least two to five times more abundant than that of BALOs; these other predators exert predation through various strategies (Jurkevitch & Davidov, ). Consequently, bacterial predation by bacteria may affect the microbiota structure and dynamics in many intricate and unknown ways. As recently described by Yang et al. (), when intraspecific variation in the predatory community is taken in account, numerical and functional responses in predator–prey systems do not maintain proportionality.
Acknowledgements
This work was supported by grants of the Israel Science Foundation (No. 1290/08 and 1583/12). The authors confirm that they have no conflicts of interest.
References
Author notes
Editor: Gary King
Israel Science Foundation
1290/08
1583/12
