Oceanographic data were accessed from the BATS ftp repository¹. Sampling, quality control and post-processing methodologies are available on the BATS website². Temperature, salinity and CTD fluorescence profiles from surface to 1000 m were graphed over time using Ocean Data View Ver. 4.8 (Schlitzer, 2017). Mixed layer depth (MLD) was defined as the depth at which density (sigma-theta; σ_θ) reached 0.125 kg m^–3 difference compared to surface (Levitus, 1982) calculated from CTD profiles. Monthly integrated primary production (as per BATS methodology; Knap et al., 1997) and total Chl-a per square meter were calculated from Dec. 2014 to March 2016. Day and night size-fractionated zooplankton biomass (as per BATS protocol and data) was also analyzed. Both total biomass as well as fractionated biomass were graphed. Sedimentary traps flux data (measuring total biomass, C, N, and P from sedimentary traps during BATS cruises; a trap array, with traps at 150, 200, and 300 m is allowed to drift for 72 h; Steinberg et al., 2001) was also retrieved.

Zooplankton samples were taken during 2015 and early 2016 onboard the RV Atlantic Explorer during BATS cruises from February 2015 to March 2016 (Table 1). Sampling methods and gear have been previously described in detail in Madin et al. (2001) and Steinberg et al. (2012). In short, zooplankton was sampled using a 1-m² rectangular, 202 μm mesh net. Oblique tows were done during the day (between 13:59 and 15:40 GMT) and night (between 01:24 and 6:19 GMT) on each cruise, with the exception of the 20314 cruise in which the day tow was canceled due to rough weather. Net tows were done at a ship speed around 1 kN, and paying out 250 m of wire at 15 m per minute both ways. Depth was recorded using a STAR ODDI DST centi-TD depth and temperature logger (effective depths ranging 113–247 m; Table 1). Once on deck, water was drained, and samples were immediately fixed in 95% undenatured alcohol. Samples were then stored at −20°C until returning to BIOS (1 to 3 days). Once at BIOS, the ethanol was replaced and samples stored at 4°C until analysis. If the ethanol was turbid after 24 h, ethanol was replaced again. Samples were stored in 500 mL of ethanol minimum. If the sample occupied more than 1/3 of the jar, it was then transferred to a 1 L jar, and filled to the top with ethanol.

Prior to extraction, each zooplankton sample was split in two subsamples using a Motoda splitter (Motoda, 1959). Half of the sample was set for extraction, meanwhile the other half was kept as archival material at BIOS. Macroplankton (>2 cm; mostly the larger fish larvae) were removed from the sample and stored in a separate jar, not being used for this study.

The Extraction method was based on the SDS-chloroform extraction in Corell and Rodríguez-Ezpeleta (2014) but using the E.Z.N.A.^® Mollusc DNA Kit (Omega Bio-Tek, Norcross, GA, United States) to recover the DNA instead of the ethanol precipitation. Three final pseudoreplicates (since they originated from the same zooplankton sample) were obtained for each sample.

Samples were filtered through a 53 μm mesh placed on a home-made sieve that allows for easy removal of the mesh without disturbing the sample, and thoroughly washed with milliQ water to remove as much ethanol as possible. The mesh was them removed and placed in a 50 mL falcon tube, with part of the mesh hanging outside the tube. The falcon tube was then tightly closed, securing the mesh to the edge. The tube was then centrifuged at 3,500 g for 10 min to remove as much water and remainder ethanol as possible. After centrifugation, the pellet was weighted, transferred to a fresh 50 mL falcon tube (carefully checking that no plankton was left on the mesh) and 15 mL of SDS buffer (Tris-HCl 10 mM, EDTA 100 mM pH 8.0, NaCl 200 mM, SDS 1%) was added if total mass was < 10 mg, and 25 mL if above. Samples were then vortexed to break the pellet away and ensure rapid mixing of the plankton with the buffer to prevent DNA degradation. The sample was then ground using a Fisher Scientific Homogenizer FSH125, equipped with a 10 mm × 160 mm steel generator. Samples were ground at 1/3 of the maximum speed (approximately 10,000 rpm; this model of homogenizer does not give precise measurement of the speed) for 5 min. This low speed was chosen since, during initial testing, it was noticed that using higher speeds (i.e., near maximum speeds), even for shorter times, often resulted in excessive DNA fragmentation and rapid heating of the sample, and fragments > 2 kb were difficult to obtain during amplification. Large gastropods, whose shells softened after contacting the buffer, did not get broken often, and were cut using dissecting scissors. After the grinding step, the generator was rinsed using 10 mL of SDS buffer to minimize the loss of sample in the generator. Between samples, the generator was completely disassembled, all parts washed with soap under running water, soaked in bleach 10% for 5 min and rinsed in MilliQ water before assembling and running the next sample. For samples < 10 mg of biomass, extraction was done using the 25 mL of buffer used (15 mL initially + 10 mL for cleaning the generator). For samples with > 10 mg, extra 15 mL of buffer were added, to make up to 40 mL of SDS buffer. Using less buffer caused DNA to degrade during extraction during preliminary tests. No heating of the sample was detected during the grinding process at low speeds.

After the grinding protocol, proteinase K (Sigma-Aldrich; St. Louis, MO, United States) was added to the buffer (0.2 mg mL^–1 final concentration), as well as 5 mL of sterilized stainless steel Shot, 3/32′′ Ellipse (Rio Grande, CA, United States). The tubes were then placed in an oven at 60°C for 4 h, vortexing at maximum speed for 30 s every 30 min. After the 4 h, the tubes were centrifuged at 3,500 g for 15 min. Three replicates of 400 μL of the supernatant were transferred to 1.5 mL Eppendorf tubes. Three additional cryovials (2 mL of supernatant each) were stored at −80°C, and the rest of the supernatant was transferred to a clean 15 mL Falcon tube and stored at −20°C. For the three Eppendorf tubes, if the supernatant still contained small particles that would not precipitate at the aforementioned centrifugation forces, the Eppendorf tubes were then centrifuged at 10,000 g for 5 min or until complete precipitation of the particles, and then supernatant transferred to a new tube. After this step, DNA was extracted following the E.Z.N.A.^® Mollusc DNA Kit (Omega Bio-Tek, Norcross, GA, United States), from the addition of chloroform:Isoamyl alcohol 24:1 step on. A modification to the standard protocol included a 14 min centrifugation step to separate the aqueous and the organic phases at this point. From that, only 300 mL of the aqueous phase were carried to the next step, to ensure avoiding the milky interface and the organic phase. All the other steps followed the manufacturer conditions. Final elution was done using 2 × 200 μL of Ambion^® nuclease-free water (not DEPC-Treated; Life Technologies, Foster City, CA, United States). DNA was then stored at −20°C until amplification after quantification (NanoDrop^TM One). DNA output ranged between 16 and 100 μg of DNA per pseudoreplicate.

One PCR reaction was done for each pseudoreplicate (total 3 PCR per sample). Amplification of the V9 region was done using the primers 1389F and 1510R (Amaral-Zettler et al., 2009) integrated in an oligonucleotide consisting of an Illumina adapter, an 8-nt index sequence and a 10-nt pad sequence identical for all oligonucleotides (Supplementary Table S1 in Supplementary File S1). Amplifications were done using the AccuPower^® ProFi Taq PCR PreMix (20 μL reaction; BIONEER; South Korea). PCR steps were performed such as to minimize amplification bias and formation of chimeras, including 300–500 ng of DNA template (to minimize the number of PCR cycles needed and the associated bias), fewer PCR cycles (12 cycles) and lower temperatures (Aird et al., 2011; Ahn et al., 2012; Gonzalez et al., 2012). The premix includes the DNA polymerase, reaction buffer, dNTPs, tracking dye, and a patented stabilizer. To the premix, 1 μL of each primer (10 μM) were added; the amplification protocol consisted in an initial denaturation at 94°C for 3 min, followed by 12 cycles of 94°C for 30 s, annealing at 56°C for 30 s and extension at 72°C for 30 s; a final extension step was done at 72°C for 15 min, and then kept at 4°C. After amplification, PCR success was tested using 3 μL of the PCR reaction in a 1.5% agarose gel in TBE buffer. For each sample, the remainder of the PCR product from the pseudoreplicates were pooled and sent for sequencing at the Microbial Analysis, Resources, and Services (MARS) unit of the University of Connecticut Biotechnology Bioservices Center. Sequencing was performed in an Illumina MiSeq using the MiSeq Reagent Kit v2 (500-cycles; 2 × 250) V2 chemistry. Sequence files are available at the NCBI site under the bioproject https://www.ncbi.nlm.nih.gov/bioproject/PRJNA540654. To test for replicability within the procedure, two samples were submitted using different indexes in separate batches.

Demultiplexed samples were analyzed in MOTHUR ver. 1.39.5. (Schloss et al., 2009). The full script, annotated, is available at https://github.com/blancobercial/BATS_zooplankton. Contigs were made allowing for trimming outside the overlapping region (therefore only reads giving full length in both directions will be retained in a later cleaning step), and bases that are compared to a gap in the other fragment were eliminated if the quality score was below 30. After pairing, all reads containing any ambiguity or shorter than 115 bp were removed. Unique sequences were aligned against the V9 region of the SILVA 128 release database (Quast et al., 2013). Sequences were trimmed to the length of the V9 region, and only those showing completeness (starting in the first base and ending in the last base of the alignment) were kept, avoiding artificial OTUs due to unfinished amplifications. After gap removal, chimeras were removed using UCHIME (Edgar et al., 2011) as implemented in MOTHUR, and unique sequences selected (since after trimming and processing a few sequences may become identical). Distances between reads were calculated, considering gaps as a single event per gap no matter the length. The effect of base by base threshold clustering was studied running an initial clustering step using the nearest neighbor algorithm, with a 0.10 cutoff and 1000 precision. After graphical identification of the clustering level corresponding to 2 bp (= 0.016), final clustering was done using the OptiClust algorithm (Westcott and Schloss, 2017) with distance set to 0.016. This distance would only remove PCR and sequencing errors, since it would only cluster the true reads with the derived reads with 1 bp-error (2 bp distance between two reads with one error each). This approach is similar in concept to the TARA oceans 1 bp swarming procedure (de Vargas et al., 2015). OTUs were taxonomically assigned to a database developed from the complete SILVA 128 release database for SSU (SILVA_128_SSURef_tax_silva; not the NR99 available from www.MOTHUR.org)³. Bash scripts used to transform the SILVA files to MOTHUR format are available in the repository containing all scripts used in this manuscript⁴. Using the NR99 would have worsened the taxonomic assignment, since the 18S is a highly conserved gene, and it was observed that the same V9 region sequence was often shared between components belonging to the same Superfamily. Based on this classification, only metazoan data were selected for all downstream analyses. For identification of taxa (OTUs) flagged by ecological analyses (see below), each of the OTU were BLASTed against the GenBank and results individually analyzed, to take into account the limits of taxonomic assignment of the V9 region and the existing references.

The observed seasonality of the water column has been previously described for the BATS station in several publications (e.g., Steinberg et al., 2001; and references therein; Lomas et al., 2013). Briefly, temperatures showed deep convective mixing early in the year (winter and early spring), moving to a strong summer stratification in surface for July, August and September, with deepening of the mixed layer later in the year (November and December) leading to the deep mixing during winter (Figure 1). During October 2015, Hurricane Joaquin passed near Bermuda (missing data in the graphs) and likely produce surface mixing. No data was registered from the BATS program due to the associated rough weather conditions, and profiles were obtained from Mid-Atlantic Glider Initiative & Collaboration (MAGIC⁵) program (R. Curry, unpublished), revealing abrupt mixing of the column after the passage of Joaquin (Figure 1). Data gaps early in both years correspond to scheduled shipyard repairs for the RV Atlantic Explorer. CTD fluorescence showed higher concentrations during spring bloom in late February (i.e., winter), with some intermittent smaller peaks during the year between 70 and 100 m depth. When considering all casts made during each BATS cruise, the fluorescence peaks roughly follows the measured Chl-a (Figure 2) or Phaeopigments (Turner fluorometer; see BATS methods). MLD showed a patter similar to those already described for BATS (e.g., Church et al., 2013), with a deeper layer in winter, abruptly shallowing in spring and summer and gradually deepening during the fall despite the hurricane (Figure 2). Primary production showed maximum peaks in late winter/early spring and August, with a minimum during the fall (Figure 2) also in agreement with previously published observations for BATS (e.g., Steinberg et al., 2001).

Zooplankton biomass showed maximum values in March, May and August, with night values being higher for monthly pairs with the exception the May and September (Figure 2 and Supplementary Table S2 in Supplementary File S1). In general, each fraction followed the same common trend, with similar maximum and minimum relative values.

All samples succeeded in producing DNA, and amplification success was 100%. A total of 7,199,536 reads were produced, of which 5,144,471 passed quality control and chimera search, with 141,553 unique sequences. After clustering at 2 bp difference, 20,114 metazoan OTUs remained, accounting for 3,781,379 reads. The remaining reads belonged to other Eukaryotes OTUs (SAR, unclassified, etc.) and were not included in downstream analyses. Raw data is provided in an excel file including OTU abundances per sample, OTU taxonomic assignment and OTU representative sequence (Supplementary File S2).

All samples were reduced to 36,471 reads per sample and only core BATS cruises were considered. Taxonomically, calanoid copepods overwhelmingly dominated the community in all samples, representing 45–65% of the community (Figure 5). No other of the taxa groups considered reached values over 13% in any case. At this major taxa group level, no clear seasonal pattern was appreciated, with proportions remaining more or less constant.

After building the Bray–Curtis similarity matrix, a nMDS analyses on the non-transformed data showed little seasonality pattern, however day/night samples were separated to a certain degree, showing statistically significant grouping by this factor despite the partial mixing (ANOSIM R = 0.289, p < 0.005; Figure 6). The taxa that contributed most (up to 25%) to the differences between day and night samples, were typically OTUs with a higher number of reads, in particular OTUs 1, 2, 3, 5, 11, and 7 (Table 2). They contributed more than 4% each to the dissimilarity, accounting for up to 35% of the total dissimilarity. GenBank BLAST analyses indicated that their identity would correspond to (1) the Family Calanidae (Copepoda, Calanoida; several genera and species of the family shared the same V9 region), (2) Pleuromamma spp. (Gaussia princeps shares this sequence, it is however unlikely due to its meso-bathy-abyssopelagic habitat; Metridia does not share this sequence), (3) the Euchaetidae or Aetideidae; (5) unknown members of the Superfamily Clausocalanoidea, (11) the Superfamily Eucalanoidea (sensu Blanco-Bercial et al., 2011; again, the V9 region is common to several genera within this superfamily, across families) and (7) the Family Euphausiidae (Order Euphausiacea; several genera of the shared this sequence) respectively (Table 2). Contribution to the dissimilarity progressively decreased for the other OTUs, and only 20 contribute more than 1% (Supplementary Table S3 in Supplementary File S1).

After standardization by the total and square root transformation, K-R Clustering showed that the maximum improvement in R occurred between 3 and 4 clusters (Figure 7A). The defined groups coincided with the four defined clusters by hierarchical group average method (Figure 7B). The nMDS analyses reflected a seasonal cycle based in those four clusters (Figure 7C), corresponding to seasonal hydrographic regimes: winter deep mixing and bloom (February to April; spring bloom happened in February, Figure 1), spring (May to July; coincident with the largest zooplankton biomass), summer stratification (August and September) and fall mixing event (November and December; post-hurricane Joaquin). Environmental variables (MLD, PP and Chla) showed high correlation with the nMDS ordination and grouping (Pearson ρ = 0.458; p < 0.005). ANOSIM analyses showed high significance for the four groups (p < 0.001) meanwhile it resulted in non-significant grouping of Day/Night (p > 0.25). A mantel test comparing the community composition with measured vertical particle flux measured at BATS (RELATE procedure) showed a statistical correlation between both similarity matrices (ρ = 0.348; p < 0.001). A direct correlation between zooplankton biomass (measured during the day, during the night, or the average of both) did not result in significant correlations when compared with the exported total mass. The raw data for the vertical flux can be found in the Supplementary File S1 in Supplementary Table S4.

SIMPER analyses between the four groups (Table 3) indicate that the transitions between seasons were due to many OTUs, each contributing at the most 1 to 2% to the dissimilarity between seasons (detailed results in Supplementary Table S5 in Supplementary File S1). Taxonomic assignment by BLAST gave a wide range of IDs, from several phyla, covering both holoplankton to meroplankton. The deep mixing and spring bloom to late spring transition was marked by the increase of Eucalanoidea (2 OTUs) and Calanidae copepods, with a decrease of winter taxa such as appendicularians, euphausiidae and Pleuromamma sp. Toward the summer, Poecilostomatoid copepods, and Acartia longiremis increased their relative abundances, meanwhile during the fall there is an increase of an Eumalacostraca and an OTU belonging to the Euchaetidae/Aetideidae clade. The transition to the winter is again characterized by increases in an appendicularian, Pleuromamma sp. and Euphausiidae, and a decrease in the groups dominant in summer and fall. Similarity within groups were caused mostly by the most abundant OTUs (especially 1, 2, 3, 5, 6, 7, 9, and 11). Some other taxa were, however, typical of only one or two of the groups, such as Poecilostomatoid (OTU 18) for the summer and the unknown Eumalacostraca (OTU 15) in the fall. In case of Eucalanoidea, it was possible to appreciate a gradual decrease in their relative contribution to the groups, with a maximum in spring, decreasing toward the summer and fall.

Grouping the OTUs by main taxa groups following the classification used from images carried out from BATS data (Ivory et al., 2018) resulted in a lack of ability to detect seasonality. Samples analyzed this way were clustered into three groups, one containing summer (June to September) samples and the other two the rest of them. Day/night pairs were, however, mixed within groups, and even same-month day/night samples fell in separate groups (Supplementary Figure S2 in Supplementary File S2).

When performing any study, the limitations of the approach should be always carefully considered, and metabarcoding is not an exception (Creer et al., 2010; Pilgrim et al., 2011; Bucklin et al., 2016). The most commonly questioned or criticized facet of the metabarcoding is how quantitative metabarcoding is (if at all). Reviewed recently by Lamb et al. (2019), there is in general a correlation between number of reads and biomass, although with large uncertainties. In the present study, the top 10 OTUs would correspond to taxa or clades that are known to be among the most dominant at the BATS site (Deevey, 1971; Deevey and Brooks, 1977; Steinberg et al., 2000), however no morphological counts were done for the same samples. It is also important to remember that metabarcoding only provides an equivalent to relative abundances/biomass. In our case, if weighing the relative abundances by the biomass for each month and day/night, it was observed that the total maximum peaks for several of those taxa would correspond with the sample of maximum biomass (May; 10314 cruise), and not the maximum of relative abundance (figure not shown). Finding the best way to transform these metabarcoding counts into real abundances/biomass should be one of the priorities for the whole scientific community.

In the pre-processing steps, the most critical step appeared to be the rarefaction to the minimum number of reads. This process has the associated problem of loss of information (Colwell et al., 2012; Chao et al., 2014), however it increased the differences between samples, biasing diversity toward those with more original reads. If this is an effect of artificial inflation during sequencing is something to be studied in the future. The removal of singletons had a stronger effect within the extrapolated data, while the additional removal of doubletons was minimal, especially in the interpolated dataset. All these results suggest that number of taxa were not be inflated due to methodological errors (that would have been marked by a more dramatic decrease in OTUs when removing singletons and doubletons). The approach in this manuscript is very conservative, and our clustering (theoretically) just removed PCR errors, and was equivalent to a 100% similarity matrix after de-noising (some OTUs are still 1 bp different, since they were the centroids of their respective error-derived cluster of reads). It did not, however, result in an excessively high final number of OTUs for the region, considering existing literature. One of the advantage of the V9 region is its short length (∼130 bp), therefore a complete overlapping between F and R reads was obtained. Discarding those pairs or reads that did not completely agree with a Q > 30 might have helped removing PCR and sequencing errors, and likely dropped chimeras. Both are known to be factors that might be causing the apparent inflation in the number of taxonomic units in metabarcoding studies.

Another limitation in this particular study is the taxonomic resolution of the V9 region. Despite the fact that it has been widely used for community characterization in the zooplankton (Amaral-Zettler et al., 2009; de Vargas et al., 2015; Abad et al., 2017), its ability to resolve species is somewhat limited in zooplankton (Wu et al., 2015). For the V9 alone, even identification at the genus level would have only an 80% success for Copepoda (Wu et al., 2015). A similar problem was also detected for euphausiids in the present manuscript (where several species and genera shared the same sequence), and for pteropods (Stewart et al., 2018), suggesting it would likely be for all taxonomic groups, if studied in detail. For Copepoda, it can also be observed that the resolution is not equal across the phylogenetic tree. Within the order Calanoida, the oldest planktonic families, Arietelloidea (= Augaptiloidea Blanco-Bercial et al., 2011) and Diaptomoidea showed good resolution to genus and even to the species level (e.g., the several Pleuromamma OTUs; or the Centropages sp. detected). If considering at the most modern superfamilies, whose planktonic radiation happened more recently (Bradford-Grieve, 2002), the same V9 is shared by whole families (Calanidae) or even across several families within the most modern superfamily, the Clausocalanoidea. This is a widespread issue (Albaina et al., 2016; Abad et al., 2017), and this asymmetry of the resolution would influence estimates of diversity, with a closer accuracy in areas where the Arietelloidea or Diaptomoidea are dominant groups (e.g., the mesopelagic or coastal environments, respectively) compared to niches where the other superfamilies are more common (e.g., the epipelagic open ocean). In the BATS region, Pleuromamma spp. and Lucicutia spp. (Metridinidae and Lucicutiidae; Arietelloidea) are significant contributors to the zooplankton biomass (Deevey and Brooks, 1977; Schnetzer and Steinberg, 2002), however the other superfamilies also contribute significantly to the community (pers. obs.; Deevey and Brooks, 1977). In terms of total number of species, 241 copepod OTUs passed the different thresholds in this study. Compared to morphological analyses, 326 species were found by Deevey and Brooks (1977) from the surface to 2000 m depth, 128 of those in the epipelagic, although they acknowledge many more (especially small ones) were not identified and were left at the genus level. A similar number of morphological species has been obtained (159 species) during summer 2016 in the upper 200 m at BATS (Blanco-Bercial and Maas, unpublished), where likely a higher species-level resolution was attained compared to Deevey and Brooks (1977). The higher numbers obtained in this metabarcoding analysis could be attributed to the fact that a whole year was sampled. Furthermore, in the metabarcoding approach, all of the sample is analyzed, while in the morphological counts a subsampling of some kind is needed in order to make the counting feasible. It would be also possible that, in fact, metabarcoding would be detecting more of the rare community, which in the Sargasso Sea it makes a significant fraction of the zooplankton. Considering the difficulty in attaining species-level ID between congeneric groups in this oceanic region, the expected number of species could easily double that number (around 500 species). That would still not be out of scale, considering the high diversity of this area and the understudied non-calanoid copepods. For the subtropical NW Atlantic, at least 713 named copepod species have been indicated (Razouls et al., 2005–2019). Although many would not be distributed in the Sargasso Sea, most of them could potentially be present in the samples analyzed here. More accordance was found in pteropods (but a very complete reference was available from Burridge et al., 2017) or in Ostracods, when metabarcoding identified 110 OTUs for ∼120 species described for the subtropical North Atlantic (Angel et al., 2008). Considering these three groups, this dataset does not seem to have inflated the existing diversity, and might be relatively well-representing the community composition.

Understanding the limits (and the capabilities) of metabarcoding is essential to make use of its full potential, while avoiding over interpretation of the results. The final goal of the metabarcoding should not to be viewed as the complete and perfect description of the community. Lacking complete databases, and working with a marker that is neither able to discriminate among nor efficiently amplify all species, among others, are limitations that, at least at present, seem unavoidable. These limitations, and many others, have been discussed in several reviews (Creer et al., 2010; Pilgrim et al., 2011; Bucklin et al., 2016). In the present study, however, metabarcoding effectively revealed community regimes, with correlation to environment and biogeochemical processes, providing a step toward a better understanding of ecological succession and processes in the Sargasso Sea zooplankton, and allowing for the development of new working hypothesis regarding the zooplankton community at the BATS site. Many potential studies are furthermore possible based on this dataset; e.g., biological interactions (Lima-Mendez et al., 2015; Stewart et al., 2018), in depth studies in relation to biogeochemical cycles (Guidi et al., 2016), among others. The equivalent analyses of the zooplankton community, done by classical methods (morphology), would have needed a much larger effort (time-wise) and several experts while the faster imaging (zooSCAN or similar) analyses do not provide the same kind of information. It is time to incorporate metabarcoding as a standard procedure in time-series, environmental monitoring, and regular research, not as a substitute of, but as a complement to, the present methods used in biological oceanographic research.