Changes in Microbiome Activity and Sporadic Viral Infection Help Explain Observed Variability in Microcosm Studies.

The environmental conditions experienced by microbial communities are rarely fully simulated in the laboratory. Researchers use experimental containers ("bottles"), where natural samples can be manipulated and evaluated. However, container-based methods are subject to "bottle effects": changes that occur when enclosing the plankton community that are often times unexplained by standard measures like pigment and nutrient concentrations. We noted variability in a short-term, nutrient amendment experiment during a 2019 Lake Erie, Microcystis spp. bloom. We observed changes in heterotrophic bacteria activity (transcription) on a time-frame consistent with a response to experimental changes in nutrient availability, demonstrating how the often overlooked microbiome of cyanobacterial blooms can be altered. Samples processed at the time of collection (T0) contained abundant transcripts from Bacteroidetes, which reduced in abundance during incubation in all bottles, including controls. Significant biological variability in the expression of Microcystis-infecting phage was observed between replicates, with phosphate-amended treatments showing a 10-fold variation. The expression patterns of Microcystis-infecting phage were significantly correlated with ∼35% of Microcystis-specific functional genes and ∼45% of the cellular-metabolites measured across the entire microbial community, suggesting phage activity not only influenced Microcystis dynamics, but the biochemistry of the microbiome. Our observations demonstrate how natural heterogeneity among replicates can be harnessed to provide further insight on virus and host ecology.

Introduction

Research in environmental microbial ecology is often driven by observations and survey-style data where samples are gathered from natural systems during transects or time series to provide important insight into the spatial and temporal diversity and function of microbial communities (e.g., Venter et al., 2004). Yet, in many cases, the ability to test hypotheses (Popper, 1959) requires experimental manipulation of the entire community under natural conditions. To do this, aquatic researchers often employ “bottle experiments” where water samples are collected and distributed across a series of incubation vessels. Bottle experiments have been used for many decades (Schelske, 1984) and are advantageous because the duration of the experiments is similar to the generation time of the organisms being studied (Sterner, 2008). Incubation vessels, or bottles, are manipulated experimentally by introducing biological, chemical, or physical variation to a subset of the bottles that can be compared to unamended “control” bottles. Many experiments of this nature have been completed and used to document processes that range from the role of iron (Fe) as a growth-limiting element in marine systems (Martin and Fitzwater, 1988; Hutchins et al., 1999) to the rate at which virus populations are actively produced within microbial communities (Wilhelm et al., 2002). In fresh waters, these types of experiments have been used extensively to resolve the effects of nutrient additions on planktonic communities (Wilhelm et al., 2003; DeBruyn et al., 2004; Xu et al., 2010). Despite their detractors (Carpenter et al., 1995; Schindler, 1998), incubation experiments remain fundamental to the aquatic microbial ecologists’ toolset.

The process of executing a well-designed experiment in an enclosure, such as the common bottle-based experiments, is challenging and requires that attention be paid to spatial and temporal scales that could influence the outcome of the experiment and the interpretation of the results. The simple act of placing water containing a natural phytoplankton community into an enclosed bottle can result in differences which are often considered to be “environmental uncertainty” or “bottle effects” (Stewart et al., 2012). Variability due to the limited volume of a container has been well documented: the simple act of controlling temperature and light across different container sizes is challenging (Spivak et al., 2011). Enclosed bottle experiments allow for independent successional patterns and levels of heterogeneity, even amongst replicate bottles of the same treatment. Additionally, unlike natural systems, regenerative processes, inputs, and sinks for nutrients and dissolved gases are often artificially restricted in a microcosm (Zobell and Anderson, 1936), with cascading nutrient limitations as a result (Davidson and Howarth, 2007). The presence of a few large grazers in bottle experiments can also introduce noise into the data (e.g., Vanni and Temte, 1990; Urrutia-Cordero et al., 2015). Hypothetically, when water containing natural plankton communities are incubated in several bottles under identical conditions, the response metric measured at the end of the experiment should be consistent among replicates; however, variation is always observed, but not always explained (e.g., see references above).

While many researchers focus only on phytoplankton responses in bottle experiments, it is well-established that heterotrophic bacteria in the phycosphere profoundly impact phytoplankton (Cook et al., 2020; Jankowiak and Gobler, 2020). Understanding the responses and variabilities associated with heterotrophic bacteria placed in a bottle may help explain variation, because slight changes to the microbiome can lead to notable differences in microbial diversity and function over time: even when an experiment starts from a consistent microbiome, it can quickly change identity and function (Pound et al., 2021a).

Perhaps more challenging is resolving the influence of viruses, which are ubiquitous in natural systems and key drivers of processes at the biogeochemical, population, and ecosystem levels (Wilhelm and Suttle, 1999; Breitbart, 2012; Sullivan et al., 2017). Unlike copepods and other grazers, viruses are often “invisible” to researchers (Roux et al., 2015) and thus are overlooked as possible drivers of biological variability (Pound et al., 2020). Yet when viruses in a system actively infect a population, the resulting cell lysis (Fuhrman, 1999; Wilhelm and Suttle, 1999) and/or viral-induced reshaping of gene expression can have strong community scale effects (Lindell et al., 2005; Forterre, 2011). Viruses have been shown to be constant companions in Microcystis-dominated communities (Yoshida M. et al., 2008; Yoshida et al., 2010) and to play potentially important roles in system ecology (Steffen et al., 2017).

To investigate the effects of short-term nutrient pulses during Microcystis blooms in the Laurentian Great Lakes, we conducted in situ bottle incubations comparing unamended controls to nutrient amendments over 48 h. We examined the transcriptome and metabolome of the entire community in independent replicates with the aim of identifying community responses to episodic nutrient introductions. During our analyses, we noted variability between treatments as well as between replicates within treatments. We examined this variation in an attempt to categorize outlier observations as due to either “unexplained random error” or “biological activity.” Our observations demonstrated not only how nutrient additions influenced transcriptional and metabolic landscapes of a toxic cyanobacterial bloom, but also how the confluence of experimental effects and natural variation can be resolved, including the variation arising from the influence of viruses on microbial communities.

Materials and Methods

Experimental Setup for Incubations and in situ Measurements

Lake water was collected in 20 L carboys from the surface of a Microcystis spp.- dominated bloom in Lake Erie, (United States) on July 21st, 2019 (41° 44.946′ N 83° 06.448′ W). Water column physiochemistry was recorded prior to sampling using an EXO multiparameter sonde (YSI xylem). Homogenized water was aliquoted into 18 × 1.2 L polycarbonate bottles and subject to four nutrient addition treatments in biological triplicate: 1 μM phosphate as KH2PO4:K2HPO4 (Belisle et al., 2016), and 180 μM nitrogen as either stable-isotope-labeled nitrate (15NO3^–), stable-isotope-labeled ammonia (15NH4^+), or stable-isotope dual-labeled urea (13CH415N2O). Non-replicated bottles of non-labeled nitrate, ammonia, and urea additions were also included to act as standards for metabolome normalization (see metabolite methods below). Six bottles were included as controls with no nutrient additions, three of which were sampled immediately and considered time zero samples. The remaining three control bottles and all nutrient-amended bottles were incubated in situ in Lake Erie off docks at Stone Laboratory (41° 39.467′ N 82° 49.600′ W) for 48 h.

Samples were collected at 48 h for metabolomics analysis and RNA sequencing, as well as to measure pH, nutrients in the dissolved fraction, chlorophyll a, and toxins. Total sample collection time took ∼ 90–120 min. Approximately 180 mL of water from each bottle was collected on 0.2-μm GTTP filters (EMD Millipore Corporation, Burlington, MA, United States), stored in 2-mL cryovials and flash-frozen in liquid nitrogen for metabolomics analysis (see below). Approximately 150 mL of water from each bottle was filtered through a 0.2-μm nominal pore-size Sterivex™ filter (EMD Millipore Corporation, Burlington, MA, United States) and flash-frozen in liquid nitrogen for subsequent RNA extraction and sequencing (see RNA extraction methods below). pH was measured via immediate readings of 15 mL subsamples using a pH probe (Mettler Toledo Seven Compact™ pH/Ion meter S220, fitted with a Mettler InLab Expert Pro-ISM electrode with a temperature range of up to 100°C). Approximately 50 mL of 0.2-μm-filtered water was captured from the Sterivex filtrate for dissolved nutrient analysis (nitrate, ammonia, and phosphate) and stored at 4°C until processing on a SEAL Analytical (Mequon, WI, United States) QuAAtro 5-Channel continuous segmented flow auto-analyzer (Chaffin et al., 2019). Approximately 100 mL of water was filtered onto 0.2-μm nominal pore-size 47-mm diameter polycarbonate filters, and chlorophyll a extracted in 90% acetone for 24 h at 4°C prior to measurement on a Turner Designs 10-AU Field Fluorometer (Welschmeyer, 1994). Intracellular microcystin concentrations were determined using HPLC coupled with single quadrupole mass spectroscopy and photodiode array spectroscopy as described by Tang et al. (2018) and Boyer (2020). Standards (microcystin-RR, -LR and -LF) were analyzed at the beginning and end of each run to ensure that the retention times did not drift, and individual toxin congeners (RR, YR, and LR) were identified based on retention time, characteristic absorbance spectrum in the photodiode array detector, and characteristic molecular ions.

RNA Extraction, Processing, and Sequencing

RNA extraction, sequencing, and quality control details can be found in previous publications (Tang et al., 2018; Pound et al., 2020). A step-by-step protocol describing the acid-phenol-based RNA extraction method with an added DNase treatment to remove any residual DNA is available (Pound and Wilhelm, 2020a). RNA quality was checked using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Walthan, MA, United States) and quantified using the Qubit hsRNA assay (Invitrogen, Waltham, MA, United States). Extracted RNA was processed using an Illumina ® Stranded Total RNA Prep, Ligation with Ribo-Zero Plus (for 96 Samples) and then 50-million 100-bp paired-end reads were generated on the Illumina NovaSeq platform at Hudson Alpha Discovery Life Sciences (Huntsville, AL, United States). Raw reads are available on the NCBI SRA database under BioProject number PRJNA737197. Sequence processing and assembly details can be found in an online, step-by-step protocol (Pound and Wilhelm, 2021b). Briefly, sequences were quality controlled and trimmed in the CLC Genomics workbench version 20.0.4 (Qiagen, Germantown, MD, United States). Lingering ribosomal rRNA reads were removed in silico using SortMeRNA version 4 (Kopylova et al., 2012). The non-ribosomal trimmed reads from all samples were jointly assembled in MegaHit version 1.2.9 (Li et al., 2015). As previously described (Pound et al., 2020), this was done to reduce redundancy among identical sequences in various samples.

Characterization of Microbiome in Metatranscriptome Assemblies

Gene coding sequences were predicted from nucleotide and translated amino acid sequences using MetaGeneMark version 3.25 (Besemer and Borodovsky, 1999; Zhu et al., 2010). Amino acid sequences were functionally and taxonomically annotated via KEGG orthology using the prokaryote + eukaryote + virus database in GhostKOALA (Kanehisa et al., 2016; Pound et al., 2021b). KEGG orthology numbers (KOs) were assigned to predicted coding sequences and transcript expression was estimated by recruiting trimmed reads to the coding sequences using a 90% similarity fraction over a 90% length fraction in CLC Genomic Workbench. Reads that could be recruited to more than one contig with equal identity were randomly assigned. The KEGG orthology annotations provided by GhostKOALA allow for estimates of individual taxa and/or individual annotated genes or functional groups. DNA-dependent RNA polymerase (rpoB) expression was used to quantify active transcription as a proxy for organism abundance. Differential expression analysis of Microcystis-specific functional genes annotated with KOs was performed using DESeq2 (Love et al., 2014).

Characterization of Virus Signatures in Metatranscriptome Assembly

Identification and transcript quantification of viruses of interest was performed using a BLASTx method previously described (Pound and Wilhelm, 2019, 2020b). Phage protein databases were established for the following genes of interest: terminase, major capsid protein, and tail sheath with all reference sequences downloaded from UniProt KB database version 2020_03 (see data availability). The assembled contigs from all samples were then queried using the command-line BLASTx against each database specific to a gene of interest, retaining sequences with a minimum contig length of 300 bp and an e-value < 1 × 10–20 (Camacho et al., 2009). Contigs with BLASTx hits were considered “candidates” and only the aligned portion of the gene sequence was used for subsequent expression analyses (Pound and Wilhelm, 2019). Non-ribosomal trimmed reads from each sample were recruited to aligned portion of each candidate contig in CLC Genomics workbench version 20.0.4 (Qiagen) using a 90% length fraction and 90% identity (Pound and Wilhelm, 2020b). Non-specific recruitments (reads that could be recruited to more than one contig with equal identity) were randomly assigned, to ensure all reads were quantified. Phylogenetic trees were established with maximum likelihood phylogenies based on reference proteins from isolated organisms using PhyML version 3.0 (Guindon et al., 2010). To determine taxonomy, candidate sequences were placed on reference phylogenetic trees using the pplacer algorithm (Matsen et al., 2010). Transcript expression for the entire Microcystis phage genome was estimated by recruiting reads to the Microcystis phage Ma-LMM01 genome (Accession AB231700.1) using a 90% length fraction and 90% identity (Yoshida T. et al., 2008; Pound and Wilhelm, 2020b).

Metabolite Extraction and Assessment

Samples for metabolomics analysis were processed at the Biological and Small Molecule Mass Spectrometry Core (BSMMSC), University of Tennessee, Knoxville, TN (RRID: SCR_021368) in a manner similar to Krausfeldt et al. (2019). In brief, filters with biomass were exposed to pre-chilled extraction solvent (1.3 mL of 40:40:20 HPLC grade acetonitrile/methanol/water with 0.1 M formic acid) at −20°C for 20 min before water soluble metabolites were extracted from filters at 4°C. Extracted metabolites were dried under a stream of N2 and resuspended in HPLC grade water. Samples were immediately placed in an Ultimate 3000 RS autosampler (Dionex, Sunnyvale, CA, United States) and an injection volume of 10 μL was separated through a reverse-phase C18 (Synergi 2.5 μm Hydro- RP, 100 Å, 100 × 2.00 mm) liquid chromatography column (Phenomenex, Torrance, CA, United States) kept at 25°C. Solvent A consisted of 97:3 water:methanol, 10 mM tributylamine, and 15 mM acetic acid. Solvent B was 100% methanol. The solvent gradient from 0 to 5 min was 100% A: 0% B, from 5 to 13 min was 80% A: 20% B, from 13 to 15.5 min was 45% A: 55% B, from 15.5 to19 min was 5% A: 95% B, and from 19 to 25 min was 100% A: 0% B with a flow rate of 200 μL/min. The mass spectrometer was operated in full scan mode following an adopted protocol (Lu et al., 2010). The chromatographic eluent was ionized via an electrospray ionization (ESI) source in negative mode and coupled to an Exactive Plus Orbitrap mass spectrometer (Thermo Fisher Scientific, Waltham, MA, United States) through a 0.1-mm internal diameter fused silica capillary tube. The samples were run with a spray voltage of 3 kV, N2 sheath gas of 10 (arbitrary units), capillary temperature of 320°C, and automatic gain control (AGC) target set to 3 × 106 ions. Samples were analyzed at a resolution of 140,000 and a scan window of 85–800 m/z from 0 to 9 min, and 110–1000 m/z from 9 to 25 min. Files generated by Xcalibur (RAW) were converted to the open-source mzML format (Martens et al., 2011) via the open-source msconvert software as part of the ProteoWizard package (Chambers et al., 2012). MAVEN (mzroll) software (Melamud et al., 2010; Clasquin et al., 2012), which uses a peak grouping algorithm and retention time alignment, was used for visualization of extracted ion chromatograms. Metabolites, validated from an in-house library of approximately 300 compounds, were first manually identified using samples from the unlabeled nutrient bottles and integrated based on exact mass (± 10 ppm mass tolerance) and retention time (Δ ≤ 1.5 min). Peak intensities for a given sample were normalized by the sum of all metabolites detected for that sample. However, only metabolites that showed stable isotope incorporation were considered in further analysis, ensuring that the metabolites considered were being actively used or produced by the cell at the time of sampling. The metabolomics profiling has been submitted to MetaboLights database[1] under study MTBLS3715.

Statistical Analyses

All statistical analyses were carried out in R studio (Rstudio_Team, 2019). Difference between treatments for chlorophyll a, toxin concentrations, transcript expression, and metabolite peak intensities was determined using ordinary one-way ANOVAs with a Tukey’s HSD test. A significance level of α = 0.05 was used in all analyses. Shifts in microbial communities for phototrophs and heterotrophs were characterized using canonical correspondence analyses (CCA) that incorporated dissolved nutrient metadata. To better characterize the co-occurring community, Microcystis was not included in this analysis. Phylogenetic groups were classified as primary or secondary producers manually. The CCAs were performed using the Vegan package (v. 2.5-7) with the species scores scaled by eigenvalues (Okansen et al., 2020). Pearson correlation coefficients were established the Hmisc package (v. 4.4-2) (Harrell, 2020) with Benjamini–Hochberg corrections for multiple comparisons using the R stats package (v. 3.6.1).

Results

Basic Bloom Characterization

Water column physiochemistry at the onset of the experiment included dissolved oxygen = 8.73 mg L–1; pH = 8.89; turbidity = 0.63 NTU; temperature = 26.2°C; and chlorophyll a = 0.23 RFU. After 48 h of incubation, chlorophyll a concentrations were higher in phosphate-amended bottles (p-value < 0.001) than in the control bottles (Figure 1 and Supplementary Data 1). Dissolved nutrients revealed a significant drawdown in nitrate in the phosphate-amended bottles (p-value < 0.001), as compared to the control bottles (Supplementary Data 1).

FIGURE 1

Chlorophyll a concentration across experimental manipulations. Individual bottles are denoted in green and the average is denoted in black. Error bars represent standard deviation and star (*) indicates significant difference from control (p < 0.05).

Microcystis Nutrient Transport Gene Expression

As expected, Microcystis was the most represented genus in the sequencing dataset, with normalized DNA-dependent RNA polymerase (rpoB) expression ranging from 5.2 to 27 per 107 base pairs (Supplementary Figure 2). Overall, ammonium (amt) and urea transporters (urt family) were significantly upregulated relative to the control in the phosphate-amended bottles, while phosphate transporters (pst family) were significantly downregulated (p-value < 0.01) (Supplementary Data 1). Relative to control, nitrate/nitrite transporters (nrt family) were downregulated in the ammonia- and urea-amended bottles (p-value < 0.001). Within control bottles, there were significant differences in nitrate and urea transporter expression relative to the T0 samples. Both transporters were significantly upregulated in the control incubations relative to T0 (p-values < 0.001) (Supplementary Data 1).

Co-occurring Host Community Gene Expression

The sum of rpoB expression for all other co-occurring, non-viral organisms in the community ranged from 5.6 to 13 per 107 base pairs (Supplementary Figure 2). In most treatments, Microcystis alone showed higher transcript expression levels than the rest of the community combined, particularly in the phosphate treatments. Ammonia was the only treatment in which the combined community transcript expression of rpoB was higher than Microcystis. Changes in the co-occurring community appeared to be strongly driven by shifts in transcription within the heterotrophic community. This was particularly apparent in the transcriptional shifts in the heterotrophic community observed between T0 and controls (Supplementary Figure 3). We note that this shift was not observed in the ammonia bottles, with the heterotrophic community demonstrating a transcriptional composition more similar to the T0. Characterization of the taxonomy within these shifts confirmed that they primarily occurred in the heterotrophic community, with lesser change in the cyanobacterial species. However, there was notable changes in cyanobacteria between T0 and P-amendment bottles (see below). The majority of active expression (75–85%) in the co-occurring community, based on the sum of all samples, originated from Actinobacteria, Bacteroidetes, Proteobacteria, and other Cyanobacteria species (Figure 2). Excluding Microcystis, Bacteroidetes showed the highest levels of active transcription in the T0 bottles, while Betaproteobacteria was the most transcriptionally active in the ammonia and nitrate bottles. Within the cyanobacterial group (excluding Microcystis), Pseudanabaena was the most transcriptionally active in all but the ammonia bottles, where Synechococcaceae were the most represented (Supplementary Data 1). Variations in dominant co-occurring taxa occur not only between nutrient treatments, but among bottles within the same treatment. For example, there is notable variability in Betaproteobacteria between replicates. Even when Microcystis remained fairly constant between replicates of a treatment (e.g., ammonia, nitrate, and urea additions), the co-occurring microbial community did not (Figure 2).

FIGURE 2

Heatmap of normalized expression of the DNA-dependent RNA polymerase (rpoB) of dominant phyla (top three) and classes (bottom two) in co-occurring heterotrophic microbial community. Expression was normalized to the library size for each sample. Histograms at the top describe the transcriptional response of Microcystis spp. rpoB for comparison.

Viral Community Characterization

An evaluation of the viral community using a BLASTx approach revealed that Microcystis phage were the most transcriptionally active viruses and that viruses infecting other community members were either largely absent or were not actively transcribed. Normalized expression computed from either Microcystis phage contigs (data not shown) or from the entire Microcystis phage Ma-LMM01 genome indicated low levels of phage activity (2.2–57 per 108 base pair) in most bottles, excluding the phosphate amendments. However, across phosphate treatments, phage expression was highly variable, ranging from 4.7 to 48 per 107 base pair (Figure 3). The third phosphate replicate had a 10-fold increase in phage transcript expression compared to the first phosphate replicate. These increases were not associated with changes in host expression, as Microcystis rpoB transcription remained unchanged between +P and other treatments.

FIGURE 3

Heatmap of normalized expression for Microcystis phage Ma-LMM01 genes in individual experimental bottles. Histograms describe the sum of normalized expression for each bottle (horizontal histogram) and each gene (vertical histogram). Expression was normalized to the library size of each sample and the length of each gene.

The pattern of gene expression across the Microcystis phage Ma-LMM01 genome held fairly constant, with most genes following a similar dramatic increase in the third phosphate replicate. This includes genes that serve as typical markers for lytic activity, such as the tail sheath gene (gp_91), the major capsid gene (gp_86), and the terminase gene (gp_118) (Stough et al., 2017; Figure 3). A hypothesized auxiliary metabolic gene thought to be induced by phosphate starvation, phoH, also showed the dramatic increase in transcription in the third phosphate replicate relative to all others. In contrast, expression of genes associated with lysogeny, such as the integrase/recombinase (gp_136) and transposase (gp_135), displayed consistent and relatively low transcription across all bottles, including all of the phosphate bottles (Figure 4).

FIGURE 4

Average (large symbols) and individual bottle (small symbols) normalized expression of a Microcystis phage marker of lytic activity (gp91 tail sheath, black circles) and a marker of lysogenic activity (gp136 recombinase/integrase, open triangles). Expression was normalized to both library size of each sample and the length of each gene.

Overall, the BLASTx analysis indicated that transcription of non-Microcystis virus genes, either prokaryotic- and eukaryotic-infecting viruses, were nearly absent (data not shown). An exception included a Synechococcus-like myovirus terminase gene that was more transcriptionally active than Microcystis phage in the T0 and ammonia-amended bottles (Supplementary Figure 4).

Correlation Between Phage Gene Expression and Metabolites

A total of 58 identified metabolites were detected in the community, all of which were being actively cycled at the time of sampling (Figure 5A). Almost half (24 of 58) of the metabolites were positively correlated to expression of the Ma-LMM01 tail sheath gene (gp_91), an indicator of lytic infection activity [false discovery rate (FDR) corrected p-value < 0.05] (Figure 5B). One metabolite (orotate) was inversely correlated to gp_91 expression. Many of the significantly correlated metabolites have been linked to phage expression in studies of other systems, including UDP sugars and metabolites involved in the pentose phosphate pathway (Thompson et al., 2011; Ankrah et al., 2014).

FIGURE 5

Heatmap of metabolite peak intensities (A). Color scale has been scaled to each metabolite (z-score). Histograms describe the sum of normalized expression for Microcystis phage Ma-LMM01 in each experimental bottle. Expression was normalized to the library size of each sample and the length of each gene. Scatter plot shows the Pearson’s correlation coefficient between each metabolite detected and Microcystis phage Ma-LMM01 tail sheath normalized expression (B). Red dots and metabolite labels indicate a significant correlation (p < 0.05), while black dots and metabolite labels indicate correlations deemed to be non-significant statistically (p > 0.05).

It is unclear what proportion of the metabolite pools can be attributed to Microcystis activity. To gain insight into linkages between Microcystis and detected Microcystis phage, we correlated all Microcystis expression to the phage tail sheath (gp_91) expression across all of our samples. Supplementary Figure 5 shows the expression of each KEGG orthology group annotated as Microcystis, and many show the same pattern of increased expression in the third phosphate replicate. In fact, 497 of the 1,414 detected (35.1%) Microcystis spp. KEGG orthology functional groups were significantly correlated to the Microcystis Ma-LMM01 tail sheath gene (FDR p-value < 0.05). All of these correlations were positive, with none of the KEGG orthology groups that were negatively correlated (429) meeting statistical significance.

Discussion

Experiments in microbial ecology often use container-bound populations to explore responses to perturbation. These experiments involve replicate bottles or micro/mesocosms to increase confidence in the observations by reproducing results. When identical treatment bottles, or replicates, do not show identical responses, many scientists consider the discrepancies to be “experimental error” or “unexplained variability.” Variability can come in many forms: for example, in the current study it took up to 2 h to process samples after the 48-h incubation period, which provided for an extended (albeit minimal) incubation time for some samples. Unfortunately, this is a logistical limitation that likely effects all experiments with a well-replicated design. In the current experiments, however, at least some of this variability can potentially be attributed to subtle shifts in the microbiome of the target community, which can occur in either natural systems or in non-axenic lab cultures (Pound et al., 2021a). Moreover, as recently shown in other systems (Kuhlisch et al., 2021; Vincent et al., 2021) we can take advantage of the variability between replicates to learn more about our system.

Cyanobacterial blooms are often characterized by the amount of chlorophyll a present. While chlorophyll a concentration can provide a useful baseline for bloom biomass, it serves as a non-specific estimate of phototrophs in the community. It does not provide information on the diversity of phototrophic community nor the heterotrophic community. In this study, measurements of chlorophyll a response to nutrient additions revealed that the phototrophic community was constrained by available phosphorus. RNA sequence analyses confirmed that transcription in the phototrophic community was dominated by Microcystis, particularly in the phosphate-amended bottles. Phosphate constraints on Microcystis biomass has been well-documented in many natural bloom systems (Davis et al., 2010). Microcystin concentrations were not different between treatments (p-value > 0.1), nor were there differences in congener profile (Supplementary Figure 1). This is in many ways not surprising as it was recently demonstrated that up-regulation of the genes associated with toxin synthesis in Microcystis can take 48 h to result in a measurable increase in toxin quota (Martin et al., 2020). There was no significant change in chlorophyll a associated with nitrogen amendments, indicating that the phototroph biomass was not constrained by available nitrogen at the time of sampling.

However, as hypothesized nearly two decades ago (Wilhelm et al., 2003), nitrogen-amendments did cause a change in the community transcription. Our canonical correspondence analysis revealed that the microbial community transcriptional shifts occurred primarily within the heterotrophic bacterial community. This was particularly evident in the control and ammonia-amended bottles. The initial heterotrophic bacterial community was dominated by Bacteroidetes transcription, but after 48 h of incubation, the active heterotrophic microbial community in the control bottles had shifted, while the phototrophic community remained largely unchanged. The addition of ammonia appeared to maintain the transcriptional dominance of Bacteroidetes and Betaproteobacteria. This observation suggested that ammonia recycling was decoupled when the in situ heterotrophic community was constrained to a bottle. To phrase this another way, it suggests that at the time of sampling, the in situ bacterial community structure may have been dependent on ammonia as a main N source. Heterotrophic bacteria generally have faster growth rates than many planktonic phototrophs, which could allow for a more rapid change in that community because of changing environmental pressures (DeBruyn et al., 2004). This shift in community activity in the control bottles is one example of how bottle effects can alter community dynamics before consequences of experimental treatments are ever considered. We note that this response is not unexpected: recently we (Pound et al., 2021a) demonstrated that a single batch culture of non-axenic Microcystis aeruginosa undergoes dramatic shifts in the transcriptional activity of co-occurring heterotrophic bacteria after only three batch culture transfers in growth media containing different nitrogen sources.

The hypothesis that the in situ community was utilizing ammonia regeneration as a nitrogen source, even in high concentrations of nitrate (initial nitrate concentration was 30 μmol/L), finds support within the Microcystis functional gene expression data. The decrease in expression of nitrate/nitrite transporters in the ammonia- and urea-amended bottles suggested these communities were not relying on nitrate or nitrite as a primary nitrogen source. This is further supported by the increase in urea and ammonia transporters observed in the controls, indicating that Microcystis was preferentially using those nitrogen sources when enclosed in a bottle. When phosphate was added to some bottles, we saw a similar increase in urea and ammonia transporters, indicating a phosphate-constrained community that showed preferential uptake of urea and ammonia as a nitrogen source. This conclusion is supported by several recent studies. The utilization of regenerated ammonium has been shown to support non-N-fixing cyanobacterial blooms, such as Microcystis, in Lake Erie (Hampel et al., 2019), and a stable isotope study showed that Microcystis preferentially assimilates low concentrations of ammonium in the presence of high nitrate concentrations (Chaffin and Bridgeman, 2014). A preference for urea has been observed in many other Microcystis bloom events: urea is a common agricultural-based anthropogenic input, serving as both a nitrogen and carbon source to bloom organisms (Davis et al., 2010; Krausfeldt et al., 2019).

Characterization of the microbial community, based on rpoB transcriptional levels, revealed that bottle effects not only caused shifts between treatments, but can also cause shifts within treatments, resulting in variability between replicate bottles. Once a community is isolated in a container, it undergoes succession independent of other communities/containers regardless of treatment. Heterogeneity in natural microbial systems provides initial variation in bottled communities, and through independent succession, communities diverge further when contained. This phenomenon is apparent in several of the observations made from our experiments. On the simplest level, we recorded variable chlorophyll a concentrations between bottles of the same treatment. While such variation is common, the cause is rarely considered nor explained. It could be the result of a physical or chemical change, such as a shift in light levels, or a biological change in the non-target (i.e., in this case the heterotrophic microbial) community. It could either mean a shift in the type of phototroph present, or a shift in the photosynthetic activity of a single taxa, or some combination thereof. While there are a variety of potential sources of variability, shifts in the microbial community, either in function or taxonomy, could play a crucial role in observed variation. We also see bottle effects influencing intra-treatment variability in community composition. We documented variations in the heterotrophic community between bottles of the same treatment, where expression levels of the same group of organisms varied dramatically between bottles.

The most dramatic example of within-treatment variability was observed in transcriptional activity of Microcystis-infecting phage in the phosphate-amended bottles. The three replicate bottles showed different levels of phage gene expression, with one replicate displaying 10-fold higher transcriptional activity of Microcystis-infecting phage relative to another. These viral expression levels in the third phosphate replicate suggested a significant ongoing lytic infection during incubation. Consistent with this idea, Microcystis abundance (as measured by rpoB transcription) was ∼25% lower in this sample compared to the other P-amended replicates. In contrast, expression of the genetic marker for lysogeny was consistent and relatively low across all bottles, including the bottle with high transcriptional levels indicating lytic activity. This suggests that the possible lytic event was not the result of lysogenic induction (or at least that any lysogenic induction was incomplete). Such variable virus expression has recently been observed in other experimental systems, including an Emiliania huxleyi bloom induced in mesocosms, where researchers noted highly variable virus expression between replicate samples (Kuhlisch et al., 2021). In a likewise manner, their observations were tied to variability in host metrics, including chlorophyll a and metabolites.

The occurrence of a large-scale virus-mediated lytic event, even when a host population does not fully lyse, can have massive effects on host dynamics. We have shown that ∼35% of KO functional groups in Microcystis were significantly correlated to the phage transcript expression, indicating that the cells were responding to the viral activity. Many of these KO functional groups seemed to be associated with nucleotide recycling and urea utilization. Given that nucleotide recycling is a major component of virus replication, this seems like a rational relationship. We have also shown that ∼45% of the detected cellular metabolites from the entire community were significantly correlated with virus marker expression. The observed high levels of UDP sugars have been previously documented during viral infection of a marine heterotroph (Ankrah et al., 2014). Those authors suggested that these sugars may reflect changes in cell wall integrity as result of virus activity. Other metabolites, such as sedoheptulose bisphosphate and fructose-1,6 bisphosphate, are involved in the pentose phosphate pathways which has been hypothesized to provide extra energy and nucleotide precursors during viral infections of other cyanobacterial species (Thompson et al., 2011). Although we cannot determine the abundance of metabolites produced solely by Microcystis, it is possible that the community, including Microcystis, is responding to the viral-activity-generated “virocell” state (cf. Forterre, 2011) within this cyanobacterium: as Microcystis biochemistry changes due to infection the co-occurring community undergoes some parallel change. We note that not all viruses detected within these experiments demonstrated this variability between replicates. For example, myoviruses infecting Synechococcus provided very consistent, reproducible data across replicates, with elevated transcriptional activity observed in the ammonia-amended bottles that also had a higher Synechococcus host transcriptional activity (Supplementary Figure 4).

Without the analysis of viral transcription and metabolite pools, the presence of and the variation within the Microcystis viral infection could have been overlooked, with the resulting variability going unexplained. Indeed, any resulting phage influence on Microcystis transcription might be misconstrued, leading to invalid conclusions on the dynamics of the effect of phosphate. While not as dramatic in the current study, the same conclusions could be erroneously arrived at due to shifts in the function of the heterotrophic community. Functional relationships between Microcystis and its associated heterotrophic counterparts are still unclear, but many studies acknowledge a possible mutualistic relationship that can alter Microcystis growth dynamics (Kim et al., 2019; Cook et al., 2020). The “minor” community shifts that occur due to bottle effects have the potential to cause functional changes in Microcystis, introducing variability that might otherwise be ignored or left unexplained. Indeed, the apparent reliance on ammonium regeneration within the heterotrophic community could be a source of such variability within a 48-h timescale. To this end, researchers need to carefully consider the dynamics of any incubation experiments they undertake, as it is highly likely that this divergence due to “unintended biology” increases with length of time over any incubation period.

One surprising observation in this data set is the lack of transcripts assigned to non-Microcystis-infecting phage. Viruses are common components of metagenomic (i.e., DNA) sequencing efforts, with phage commonly having significant representation. This is expected given estimates of 105–107 phage per milliliter of water (Wilhelm and Matteson, 2008; Wigington et al., 2017). Yet in transcriptional analyses of several systems (e.g., Stough et al., 2018; Pound et al., 2020), we have observed that transcripts that can be assigned to phage appear more limited than expected. Given that our marker gene approach can find phage in DNA (genomic) samples, it is surprising that the same marker gene approach does not detect a higher expression of phage transcripts. In the current study, we again see phage contributing a surprisingly limited part of the transcriptome: while transcripts from phage infecting Microcystis are abundant, all other phage transcripts account for only a minimal component of library reads. While many scenarios might explain these observations, the most parsimonious would be that phage infection of the co-occurring community remains rather limited. While the idea of phage infecting the most abundant host is not new (Murray and Jackson, 1993; Thingstad and Lignell, 1997; Wilhelm et al., 1998), the extremely low abundance of phage transcripts would seem to fly in the face of past ideas on the importance of viruses in systems ecology (Wilhelm and Suttle, 1999) and perhaps requires further examination.

The complexity of biological systems cannot be overstated. Each organism (including viruses, heterotrophic bacteria, cyanobacteria, eukaryotic microalgae, and fungi) can play a pivotal role in cyanobacterial bloom dynamics. Going forward, this work demonstrates a need to broadly consider the influence of each entity in the ecosystem under study. Moreover, casting light on variability within basic experimental methods will not only provide a better understanding of the community of interest, but of biology (and all its possibilities) in general.

Data Availability Statement

The datasets presented in this study can be found in online repositories. Repository names and accession number(s) are provided in the article and the Supplementary Material.

Author Contributions

HP and SW designed the experiments. HP, BZ, RM, and JC completed the sample collection and the experimental incubations. HP and RM processed the samples. CC, SH, HC, and SC performed the mass spectrometry for the metabolomes. HP analyzed the transcriptomes and metabolomes. All authors participated in development of the manuscript and manuscript revision.

Conflict of Interest

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.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.