The experiments conducted here used a community of microbes sampled in the open water of the Long Island sound in front of the University of Connecticut – Avery Point Campus (41.30° N, 72.06° W) on 15 March 2016. We collected in situ water in four 20 L carboys, filtering first through an 80-μm mesh to remove predators and other large organisms, enabling us to focus on the protist community.

To perform our microcosm experiment, we used four cultures: three phytoplankton (prey) cultures Tetraselmis chui, Isochrysis galbana, and Phaeodactylum tricornutum; and a culture of Acartia tonsa (copepod, predators). The species of phytoplankton used are common in aquaculture application, represented various size, shape, pigment composition, and lineages, and have been used successfully in the past (McManus et al., 2004, 2012; Grattepanche et al., 2019). All phytoplankton were from the culture collection of the National Marine Fisheries Service Laboratory in Milford, CT, United States, and were grown in F/2 medium prior to experiments. The copepods were picked from cultures maintained at the University of Connecticut Department of Marine Sciences. In order to achieve the final concentration of phytoplankton for our microcosm experiment, we concentrated the combined phytoplankton by centrifugation at 3.9 × g to create a stock concentration of 5 × 10⁷ cells mL^–1 with the three algae within the same order of magnitude (P. tricornutum was added at 1/5 of the concentration of the two other species). We estimated the culture concentrations under a light microscope.

We used dialysis tubing to perform our microcosm experiment and incubated samples in a sea table with circulating water from collections sites to minimize impact of isolation, following the approach of Grattepanche et al. (2019). The dialysis tubing is composed of a cellulose membrane that is pervious to molecules <12,000 molecular weight (product D9402, Sigma), allowing for exchange of nutrients during incubation in sea tables. Each microcosm (tied-off dialysis tubing) contained 1 L and was floated in a sea table with continuous in situ seawater circulation as recommended by Capriulo (1982).

To test impact of varying levels of predators and food resources we first pooled our four in situ water samples together, and then divided them back into four carboys of 20 L to which we added: (1) nothing (control; p0), or phytoplankton at (2) 5 × 10² cells mL^–1 (p1), (3) 5 × 10³ cells mL^–1 (p2), and (4) 5 × 10⁴ cells mL^–1 (p3) in final concentration. For each of the carboys, we filled nine microcosms with 1 L of the seawater with (or without) phytoplankton added. For each of the phytoplankton abundance conditions, including controls, we added zero (z0), 5 (z1), or 10 copepods (z2) to triplicate microcosm for each treatment (three sets of dialysis tubing per phytoplankton/zooplankton combination). Based on insights from a preliminary study (Grattepanche et al., 2019), the microcosms were incubated for 3 days and then the total contents of each bag were collected for nucleic acids extractions. We used a 3-day incubation to minimize cascade effects and maximize nano- and microsized plankton impact based on insights from Grattepanche et al. (2019). We also collected duplicates of the initial in situ water sample and the in situ sample plus high phytoplankton abundance added (t0p3) at the beginning of the experiment.

Prior to isolation of nucleic acids, each sample was passed through serial filtration at 80 μm (to remove predators and larger organisms to focus on small protists and avoid some PCR inhibitors) and then collected on 10 (referred to throughout as the microsized fraction) and 2 μm (referred to throughout as the nanosized fraction) filters. We isolated organisms of the entire content of each dialysis bag on polycarbonate filters after size-fractionation, cut filters in two, and isolated DNA and RNA using appropriate kits. The filters do not exactly line up with current definitions of nano- and microplankton, as some have argued picoplankton are up to 3 μm in diameter (Vaulot et al., 2008) while nanoplankton can include species up to 20 μm (Sieburth et al., 1978). We use the terms “microsize” and “nanosize” as the bulk of the lineages caught on our filters likely fall in these sizes classes, the plasticity of species and the irregularity of some body plans all confound the efficiency of filtering.

For the nucleic acid extraction, we cut the resulting filters in half, storing one half for DNA (DNA prep buffer; 100 mM NaCl, Tris-EDTA at pH 8, and 0.5% of SDS) and the other half for RNA [RLT lysis buffer (Qiagen, Germany) plus 10 μL of beta mercaptoethanol per milliliter of RLT lysis buffer as recommended by the manufacturer]. DNA and RNA were extracted using Zymo Research soil extraction kit (Zymo Research, CA, United States) and Qiagen RNeasy (Qiagen, Germany) kits following the manufacturers’ instructions. For the RNA prep, residual DNA was removed using the Turbo DNA-free kit (Invitrogen, CA, United States). This approach (both filtration and extraction) was tested on previous samples and resulted in good DNA/RNA quality (Tucker et al., 2017; Sisson et al., 2018; Grattepanche et al., 2019).

The amplicon analyses followed protocols from Grattepanche et al. (2016b) and Sisson et al. (2018). In short, we used a primer set specific of the SSU-rDNA gene of the SAR lineages, which amplified a 150-bp fragment of the hypervariable region V3, to amplify DNAs extracted from filters. As discussed in Grattepanche et al. (2019), these primers amplify almost all known SAR lineages, excepted the highly divergent Foraminifera, which are also removed by our prescreening on 80 μm mesh). The amplification was done with the Q5 polymerase enzyme (NEB, MA, United States) following manufacture’s protocol. PCR products were cleaned using Agencourt AMPure XP beads (Beckman Coulter, CA, United States) and sent to University of Rhode Island Genomics and Sequencing Center for sequencing on an Illumina MiSeq sequencer (2 × 150 cycles). This resulted in a dataset of 7,207,909 reads. The raw reads are available from NCBI under the BioProject PRJNA550423 and the Sequence Read Archive SRP212194.

The Amplicon dataset was analyzed using a pipeline combining third party tools and custom python scripts (Sisson et al., 2018; Grattepanche et al., 2019). We first generated paired-end reads using Paired-End reAd mergeR (Zhang et al., 2013), refined the sequences (removed unpaired sequences; sequences without primers), created an OTU library using SWARM (v 2.1.9; Mahe et al., 2015), refined the OTUs (removed chimeric OTUs, outgroup OTUs, and OTUs with <10 reads), and then assigned taxonomy by tree using a curated SAR SSU-rDNA gene database¹. The final step corresponds to the rarefaction, i.e., subsampling a fixed number of reads for each sample to enable comparison. This resulted in a total of 1,148 OTUs corresponding to 3,400,000 rarefied reads (50,000 reads per sample).

To compare the effect of experimental conditions on diversity (i.e., taxonomy), only the OTUs present in two replicates of the same treatment and the same size fraction, and with >5 reads were kept. The reported number of reads for each treatment and size corresponds to the average of the read number among the three replicates, resulting in a total of 277 OTUs. We used this set of OTUs to assess changes between the control and the treatments.

To characterize expressed eukaryotic genes (i.e., those with polyA tails) we constructed mRNA libraries using TruSeq Stranded mRNA Sample Prep (Illumina, CA, United States) following manufacturer’s protocol. In summary, we isolated RNA using oligo-dT beads, fragmented the transcripts, and synthesized first and second strands of DNA (cDNA). The 3′-end of the fragments was adenylated to avoid ligation to another fragment, the adapters were ligated, and the fragments were amplified using specific PCR (see TruSeq protocol). DNA purifications were performed between each of these steps using Agencourt(r) AMPure XP beads (Beckman Coulter, CA, United States). The remainder of the protocol was performed at the University of Maryland, Baltimore Institute for Genome Sciences. This included library quantification and quality checking with an Agilent Technologies 2100 Bioanalyzer, pooling, and sequencing with an Illumina HiSeq 4000. We used a similar depth of sequencing for all our samples. While this depth of sequencing does not allow access to the lowly expressed transcripts, the use of similar depth of sequencing across samples allowed us to better compare these different samples. This resulted in a dataset of 783,472,526 reads, which are available from NCBI under the BioProject PRJNA550423 and the Sequence Read Archive SRP212194.

The metatranscriptome libraries were assembled using rnaSPAdes (version 3.10.1 with default parameters; Bankevich et al., 2012). We assembled the libraries in three ways: (1) each library independently (replicate set), (2) each triplicate by size fraction (nano- and microsized fractions, separately; treatment set); and (3) all libraries together in order to create a reference database for read number and RPKM calculation (Reads Per Kilobase Million = Reads Per Kilobase of transcript, per Million mapped reads, which is a proxy for gene expression).

As the reference genomes of microbial eukaryotes are still rare (del Campo et al., 2014), we looked at the pattern of conserved gene families (GFs) using PhyloToL (Cerón-Romero et al., 2019), which relies on the classification of GFs determined in OrthoMCL (Li et al., 2003; Chen et al., 2006). The resulting 1,485,323 transcripts were then refined using the “GF assessment for taxa” part of the PhyloToL pipeline (Cerón-Romero et al., 2018, 2019; Maurer-Alcala et al., 2018). In summary, we removed transcripts <200 bp (77,230 transcripts <200 bp), ribosomal DNA (5,699 transcripts matching rDNA gene references), and bacterial (38,305 transcripts matching reference bacterial genomes) contamination, and then binned the remaining sequences in orthologous groups (hereafter OG; proxy for GFs) using OrthoMCL (Li et al., 2003; Chen et al., 2006) as reference (990,632 transcripts did not match our reference GFs).

At this stage, the dataset was composed of 373,457 transcripts. We considered only the transcripts matching an OG with an E-value cutoff of 1e^–50 and a coverage of 10 (coverage as calculated by rnaSPAdes; 18,599 transcripts were discarded). This results in a reference dataset composed of 86,951 transcripts representing 5,258 conserved eukaryotic GFs (i.e., OGs from OrthoMCL; Chen et al., 2006) plus a large number of lineage-specific genes. We then decided to remove transcripts from the phytoplankton we added. For this purpose, we isolated transcripts from the added phytoplankton: we removed transcripts that were present in the in situ sample plus high phytoplankton abundance (t0p3) but absent from the in situ samples without added phytoplankton (t0). In this way, we identified 5,159 transcripts from the phytoplankton we added. We also removed the transcripts that were 10 times more expressed in t0p3 than in t0, leading to an additional 158 transcripts. Of the original reference transcripts (generated by combining mRNAs from all treatments), we removed 5,317 transcripts that represent genes expressed by the added phytoplankton, resulting in a refined reference dataset of 81,634 transcripts. We then used Blast2GO to assign gene ontology (GO; Conesa et al., 2005) and eggNOG (Huerta-Cepas et al., 2018) to assign clusters of orthologous groups (COG; Tatusov et al., 2000) and Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology (KO; Kanehisa et al., 2015).

To assess gene expression for the two sets (replicate and treatment sets), we mapped back the reads of each metatranscriptome to the refined reference (all libraries pooled together and refined using PhyloToL) using Seal (Sequence Expression AnaLyzer, version 35.92) implemented in BBmap (Bushnell, 2016). In short, each read was mapped to our reference using a Kmer of 31, and the count was transformed in RPKM to take into account the depth of sequencing for each sample. This software and a custom python script produced a table with the RPKM (read number or gene expression) of each reference transcript present in each of our samples (similar to an “OTU table”).

To assess the effect of each treatment on the community composition, we used the Unifrac dissimilarity index (Hamady et al., 2009; Lozupone et al., 2011) and principal coordinates analysis (PCoA). The analyses were performed in R using the Phyloseq (McMurdie and Holmes, 2013) and vegan packages (Oksanen et al., 2007, 2016). We tested the significance of the PCoA axis and factor (phytoplankton or predators) using envfit implemented in vegan, in which the data were randomly permuted 999 times (total of 1,000 datasets) and the random data were compared to the experimental data generated through a fitted regression model. Under the null hypothesis of no relationship between the ordination “axis” scores and the environmental variable, the observed R²-value should be a value among the permuted R²-values. However, if the observed R² is extreme compared to the permutation distribution of R² then we can reject the null hypothesis. The proportion of times a randomized R² from the distribution is equal to or greater than the observed R² is a value known as the permutation p-value. The same analysis was performed for the metatranscriptomics data using Canberra dissimilarity index (Lance and Williams, 1967).

Differential gene expression was calculated as (1) difference of read number between control and treatment or (2) ratio of RPKM between control and treatment. If the ratio was >1, the transcript was upregulated in the treatment compared to the control; if the ratio was <1, the transcript was downregulated, and we inverted the ratio (e.g., a ratio of 0.5 becomes −2). To be conservative, we only consider up/downregulation when the ratio was larger (or smaller) than 2 (or −2 for downregulation). In addition, we also assessed the OTUs and transcripts significantly up- or down-regulated (log2 fold change with p < 0.05) using DESeq2 package (Love et al., 2014). The transcripts significantly up- and downregulated assessed with DESeq2 matched with our estimates.

The diversity of SAR lineages remained relatively constant despite increasing predation pressure and food availability (Figure 1 and Supplementary Figure S1). Because the three replicates were similar overall (Supplementary Figure S1), we pooled replicates to evaluate the impact of treatment and considered only the OTUs (i.e., species) as described in the section “Materials and Methods.” We observed that the number of abundant OTUs (i.e., OTUs with more than 10 reads) varied between 63 and 125 (Figure 1A), and up to 239 OTUs when including the rare (Figure 2). The richness of abundant OTUs showed no clear pattern in response to phytoplankton abundance or copepod density (Supplementary Figure S1). In other words, the perturbations (adding prey and/or predators) did not have a strong impact on the overall diversity of SAR lineages.

While the overall diversity did not change within our microcosm, the community composition did respond to the various treatments (i.e., food availability and predation pressure). Overall, up to 10% of the nanosized OTUs significantly decreased in read number with increasing phytoplankton abundance, and up to 5% of the OTUs significantly increased in read number with increasing predation pressure (copepods, Figure 2). By comparison, the microsized plankton community showed a more neutral/stochastic response, with no clear relation with increasing phytoplankton and/or copepods (Figure 2). We also observed that 81% of the community variability is explained by differences between the nanosized and microsized plankton (R² = 0.7894, p < 0.0001; Supplementary Figure S2).

Among the nanosized community, the contribution of ciliates relative to stramenopiles increased with increasing phytoplankton (i.e., p0z0 vs. p1z0 vs. p2z0 vs. p3z0 showed an increasing contribution of ciliates; Figures 1A,B). The pattern among microsized species was more complex. For instance, ciliates, especially within Spirotrichea, increased with increasing phytoplankton in the absence of predators (Figure 1). However, this effect appears buffered by the increasing copepods as the ciliate contribution was almost constant when incubated with high predation pressure (i.e., 10 copepods L^–1) regardless of the concentration of phytoplankton (Figure 1). This suggests that while microsized ciliates responded to the increase of food, the copepod predation was able to keep them at a constant abundance. By comparison, we only observed a slight impact of copepods on nanociliates (Figure 1). In addition, we observed a slight decrease of phytoplankton (mainly Dictyochophyceae and Coscinodiscophyceae; Figure 1B) with increasing phytoplankton and a more mixed response with increasing copepods (Figure 1 and Supplementary Material).

Only a few OTUs, our proxy for species, responded repeatedly and similarly to all treatments. The control [no copepod and no phytoplankton added (p0z0)], the high predation pressure [high copepod abundance and no phytoplankton added (p0z2)], and the phytoplankton bloom incubations [high phytoplankton abundance and no copepods added (p3z0)] shared a quarter of the total OTUs (56 of the 212 OTUs observed in p0z0, in p3z0, and in p0z2; Supplementary Figure S3). The dominant OTUs did not show clear response to our treatments (see Supplementary Material; Supplementary Figure S2). A few OTUs increased or decreased significantly with food availability including OTU302 (closely related to the stramenopile Rhizochromulina) and OTU74 (likely an oligotrich ciliate) or with both food availability and predation pressure such as OTU94 (closely related to the euglyphid Cyphoderia major; Supplementary Figure S4). In the same way, only a third of the OTUs are specific to a size fraction and treatment (70 OTUs; Supplementary Figure S3), and these OTUs represented <1% of the community based on read number.

To assess changes in the overall community, we used PCoA with UniFrac index, which considers phylogenetic relationships among OTUs. Overall, these analyses indicate that predators impacted the microsized SAR community while the phytoplankton changed the nanosized SAR community. Adding copepods significantly impacted the microsized community (phytoplankton R² = 0.3689, p = 0.2238; copepods R² = 0.3726, p < 0.05 by envfit test, Figure 3), while phytoplankton addition had a significant impact on nanoplankton (phytoplankton R² = 0.6446, p < 0.05; copepods R² = 0.2976, p = 0.1638; Figure 3). Looking at lower taxonomic levels, we observed that Ciliophora, Stramenopila, and within Stramenopila, Dictyochophyceae in the nanosized fraction showed a significant response to phytoplankton treatments, while none of the microsized fraction showed significant response with increasing phytoplankton (Table 1). The copepods seemed to have only impacted the Bacillariophyceae from both nano- and microsized fractions (Table 1).

Our metatranscriptomics analysis focused on a dataset composed of 81,634 transcripts that match conserved eukaryotic GFs, and represent 13,029 GO terms, 5,176 conserved eukaryotic GFs, and 3,699 KO terms. We evaluated the impact of our treatment on community function using the number of conserved GFs and transcripts at three levels of expression: present, expressed, and highly expressed (>0, >10, and >1,000 as gene expression estimated by RPKM, respectively). The nano- and microsized fractions were composed of a similar number of expressed transcripts (40,667 and 41,369 transcripts, respectively) with half shared in both sizes (21,940 transcripts; see Supplementary Material).

Differential expression profiles showed a clear contrast between the two size fractions, and more intriguingly, a large number of transcripts were downregulated in the nanosized community incubated with phytoplankton at bloom conditions (Figure 4). The nanosized plankton showed up to seven times more downregulated transcripts when incubated in phytoplankton bloom conditions (Figure 4). In addition to transcript number, we also evaluated the impact of our treatments on the community function using three metrics: the number of conserved GFs, GO, and KEGG orthology (KO) terms associated with up- or down-regulated transcripts (Figure 2 and Supplementary Figure S5). Using all of these measures, we again saw a substantial change in nanosized plankton functions (up to 70% of transcripts and conserved GFs were downregulated in bloom conditions; Figures 2, 4 and Supplementary Figure S5), while the impact on microsized plankton was less marked (Figures 2, 4 and Supplementary Figure S5).

In comparison to the impact of prey abundance (i.e., phytoplankton), the impact of predator density on function was in agreement with our expectations: negligible for the nanosized fraction and leading to an increase in upregulated function in the microsized fraction (Figures 2, 4 and Supplementary Figure S5). Looking in detail at the interactions between phytoplankton bloom and high predation pressure, the increase in copepods density tended to increase the proportion of differentially expressed transcripts (Figures 2, 4 and Supplementary Figure S5). The number of upregulated functions in the microsized plankton increased with the increasing copepod abundance (Figure 2 and Supplementary Figure S5). For both sizes (nano- and microsize), the copepods tended to reduce the proportion of downregulated transcripts and GFs during the phytoplankton bloom incubations (Figures 2, 4 and Supplementary Figure S5). In other words, the impact of copepods seemed to counterbalance the effect of a phytoplankton bloom on the community function and vice versa.

Overall patterns of transcripts assessed by PCoA with a Canberra distance metric showed the same pattern as the SAR amplicon sequencing, i.e., samples grouped firstly by size fraction (nano- and microsized) and then by phytoplankton abundance for the nanosized community (35% of the variance expression; Figure 3). The expression within the microsized community was similarly impacted by phytoplankton and copepods (18 vs. 16% of the variance, Figure 3). This again suggested that while the nanosized plankton was primarily impacted by the change in phytoplankton abundance, the microsized plankton expression profile was a result of interactions between phytoplankton and copepod abundances.

To further assess community function, we assigned transcripts into functional categories using COG. One-third of the gene expression (as measured by RPKM) and one-fifth of the number of transcripts were involved in translation (category COG J). Another 40% of the transcript expression and transcript number was represented by transcripts from post-translational modification and protein turnover (COG O), unknown function (COG S), cytoskeleton (COG Z), energy production and conversion (COG C), and carbohydrate metabolism and transport (COG G). The unknown functions (COG S) included transcripts related to fucoxanthin chlorophyll a/c (55% of the reads for this study) and chlorophyll a/b binding proteins (6% of the read) by similarity.

To evaluate the main functions impacted during our experiments, we pooled the transcripts by COGs and estimated the differential expression in each treatment relative to the control. We observed a downregulation in cytoskeleton mRNA (COG Z in Figure 5) and a slight upregulation in mRNAs involved in translation in the nanosized communities incubated in phytoplankton bloom condition (COG J in Figure 5). The microsized community incubated with many copepods experienced an upregulation of transcripts of unknown function (mainly transcripts related to fucoxanthin chlorophyll a/c by similarity; COG S in Figure 5) and a downregulation of transcripts involved in translation (COG J in Figure 5). The other treatments showed only slight differences and no pattern was discernable between COG and incubation condition.

Phylogenomic analyses indicated that transcripts strongly impacted in our microcosm were from SAR lineages and other phytoplankton (e.g., chlorophyte, cryptophyte, and haptophyte). We generated phylogenies for 11 exemplar genes that showed high diversity (within the top 10 genes for the number of transcripts) and/or high expression (among the top 10 gene for the RPKM), and were annotated in reference databases (e.g., three ribosomal proteins, α and β tubulins, actin, HSP70, chlorophyll a/b binding proteins, and glyceraldehyde-3-phosphate dehydrogenase; Supplementary Figure S6). For each of these genes, SAR represented almost 70% of the transcripts; Stramenopila, Alveolata, and Rhizaria representing approximately 45, 40, and 10% of the highly expressed transcripts, respectively. Other algae (e.g., chlorophytes, haptophytes, and cryptophytes) comprised ∼25% of the highly expressed transcripts, while Opisthokonta, Amoebozoa, and Excavata were represented by <6% of the highly expressed genes (Supplementary Figure S6).

While ciliates clearly responded to increased phytoplankton abundance in amplicon analyses (Figure 1), our taxonomic analyses of some GFs did not find evidence of their impact on the community function. Ciliates contributed up to 30% to some genes including eEF1α (Supplementary Figure S6). However, ciliates did not show any strong response to treatments for this set of genes except a more important contribution to the actin expression in our phytoplankton bloom treatment (Supplementary Figure S6). Other cytoskeleton proteins, such as α-tubulin and β-tubulin, did not show this pattern. We looked at other genes involved in phagosome formation such as ARP complex and RAC as a proxy of phagotrophy (Yutin et al., 2009) and the proteins identified by Burns et al. (2018), but again we did not observe a clear pattern across treatments (Supplementary Figure S7). While this was surprising, it is important to remember that many transcripts (almost 1 million) were not considered here as they were not annotated in databases and instead might represent lineage-specific genes. In other words, the change in function of ciliates may be through lineage-specific genes that lack annotated homologs in current databases.

The nanosized plankton (i.e., 2–10 μm) response to increasing phytoplankton (three species ranging in size from 4 to 15 μm) was unexpected. We added small sized phytoplankton expecting to see an impact on their predators within the microsized (i.e., 10–80 μm) heterotrophs and saw a small increase of microsized ciliates with increasing phytoplankton (Figure 1). In other in situ studies, microheterotroph biomass increased with phytoplankton and nanoplankton were responding to bacteria (Irigoien et al., 2005; Grattepanche et al., 2011a). So, we did not expect to observe an increase of nanosized lineages, assessed by SAR amplicon analyses, incubated in bloom conditions (Figure 1). This suggests that either nanosized lineages (i.e., ciliates) were able to consume prey of similar to their own size or that their response is indirect. For example, it is possible that phytoplankton exudate caused increased growth of bacteria, which are in turn grazed by nanosized ciliates (i.e., more food leading to increased abundance) (Fenchel, 1987). Another possibility is that predators of nanosized ciliate fed selectively on the same size phytoplankton, and therefore relieved predation pressure on these small ciliates (i.e., fewer predators leading to increased abundance).

We evaluate the possibility of lineage-specific changes in predation among nanosized community members by looking at genes involved in phagocytosis. While we observed an upregulation of actin for the nanosized ciliates in bloom condition (Supplementary Figure S6), we did not see an increase in expression for phagotrophic genes such as the ARP complex and RAC GF identified by Yutin et al. (2009) and a longer list of genes identified in two other studies (Burns et al., 2018; McKie-Krisberg et al., 2018; Supplementary Figure S7). These analyses are consistent with an indirect response of nanosized ciliates rather than an increase in grazing in response to added phytoplankton. However, this conclusion needs to be taken with caution as the mechanism of phagotrophy for microbial eukaryotes is still poorly understood (Yutin et al., 2009) and many lineage-specific GFs are likely missing from our database (Cerón-Romero et al., 2019).

The transcript number and expression levels for communities incubated in bloom conditions strongly decreased for the nanosized fraction (Figure 4 and Supplementary Figure S5). A decrease of functional diversity is expected when just a few species dominate the ecosystem (i.e., in bloom conditions), but the observed impact on nanosized and not on the microsized plankton is surprising. One possibility is that changes in ecological niche (abiotic and biotic factors) in the phytoplankton bloom conditions might have selected for only a subset of the nanosized species present at the beginning of incubation, resulting in a decrease of expressed function. A decrease of diversity and function has been documented during bloom events for bacterioplankton (Teeling et al., 2012; Wemheuer et al., 2014), but not yet for eukaryotes.

We did see patterns in changes of genes related to photosynthesis and stress. For example, we saw a reduction in expression of transcripts related to photosynthesis (COG S, Figure 5; chlorophyll a/b binding protein; Supplementary Figure S6). In addition, stress-related proteins (e.g., HPS90, HSP70), particularly of phytoplankton groups (e.g., Archaeplastida, Supplementary Figure S6), increased with increased phytoplankton. This suggests a stress among phytoplankton in response to bloom conditions. This is consistent with the observed decrease of SAR phytoplankton (mainly Dictyochophyceae and Coscinodiscophyceae) observed with the increase abundance of added phytoplankton (Figure 1).

The addition of varying levels of phytoplankton impacted the diversity and function of both the nanoplankton and, to a lesser extent, the microsized plankton. As expected for bloom conditions, we observed an amplification of changes in both up- and downregulated transcripts with increasing phytoplankton (Figure 2). However, none of these genes are directly related to phytoplankton blooms based on findings in previous studies (e.g., proteins involved in nutrient metabolism or carbohydrate metabolism; Zhang et al., 2019). We hypothesize that the change in community function is related to changes in metabolism during the extreme conditions generated by the bloom.

The species composition (i.e., OTUs) of the microsized community was impacted by top–down control, but not the functions of the microsized community. We had expected microsized ciliates to increase with phytoplankton concentration as we selected three phytoplankton species (T. chui, I. galbana, and P. tricornutum) considered to be “good” food resources for ciliates (Verity and Villareal, 1986; Stoecker et al., 1988; Christaki and Van Wambeke, 1995; McManus et al., 2012; Schoener and McManus, 2012). Instead, we only observed top–down control of ciliates by copepods, as ciliates represent a superior food resource for copepods and may be selectively grazed (Calbet and Saiz, 2005). Our data are a direct illustration that copepod grazing on microheterotrophs (i.e., ciliates) reduced top–down controls, allowing phytoplankton to increase (Irigoien et al., 2005; Behrenfeld, 2010; Kuhn et al., 2015). While this has been hypothesized from microscopic observations of in situ samples (Irigoien et al., 2005; Leising et al., 2005), here we showed, in our closed system, that ciliates were not able to control the amplitude of a phytoplankton bloom because of the predation pressure from zooplankton (copepods).

Copepods did have a significant impact on microsized phytoplankton community function as photosynthetic activities (e.g., chlorophyll a/b binding proteins) increased with increasing copepod abundance (Figure 5 and Supplementary Figure S6). Copepods have been reported to release dissolved organic carbon and inorganic nutrients through sloppy feeding and fecal pellets (Saba et al., 2011), which may increase phytoplankton production. Together these factors (increase in nutrients and release of predation pressure by feeding on microciliates) may explain the functional changes in microsized community members observed across our experiments.