Discussion
The data presented here supports the presence of a significant seasonality in the photic community in the Subtropical North Atlantic, whose signal and strength gradually decreases with depth, becoming undetectable in the mesopelagic. The structure of epipelagic and mesopelagic protist communities in the Sargasso Sea is strongly shaped by local ocean dynamics (e.g., mixing, stratification, nutrient distributions and light penetration) and expressed in the vertical and seasonal partitioning of community composition, diversity, and function. Mixotrophic groups dominate the sunlit epipelagic layers, with autotrophic taxa gaining relevance where the base of the photic layer meets the top of the nutricline (Arenovski et al. 1995). The mesopelagic layers below are inhabited by heterotrophic communities where Rhizaria becomes the most abundant group, outside of parasitic lineages (i.e. Syndiniales). The deep communities are likely fueled by complex food webs, feeding on prokaryotes, particles and each other (Byung Cheol et al. 2000; Rocke et al. 2015).
Nutrient limitation in the epipelagic zone favors mixo- and heterotrophic protists closer to the surface, with autotrophs balancing nutrients and light at its base. This is largely a consequence of ocean stratification, which enhances or inhibits vertical exchanges. The degree of stratification at the base of the photic zone affects the rate of nutrient supply fueling the autotrophic lineages at the local DCM. Over the annual cycle, hydrographic Layers 0, 1 and 2 exhibit dramatic changes in stratification driven by local air-sea fluxes, with Layers 1 and 2 being subsumed into Layer 0 beginning in the Fall and ending at the Spring transition. This process creates the vertical partitioning of clusters 1, 2 and 3 and their seasonality observed in K=9: i.e. clusters 1 and 2 dominate Layers 0 and 1 above the DCM from Spring through the Stratified season, with cluster 1 disappearing in the Fall and Mixed seasons, while cluster 3 in Layer 2 waxes in Spring then gradually wanes through the Stratified and Fall seasons. Such seasonality is absent at ALOHA in the North Pacific where winter mixing is shallower, the base of the photic layer remains stratified throughout the year, and protist communities align with depth (Ollison et al. 2021). This difference is not unexpected, due to the more constant hydrographic conditions, primary production and particle flux measured throughout the year by the Hawaii Ocean Time-series compared to those at BATS (Churchet al. 2013). The more homogeneous conditions would then favor the presence of a dominant community all year long at ALOHA, showing only depth-controlled layering.
The pronounced seasonality observed in the Sargasso Sea implies shifts in trophic ecology and energy fluxes which likely influence biogeochemical cycles over the annual cycle. Elevated abundances of mixotrophs and heterotrophs relative to autotrophs during the stratified and fall seasons indicate a complex recycling food web where small protists such as Warnowia , Telonemia , Karlodiniumor Minorisa minuta and MAST and MOCH lineages prey on picophytoplankton, responsible of most of the primary production in the Sargasso Sea (Caron et al. 1999; Cotti-Rausch et al. 2020; del Campo et al. 2013; Klaveness et al. 2005; Massanaet al. 2014; Orsi et al. 2018; Place et al. 2012; Riemann et al. 2011; Sanders et al. 2000). Both autotrophs and small mixo- and heterotrophs would also be preyed upon by larger mixotrophs or heterotrophs protists (Andersen et al. 2011; Evelyn & Michael 1998; Quevedo & Anadón 2001). Autotrophic protists (picophytoplankton and small nanophytoplankton, such asOstreococcus , Bathycoccus or Pelagomonas calceolata ) had a numerical relevance only in the mixed and spring seasons, but mixotrophs still represent a significant portion of the community. Those mixotrophs would benefit from their ability to carry out photosynthesis, while still preying on the small autotrophs (Sanderset al. 2000), likely due to the concurring C fixation by photosynthesis with the heterotrophic acquisition of N and P (Edwards 2019). Recently, Choi et al. (2020) reported elevated summer abundances of uncultured dictyophytes at BATS, and evidence that these taxa are more prevalent in the upper water common when nutrients are most depleted. These and other observations, and the findings we report, suggest that mixotrophic strategies become relatively more advantageous when production is limited by macronutrients, likely due to the repartitioning of N and P rather than the energetic benefits of mixotrophy. As such, in the oligotrophic Sargasso Sea they would be at a significant advantage compared to non-mixotrophs (Choi et al.2020; Duhamel et al. 2019; Edwards 2019), only attenuated during the strong mixing period.
The observed month-to-month depth expansion and contraction of the different protist communities influences community functionality. Cluster 1, which occupies the near surface layers throughout the spring and stratified season, disappears in fall when the surface MLD begins to erode into the DCM, and only emerges again when the ML abruptly shoals the following spring. The temporary nature of this group contrasts with the other communities, which are likely present throughout the year, albeit with periods of expansion and contraction. As such, heterotrophy, and its trophic and biogeochemical consequences, would increase in summer and fall, with the expansion of the stratified surface community. The fall is also when the Chl-a -max corresponds to a community rich in heterotrophs, although many are potential mixotrophs via endosymbionts (Bjorbækmo et al. 2020). The effect of these non-constitutive mixotrophs on the biogeochemical cycles is one major question to be addressed in the future. Although in some cases a few clusters are missing in some months, it might be an artifact of the punctuated depth sampling, which would miss their depths at that particular month. As a comparison, the bottle sampling also missed some of the density-derived layers, but they were all present in the continuous profiles. In those months, some of the deeper clusters were then missed.
It should be noted that in some cases the functionality of taxa is relatively easy to assign (e.g. Telonemia orOstreococcus ), while in many other cases it is difficult to assess with certainty. Close relatives, with different functional profiles (i.e. mixotrophy vs. autotrophy) often share a single V4 sequence. Other sources of uncertainty come from the non-constitutive mixotrophs. In our case, the main group is Rhizaria. These are known to potentially carry photosynthetic endosymbionts (e.g., Phaeocystissp.) (Bjorbækmo et al. 2020). Addressing these uncertainties is required to understand the full role of mixotrophy in the oceans, and its influence in the biogeochemical cycles (Edwards 2019; Flynn et al. 2019; Gonçalves Leles et al. 2021; Gonçalves Leles et al. 2018; Mitra et al. 2014).
In the mesopelagic zone, protist are vertically partitioned into distinct communities that occupy specific density strata. This may reflect a response to resource availability and/or particular environmental conditions (Biard & Ohman 2020), such as the OMZ, where there is an increase in particles. Within the main clades, specific lineages also occupied specific vertical layers, suggesting a preferred niche for each. An alternate explanation is that vertical biological gradients may reflect “where” and “when” the density layers were at the sea surface – i.e. the ventilated thermocline model of gyre circulation (Luyten et al. 1983). The “age” of a water parcel at BATS typically increases with depth/density reflecting longer distances traveled since ventilation. Geographic gradients in surface community composition could thus translate to vertical gradients at mesopelagic depths, and/or, if the community composition evolves as a function of time, the vertical gradients may reflect “age”.
Diversity patterns
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 in DNA metabarcoding, while Ollison et al. (Ollison et al. 2021) used an RNA-based 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 Chla-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; here light and nutrients are able to sustain a rich autotroph community (below it is too dark; nutrients are scarce above). 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 overlapping presence of light and nutrient gradients would favor the presence of more, and more diverse, ecological niches, facilitating the coexistence of more phylogenetically apart lineages.
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.
Methodological caveats
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 PR2 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, and would very likely hold true independently of the marker used. Finally, diversity indices (the raw value itself) 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 (Blanco-Bercial 2020; Santoferrara 2019). Again, despite any possible bias, it is likely that the observed patterns of variability (where diversity is higher/lower, etc.) are robust against those bias. Any comparison with other existing data would require a reanalysis with shared bioinformatics approaches.