Monthly CTD profiles of temperature, salinity, density, oxygen, and chlorophyll fluorescence paired with discrete bottle sampling, provided hydrographic context. The details of sampling, quality control and post-processing are available on the BATS website (http://bats.bios.edu/index.html). Each profile was assigned to one of four seasons (mixed, spring, stratified, fall; Table 1 ) based on retrospective evaluation of the full time series for each year. Hydrographic layer boundaries were determined based on a set of pre-defined physical criteria ( Table 1 ) on a profile-by-profile basis from the fluorometer and density criteria.

Biological samples were acquired monthly from February 2016 to December 2018. March 2016 was sampled twice (March 8th and 25th); March 2017 sampling was cancelled due to rough weather. The protocol included 12 fixed depths between surface and 1000 m (0 m, 40 m, 80 m, 120 m, 160 m, 200 m, 250 m, 300 m, 500 m, 600 m, 800 m, 1000 m). At each depth, 4 L of water were filtered using 0.22 µM Sterivex filters (MERCK, MA, USA). Filtering pressure was kept at 5 mmHg. Both ends were sealed with parafilm or a sterile plug/cap. Filters were wrapped in foil and stored at -80°C.

Sucrose lysis buffer (1 mL) was added to the Sterivex filters upon thawing. DNA was extracted using the method by Giovannoni et al. (1990). Amplification and sequencing was carried out in the Center for Genome Research and Biocomputing at Oregon State University. The 18S rRNA V4 hypervariable region was amplified using the primers 5-CCAGCA[GC]C[CT]GCGGTAATTCC-3 and 5-ACTTTCGTTCTTGAT[CT][AG]A-3 (Stoeck et al., 2010), using the KAPA HiFi HotStart ReadyMix (Kapa Biosystems; Wilmington, MA, USA). The PCR thermal cycler program consisted of a 98°C denaturation step for 30 s, followed by 10 cycles of 10 s at 98°C, 30 s at 53°C and 30 s at 72°C, and 15 cycles of 10 s at 98°C, 30 s at 48°C and 30 s at 72°C, and a final elongation step at 72°C for 10 min. After PCR, dual indices and Illumina sequencing adapters were attached using the Illumina Nextera XT Index Kit. Amplicon sequencing was conducted on an Illumina MiSeq using 2x250 PE V2 chemistry.

Quality control and read merging were carried out with MeFiT Ver. 1 (Parikh et al., 2016), using Jellyfish Ver. 2 (Marçais and Kingsford, 2011), with strict parameters to compensate for the lack of complete overlap between forward and reverse: CASPER Ver. 0.8.5 (Kwon et al., 2014) K-mer size of 5, threshold for difference of quality score 5, threshold for mismatching ratio 0.1, and minimum length of overlap of 15. The Merge and Quality Filter Parameters (using the meep score (Koparde et al., 2017)) were set at a filtering threshold of 0.1, discarding non-overlapping reads. Assembled reads passing the QC were fed into MOTHUR ver. 1.43 (Schloss et al., 2009). In short, after trimming to the V4 region (by aligning to the SILVA 138 release trimmed to the V4 region) and removal of chimeras with VSEARCH (Rognes et al., 2016), single variants (=100% OTU; Porter and Hajibabaei, 2018) were obtained using UNOISE2 (Edgar, 2016) as implemented in MOTHUR (diffs=1). Taxonomic assignment was done using the PR² database (Guillou et al., 2013) using the RDP classifier in MOTHUR. After excluding metazoans and “Eukaryota_unclassified”, samples were rarefied to 20,000 reads and all ASVs with less than two reads (global singletons) were equaled to zero. Scripts are available at the GitHub repository https://github.com/blancobercial/Protist_Time_Series.

Counts per sample and environmental metadata were loaded into PRIMER Ver. 7 (Clarke and Gorley, 2015). Diversity indices (total number of taxa S, Shannon diversity index H’, and Pielou’s evenness index J’) were calculated in PRIMER; Faith’s phylogenetic diversity PD (Faith, 1992) was calculated in MOTHUR using the same input table as for the other indices, after calculating the tree using ClearCut in MOTHUR from the aligned dataset. Samples were then standardize by the total, and square-root transformed (Hellinger transformation) (Legendre and Gallagher, 2001). A Bray-Curtis distance similarity matrix was built, and a Principal Coordinates Analysis (PCoA; Gower, 1966) was carried out. The relationship between environmental variables and the similarity between samples was analyzed with BioEnv in PRIMER, using Spearman Rank and up to five variables for the model. The variables included day of year, depth, temperature, salinity, oxygen, fluorescence, turbidity and density. A second analysis included the factor “hydrographic layer” ( Table 1 ), to assess how this partitioning compared to the protist community. All environmental variables were normalized (subtracting the mean and dividing by the standard deviation) before analyses to remove biases associated to different variable scales.

Non-hierarchical k-means clustering based on group-average ranks was used to find the best grouping reflecting the self-organization of the community, from number of clusters K=2 until two of the clusters were not statistically different (K=14). At each K level, an analysis of similarities (ANOSIM) was run (with 999 permutations) to assess the R-value at each clustering level and if the clustering level was statistically significant (both globally and between pairs of clusters). The K levels showing the maximum R-values were further analyzed. Higher R values (closer to 1) represent higher differences between clusters. Clustering was mapped to the hydrographic layers to assess the fidelity of each community cluster to the pre-defined hydrography-based layering, evaluated within the PCoA ordination, and their taxonomic composition at PR² “Super-Group” (=Phylum) level analyzed. A test of homogeneity of multivariate dispersions (PERMDISP) was run to compare the dispersion among clusters, as a measure of intra-cluster variability.

The molecular lineages (single variants) with the largest contribution to the similarity in each cluster, and to dissimilarity between adjacent clusters, were assessed using a Similarity Percentage (SIMPER) analyses. Fine taxonomic assignment was done using BLAST (Altschul et al., 1990), analyzing several of the top scores to understand the precision of assignments (species/genus/family/etc.). After selecting the highest scores (and only those with similarity >90%), the lowest shared taxonomy assignment was retained. Interpretation of the seasonal and vertical changes was contextualize within the oceanographic landscape and transitions between clusters, and with the functional profile (based in the literature) of each lineage.

PCoA analyses were further performed at each discrete depth to study the interaction between seasonality and depth. Same procedure was repeated following those samples that overlapped with the Chl-a maximum peak for each profile (which did not happen every month because a set of fixed depths were sampled each time). SIMPER analyses were run to detect taxa related to seasonal succession at each depth.

For all analyses, a Holm-Bonferroni correction was applied to the significance threshold when multiple comparisons were done (Holm, 1979).

Over the 3-year period, hydrographic sampling revealed the expected seasonal cycles of cold season mixing, warm season stratification, the short spring transition, and the vertical distribution of chlorophyll fluorescence ( Figure 1 ). The hydrographic layer boundaries fluctuated in depth over the observational period, highlighting the ability of this physical framework to transcend the temporal variations in vertical water mass structure associated with winter mixing, mesoscale eddies, location of the nutricline, and varying light penetration ( Figure 1 , Supp. Figure S1 ). The depth of the surface mixed layer (Layer 0) varied seasonally from 10 m in the stratified season to 170-212 m in the mixed. The depth of the DCM (Deep Chlorophyll Maximum; Layer 2; Figure 1C ) ranged from 70 to 130 m. The top of Layer 8, corresponding to the deep O₂ minimum ( Figure 1B ), varied between 600 and 850 m. Other hydrographic features detected were the uplifts in the deep mesopelagic layers due to the passage of large eddies ( Figure 1 ). The most prominent uplifts (of about 200 m) occurred during 2018.

Of the 408 samples (34 sampling months, 12 depths), 369 produced a library. The remaining samples failed at some point of the procedure (sampling at sea or extraction), and did not produce a useable library. It usually affected several samples from the same cast, especially the deeper layers (where DNA concentration was always lower). About ~39 million reads were generated (average ~105,000 reads per sample; 23,000 to 364,000). Of those, the strict QC retained 42% as (16.5M reads; 44,000 per sample average) “very high quality” (MEEP definition for those passing the QC settings). There was no link between reads passing QC and months or depths.

Non-hierarchical clustering based on Bray-Curtis distances showed two local R maxima (that is, a maximum in the difference among clusters), for the number of clusters (K) K=2 (ANOSIM R=0.903; p<0.001) and K=9 (R=0.904; p<0.001). These clusters represent groups of samples with similar community composition that are statistically different from the communities in the other clusters. In both cases, all pairwise comparisons between clusters were significant. The K=2 partitioning roughly divided the samples between the epipelagic (photic) and the mesopelagic zones, although when the MLD (Mixed Layer Depth) reached into in the mesopelagic, the upper community followed the MLD depth ( Figure 2A ). The K=9 clustering aligned closely with the hydrographic layering, with higher numbers of different clusters in the epipelagic, and greater homogeneity in mesopelagic ( Figure 2B ). In some cases, single clusters spanned several hydrographic layers in the mesopelagic. The clusters precisely traced hydrographic events such as the deepening of the MLD in winter, the abrupt succession to a different community following the shoaling of this layer at the spring transition, and the uplift of mesopelagic layers due to the passage of mesoscale eddies ( Figure 2C ).

In the K=9 grouping ( Figure 2B ), the epipelagic zone was characterized by a succession of different communities over the annual cycle. Cluster 1 dominated the hydrographic Layer 0 (surface MLD) in the spring and stratified seasons, but was displaced by Cluster 2 during the fall and mixed seasons. Cluster 1 was detected as soon as the MLD shoaled in each year, and remained detectable until the onset of fall when the MLD (Layer 0) gradually deepened and Cluster 2 became dominant. Cluster 3, found below clusters 1 and 2, corresponded to the lower portion of the epipelagic zone that included the DCM ( Figures 2A, C ). Clusters 4, 5 and 6 occupied the underlying water mass (upper mesopelagic), corresponding to the local deep winter mixed layer (the Winter Mode Water, WMW; Layers 3 & 4), a very weakly stratified portion of the water column occupying the approximate depth range 150-400 m, that did not strictly align with density. Clusters 7, 8 and 9, in the lower mesopelagic (500 – 1000 m), aligned with the density fields. Cluster 8 overlapped with Layer 8, the oxygen minimum zone (OMZ, with O₂ < 160 µmol kg_-1). The mesopelagic layers were uplifted during the passage of the fronts in 2018, detected by the raised density layers ( Figure 2 ). Clusters were more tightly grouped with depth (PERMDISP test for homogeneity of multivariate dispersions; p < 0.0001), indicating a more variable community in the epipelagic compared to the mesopelagic, and in the upper mesopelagic compared to the deep mesopelagic.

The PCoA ordination arranged the samples with the first axis roughly corresponding to a combination of depth and environmental factors ( Figure 3 ), with many samples departing from their theoretical position if based only on depth (also shown in Figure 2B ). The BioEnv procedure determined that the best model to relate the sample ordination with the environmental variables combined O₂, fluorescence, depth and density (ρ=0.729 p < 0.001), but since depth and density are highly correlated, the next best model included O₂, fluorescence and density (ρ=0.690, p < 0.001). The proposed hydrographic layering was the best-correlated single environmental variable (ρ=0.651, p < 0.001). Superimposing the non-hierarchical clustering onto the PCoA, the K=2 showed a sharp divide between the clusters, despite the intrusions of the upper cluster into deeper layers described before ( Figure 1 ). For K=9, separation of clusters 1 & 2 reflected the seasonal shifts in the upper epipelagic layers (mixed vs stratified conditions). The remaining clusters projected consecutively in the PCoA, with very little mixing at the boundaries. The PCoA also indicated a larger separation between cluster 7 (occupying depths 500 – 600 m) and the other communities. In contrast, clusters 8 (OMZ) and 9 were very close to each other, with minimal mixing.

The taxonomic composition showed 40-45% of the reads corresponding to Syndiniales ( Figure 4 ), a parasitic group (Guillou et al., 2008); the proportion originating from free dispersal states versus the parasitic state is unknown. At K=2, the taxonomic composition showed a shift from a more diverse community in the epipelagic, with many different clades including autotrophs, heterotrophs and mixotrophs, to a heterotrophic, Rhizaria-dominated, community at depth ( Figure 5 ). This result was more pronounced if the Syndiniales were not considered in the analyses. The variables/single variants with the largest contribution to the similarity within each cluster, identified with the SIMPER analyses, indicated that the epipelagic cluster was characterized by a large presence of autotrophs (mostly Pelagophyceae Pelagomonas calceolata) and mixotrophs (from Ochrophyta, Stramenopiles and Dinophyta) although some of the latter clades include all trophic possibilities. In contrast, heterotrophs overwhelmingly dominated the deep cluster, principally Radiolarians (Polycystinea and Acantharea), and representatives of Stramenopiles (Labyrinthulea). The dissimilarity between groups was then driven by the mixture of autotrophs (especially P. calceolata and the Chlorophyta Ostreococcus sp.) and heterotrophs/mixotrophs (Stramenopiles MAST-4A, and several Dinophyceae lineages) in the epipelagic, compared to their virtual absence in the mesopelagic, where depth-specific heterotrophic Radiolaria dominated.

The taxonomic composition at K=9 reflected a more detailed partitioning ( Figure 4 ). Alveolata dominated the near-surface community during the spring and stratified seasons even if excluding Syndiniales (cluster 1), together with a mixture of Hacrobia and Stramenopiles. In contrast, the near-surface taxa during the mixed periods (cluster 2) exhibited greater abundances of Stramenopiles and Rhizaria. Cluster 3, below clusters 1 & 2, contained the highest concentrations of Archaeplastida, although Rhizaria was the dominant free-living group. A gradual decrease of all the non-Rhizaria groups was associated with increasing depth: Hacrobia and Stramenopiles disappeared almost completely by cluster 7, with only Alveolata maintaining a significant presence. Syndiniales abundances slightly decreased with depth, but had a greater presence in Cluster 7 (500-600 m). SIMPER analyses complimented these findings, determining which taxa showed a clear affinity with cluster. There was no single taxon with a significant presence in all clusters, and most were significant only in one or two clusters. A few Rhizaria, however, showed high numbers in several mesopelagic layers ( Figure 5 ; Table 1 ). Cluster 1 (surface stratified) was characterized by mixotrophs and heterotrophs of Alveolata Dinophyta (e.g., Warnowia sp., Karlodinium sp. or Lepidodinium spp.), Hacrobia Prymnesiophyceae (Chrysochromulina sp.) and several Stramenopiles. The only autotroph was a clade of Phaeocystis sp., although this group is known to appear as free living or as a symbiont autotroph in Rhizaria colonies. Cluster 2 (near-surface, mixed periods) shows large abundances of P. calceolata (autotroph), but most of the other clades belong to mixotrophs (e.g. several Dinophyceae, and MOCH-2), or heterotrophs (e.g., MAST lineages 25, 4A, 4C and 1D, Hacrobia Pterocystida). Cluster 3 (DCM) showed several mixotrophs among the dominant lineages, but the main groups were autotrophs (the small P. calceolata, Ostreococcus sp., Bathycoccus prasinos, and the Haptophyta Phaeocystis globosa). Heterotrophs characteristic from this clade included the Hacrobia Leucocryptos marina and several Rhizaria (although these might have autotrophic symbionts such as Phaeocytis). Clade 5, below the Chlorophyll maximum, showed significant abundances of the autotroph P. calceolata, however the shift towards a heterotrophic Rhizaria-dominated community was noticeable, together with heterotrophic (Telonemia sp., Leucocryptos marina) and likely mixotrophic (Prymnesiophyceae Clade E) Hacrobia. Different Rhizaria dominated all remaining clusters, while some heterotrophic Hacrobia and Stramenopiles were still among the dominant clades in clusters 5 and 6.

The dissimilarity between consecutive clusters reflected transitions in community function. The mixotroph-dominated cluster one (surface stratified) was replaced by the autotroph-dominated community occupying the ML during the periods of deep mixing (cluster 2). In cluster 3 there was a larger increase of the autotrophs at the cost of mixotrophs from cluster 2. Deeper, differences were due to the progressive decrease in autotrophs, increase of heterotrophs, and layer to layer replacements between different heterotrophs (mostly Rhizaria).

When considering all samples together, there was no discernible seasonal signal. Seasonality, however, appeared when analyzing single depths/clusters. Seasonality was evident in the top layers (surface and 40 m especially; ANOSIM R=0.66 and R=0.62 respectively; p<0.001), but the signal faded with depth, with ANOSIM R progressively decaying (indicating lowering differences among seasons), and becoming non-significant below 300 m. There was no significant seasonality in the mesopelagic. In the lower mesopelagic, the passage of eddies was, however, detectable as a group of dissimilar samples ( Figure 6 ).

Near the surface, the highly diverse community was composed of many taxa with low relative abundances (mostly < 2% average), dominated by mixotrophs and heterotroph clades. Only mixed and spring seasons had an autotroph (P. calceolata) among the most abundant clades (only representing ~ 2% of the reads), while a single heterotroph identified as Warnowia sp. represented 8% of reads in spring. In the stratified season, mixotrophs, and a diverse community of heterotrophs, including Telonemia, Warnowia, several Hacrobia Centroheliozoa and the Rhizaria Minorisa minuta dominated the community. In the fall, mixotrophs gained more prominence compared to heterotrophs.

Cluster 3 was generally associated with the Layer 2 (broadly defined DCM) but, since the DCM varied substantially, sometimes this layer did not capture the feature. Cluster 3 community exhibited depth dependence, but no apparent seasonal cycle. When DCM samples were restricted to those acquired within 90% of the actual chlorophyll maximum, a strong and significant seasonal pattern was detected (ANOSIM R = 0.762, p < 0.001; Figure 6 ). The stratified and fall periods, despite mixing in the 2PCos representation, were statistically different (R = 0.5, p < 0.05) and were separated in the 3-PCoA ( Figure S2 ). The SIMPER analyses performed to find the main contributors to the dissimilarity between these seasons in this narrow Chl-a max, showed that the main difference happened during the spring season, when a few autotrophic clades dominate the community (Ostreococcus sp., P. calceolata and B. prasinos). Three different clades of Micromonas, and P. globosa, were also among the most abundant clades. Together, autotrophs represented over 30% of the reads of the non-parasitic community. Other main groups included Dinophyceae (likely mixotrophs) and heterotrophs such as Telonemia sp. and Radiolaria (although these might have autotrophic endosymbionts). After transitioning to the stratified period, only P. calceolata, B. prasinos and P. globosa showed relatively high abundances in the DCM; however, their prevalence was much lower (just above 5%). In contrast, mixotrophs and heterotrophs increased in relative abundance. In the fall transition, P. calceolata and P. globosa were the only autotrophs among the main taxa (with slightly higher relative abundances compared to the stratified period), while the proportion of heterotrophs (Rhizaria and Telonemia) increased. During the mixing season, there was again a strong shift to autotrophs (especially P. calceolata, B. prasinus and Ostreococcus. sp.; about 10% of the total reads combined) and mixotrophs, with no free heterotrophs among the main lineages in the community.

Diversity indices indicated a mismatch between species-based and phylogenetic-based diversity indices ( Figure 7 ). The data showed a higher number of species (S), species diversity (H’) and evenness (J’) in the upper layers of the water column, gradually decreasing with depth, reaching a minimum in the cluster corresponding to the oxygen minimum zone. These profiles also revealed an interesting feature: superimposed upon the general trend of decreasing diversity and evenness with increasing depth throughout the upper 1000 m, the indices identified three subgroups (clusters 1-3, 4-6 and 7-9) each of which exhibited its own decreasing trend. The subgroups correspond respectively to the epipelagic, upper mesopelagic, and lower mesopelagic portions of the water column. The profile for Faith’s phylogenetic diversity was distinctly different from these: it showed rising diversity from clusters 1 to 4, followed by a steady decrease with depth to a minimum in the OMZ.

Our findings, showing higher diversity in the epipelagic, peaking below the Chl-a max, and then decreasing with depth, are opposite in pattern to the study by Ollison et al. (2021). These differences can be attributed to many causes. There is great disparity in the environmental conditions, with strong seasonality and fully oxic OMZ at BATS compared to the quasi-permanent conditions and the thick and suboxic OMZ at ALOHA (Church et al., 2013). Additionally there are pure methodological differences: the present manuscript is based on DNA metabarcoding, while Ollison et al. (Ollison et al., 2021) used an RNA-based metabarcoding approach. Our data, however, coincides with the findings of Canals et al. (2020), where a similar pattern of decreasing diversity with depth and a peak around the Chl-a maximum, was found for ciliates in the Atlantic Ocean, suggesting that the environmental landscape (the hydrography) is driving the differences between ocean basins.

The clear signal of increasing phylogenetic diversity towards the Chl-a-max, with the layers below the 1% (but above the mesopelagic) showing the highest diversity, is notable. This maxima might be related to the presence of more, and more diverse, ecological niches, since here light and nutrients are able to sustain a rich autotroph community (below it is too dark; nutrients are scarce above) in addition to the heterotroph and mixotrophs, facilitating the coexistence of more phylogenetically apart lineages. This pattern contrasts with the taxonomic diversity, either H’, or SVs richness, which did not show a noticeable change in the epipelagic. This mismatch is significant since PD is often directly correlated to species richness (Barker, 2008; Voskamp et al., 2020). The pattern here likely reflects the abrupt transition from an Alveolata-dominate surface community towards one more phylogenetically diverse, despite a relatively constant species richness, owing to the availability of nutrients.

The mesopelagic communities exhibited lower diversity and tighter grouping in the PCoA, indicative of less month-to-month variability compared to epipelagic communities. These patterns may be related to seasonality in the epipelagic layers, where a more dynamic environment would favor community succession processes month to month (represented as a higher dispersal between samples) and inhibits efficient competitive exclusion processes (causing higher diversity within sample). The deeper layers, however, characterized by less environmental variability, would favor the same lineages with competitive advantages month after month, limiting the number of taxa, and resulting in more stable communities. The vertical stratification of the Rhizaria community could be due to resource availability (decreasing with depth; Pedrosa-Pàmies et al., 2018), resource structure (e.g. presence and size of aggregates, filaments, or dissolved) and/or elemental stoichiometry with depth.

Metabarcoding approaches, as any other method, are subject to many methodological limitations and bias (Bucklin et al., 2016; Santoferrara, 2019). The initial selection of which region and primers to use, the marker copy number variability between taxa, etc., are a source of biases that would affect the description of the community. Other source of uncertainty is the taxonomic assignment of the different lineages (SV or OTUs), due to both incomplete databases (where there is a lot of work to do by the community), and to the lack of resolution at the species level by most ribosomal markers. To mitigate these biases, the taxonomy for the discussed lineages was analyzed against GenBank, to understand the ambiguity of the taxonomic assignment, and corrected if needed (compared to the mothur PR² based assignment). Regarding the copy number bias, we are confident that the core of this manuscript, the seasonal and depth community transitions, are not affected by this fact, since there are strong changes in the community composition and not small changes of the same taxa, and would very likely hold true independently of the marker used. Finally, the raw value of the diversity indices should be considered only in the context of this study, since methodological questions as the indicated or the rarefaction (or extrapolation) level chosen can drastically affect diversity estimates (Santoferrara, 2019; Blanco-Bercial, 2020). Despite this possible bias, it is very likely that the observed patterns of variability (where diversity is higher/lower, etc.) are however robust against those bias. Any comparison with other existing data would require a reanalysis with shared bioinformatics approaches.

The protist community in the Sargasso Sea shows a dynamic seasonality in the epipelagic, responding to hydrographic yearly cycles. Mixotrophic lineages, able to take advantage of the smaller picophytoplankton and heterotroph bacterioplankton, dominate throughout the year; however, autotrophs bloom during the rapid transition between the winter mixing and the stratified summer. Pure heterotrophs have their peak moment at the end of summer, when the base of the thermocline reaches its deepest depth, and likely mixotrophs lose part of their photosynthetic advantage due to the almost complete depletion of nutrients near the surface. Below the photic zone, the community, dominated by Rhizaria, is depth-stratified, relatively constant throughout the year. These populations, relying on the vertical flux of organic matter, respond to local hydrographic and biological features such as the oxygen minimum zone. Together, this suggests a dynamic partitioning of the water column, where the niche vertical position for each community changes throughout the year, likely depending on nutrient availability, the mixed layer depth, and other hydrographic features. Future research should address several uncertainties needed to implement these data into biogeochemical models. The most prominent problem to solve is the lack of lineage-level quantitative data (i.e., counts), which would allow a better understanding of the energy flows. Most of the lineages we report on here are 2-20 µm in cell diameter, falling in the range where taxon-specific cell counting is most difficult, between flow cytometry and automated image analyses. A second uncertainty is clarifying the functionality (trophic modes) in those lineages (e.g. MOCH groups with autotrophs, mixotrophs, and heterotrophs), to more accurately associate these strategies with the temporal and spatial patterns in the system. Together with the present data, these would allow better characterization of the carbon flux throughout the water column, and the role of mixotrophs and heterotrophs in controlling the biogeochemical cycles in the oligotrophic oceans.