Patterns of species diversity of the many taxonomic groups of marine zooplankton are key characteristics of ocean ecosystems, determining their function, sustainability, and responses to environmental variation and anthropogenic impacts, including climate change (Sherman et al., 2002; Friedland et al., 2020). Pelagic ecosystems of the NW Atlantic Ocean have been monitored and studied over many decades, providing an invaluable time-series record of zooplankton diversity and abundance (Wiebe et al., 2012; O’Brien et al., 2013).

DNA metabarcoding of zooplankton samples has been used in association with ecosystem monitoring and management in recent years (Mohrbeck et al., 2015; Deagle et al., 2017; Djurhuus et al., 2018; Blanco-Bercial, 2020; Matthews et al., 2021). These efforts have used different molecular protocols and bioinformatics pipelines (Bucklin et al., 2021a). A subset of the studies have focused on discrimination and identification of species diversity across the zooplankton assemblage based on metabarcoding analysis using portions of the mitochondrial cytochrome oxidase I (COI) gene (reviewed by Bucklin et al., 2021b), which has been widely used to identify species of marine organisms (Bucklin et al., 2011).

The potential of COI metabarcoding has been examined from many perspectives, including accuracy of species-level resolution of biodiversity (Brown et al., 2015; Leray and Knowlton, 2017; Schroeder et al., 2021), availability of reference sequence databases resulting from DNA barcoding efforts (Andújar et al., 2018; Porter and Hajibabaei, 2018; Bucklin et al., 2021b; Singh et al., 2021), and prospects for quantitative analysis related to species abundance and/or biomass (Lamb et al., 2018; Matthews et al., 2021). Challenges and disadvantages of COI as a metabarcode include lack of universal primers and missing groups due to primer-mismatch (Deagle et al., 2014; Clarke et al., 2017; Hajibabaei et al., 2019). A number of studies have evaluated the results using COI metabarcodes based on parallel analysis using multiple metabarcode gene regions, including hypervariable regions of 18S rRNA (Djurhuus et al., 2018; Steffani et al., 2018; Giebner et al., 2020; Pitz et al., 2020; Brandão et al., 2021; Pappalardo et al., 2021; Questel et al., 2021; Zhao et al., 2021).

An important consideration for species identification based on COI metabarcoding is the availability of COI barcode sequences for identified specimens (Andújar et al., 2018; Porter and Hajibabaei, 2018; Steinke et al., 2021). The goal of taxonomically-complete COI reference databases for marine zooplankton, with search capacities for targeted taxonomic groups and geographic regions (Bucklin et al., 2021b), is a priority for an international collaborative effort sponsored by the Scientific Committee for Oceanic Research (SCOR), MetaZooGene (WG157; see https://metazoogene.org/).

There is broad interest and considerable enthusiasm for potential applications of metabarcoding for ecosystem assessment and fisheries management (Bourlat et al., 2013; Ji et al., 2013; Kelly et al., 2014; Kelly, 2016; Goodwin et al., 2017; Aylagas et al., 2018). A particular focus is the importance of rapid, accurate, and reliable species-level characterization of time-space variability of the taxonomically-complex pelagic assemblage (Andújar et al., 2018).

Time-series monitoring of the NW Atlantic continental shelf was established by the Northeast Fisheries Science Center (NEFSC) of the National Marine Fisheries Service (NMFS) in 1977 and continues into the present (Cox and Wiebe, 1979; Hare and Kane, 2012). Surveys are scheduled to sample six times each year throughout four regions ( Figure 1 ), with collection of environmental (hydrographic) data and samples at a stratified-random selection of standard station locations (see https://www.fisheries.noaa.gov/new-england-mid-atlantic/ecosystems/monitoring-ecosystem-northeast). Morphological microscopic examination of zooplankton samples from NEFSC surveys has allowed analysis and interpretation of temporal and spatial patterns of variability across the region (Kane, 2007; Kane, 2011; O’Brien et al., 2013). The time-series records have provided clear evidence that the region is experiencing rapid climate change (O’Brien et al., 2013; Friedland et al., 2020). Regime shifts (i.e., persistent changes in the structure and function of ecosystems) have been documented during the 1990s and 2000s in pelagic community structure (Pershing et al., 2005; Walsh et al., 2015; Morse et al., 2017), including zooplankton diversity (Head and Sameoto, 2007; Record et al., 2010; Johnson et al., 2011; Bi et al., 2014). NEFSC time-series records revealed a marked increase in zooplankton displacement volume on Georges Bank (GB) and the Gulf of Maine (GoM) in ~1990 (O’Brien et al., 2013; Figure 2 ). Another regime shift was evident in ~2000, when the earlier changes were partially reversed. In the Gulf of Maine (GoM), zooplankton diversity increased markedly during the early 1990s and decreased rapidly about 2000 (Record et al., 2010; Johnson et al., 2011).

Monitoring efforts in North Atlantic pelagic ecosystems (Wiebe et al., 2012; O’Brien et al., 2013) have provided an essential foundation for recognizing and understanding regime shifts (Borja, 2014; Conversi et al., 2015; Morse et al., 2017; Stern et al., 2018). Regular, standardized, and sustained analysis of zooplankton diversity and abundance based on morphological microscopic examination of zooplankton samples has also provided opportunities for evaluating the accuracy, reliability, and power of rapidly-developing approaches to molecular analysis of zooplankton biodiversity, including DNA metabarcoding (Bucklin et al., 2016; Bucklin et al., 2019).

The EcoMon zooplankton database includes records from morphological (microscopic) taxonomic counts reported as numbers per 10m² for zooplankton taxa, which are identified to species when possible (Kane, 2007; Kane, 2011). A total of 186 zooplankton species of 14 taxonomic groups of metazoan zooplankton have been detected in net samples collected during surveys of the NW Atlantic continental shelf since 1977; of these, 43 species of 9 taxonomic groups have been recorded at >1% frequency of occurrence among all samples (David Richardson, pers. comm.). An additional 237 taxa, including groups of congeneric species, genera, and higher taxonomic groups, are are listed in the EcoMon database, https://www.ncei.noaa.gov/archive/accession/0187513; accessed May 23, 2021 (NMFS/NEFSC, 2019).

Molecular analysis of EcoMon samples has been carried out since 2000, first in partnership with the international project, ZooGene (http://www.zoogene.org/), and from 2004 to 2010 with the Census of Marine Zooplankton (CMarZ; http://www.cmarz.org/), an ocean realm field project of the Census of Marine Life (CoML; Bucklin et al., 2010). The partnership, with collection and preservation of samples for molecular analysis, continues today in association with another international program, MetaZooGene (https://metazoogene.org/) Working Group 157 of the Scientific Committee for Oceanic Research (SCOR).

Bucklin et al. (2019) analyzed 27 EcoMon samples collected from 2002 – 2012 using DNA metabarcoding of the V9 hypervariable region of 18S rRNA, with sequences classified into 28 taxonomic groups of zooplankton. The conserved nature of the 18S rRNA gene allows detection of taxa across the spectrum of marine zooplankton, but does not accurately resolve or identify species (Blanco-Bercial, 2020; Govindarajan et al., 2021; Questel et al., 2021). Bucklin et al. (2019) reported significant positive correlations between V9 18S rRNA sequence numbers and microscopic counts for 7 taxonomic groups for which both types of data were available, with significant regressions for three groups: Calanoida, Gastropoda, and Chaetognatha. These results provided promising evidence that DNA metabarcoding using V9 18S rRNA can provide accurate classification and relative quantification for targeted zooplankton groups, which are important goals for applications for ecosystem monitoring (Lamb et al., 2018; Matthews et al., 2021).

Zooplankton samples for this study were collected by the NOAA Northeast Fisheries Science Center (NEFSC) Oceans and Climate Branch during surveys by the Ecosystem Monitoring Program (EcoMon) of the NW Atlantic continental shelf (Kane, 2007; Kane, 2011; Hare and Kane, 2012; Bucklin et al., 2019). Surveys are designed to sample four regions of the shelf ecosystem ( Figure 1 ). Samples for morphological taxonomic analysis were collected following a standard protocol (Richardson et al., 2010), with both day and night sampling using a 61-cm bongo net fitted with a 333-μm mesh net; oblique tows were a minimum of 5-min in duration and sampled from the surface to within 5 m of the seabed or to a maximum depth of 200 m. A mechanical flowmeter was fitted in the mouth of each net to record the volume sampled. Samples were preserved in 5% formalin and archived at the NEFSC.

Zooplankton samples for genetic analysis were collected during EcoMon survey cruises at 5 randomly-selected locations in each region. Sampling was done using a 20-cm bongo net with 165-μm mesh nets, which was attached to the same cable and deployed with the 61-cm bongo nets. Differences in the opening diameter and mesh size of the nets used for collection of samples for genetic analysis may have resulted in differences between the sets of samples, but the methods were unchanged across all years, regions, and stations, and the resulting time-series patterns of variability were evaluated with this caveat. Samples were preserved immediately in 95% undenatured ethanol, which was changed 24 hr after collection. Samples were transported to and archived at the University of Connecticut, with long-term storage in walk-in freezers (-20°C)

Sample sorting and identification was done at the Morski Instytut Rybacki Plankton Sorting and Identification Center (Szczecin, Poland). Zooplankton samples were split to an aliquot containing approximately 500 specimens; individuals were sorted, counted, and identified to the lowest possible taxon (Kane, 2007; Kane, 2011). Data recorded include abundance measured by area (conc/10m²) and volume (conc/100m³) for selected taxononomic groups and species of zooplankton, including fish larvae. The morphological species count data were downloaded from: https://www.ncei.noaa.gov/archive/accession/0187513, accessed May 23, 2021 (NMFS/NEFSC, 2019).

A total of 27 samples was selected for metabarcoding analysis, including one sample collected in each of three regions, Georges Bank (GB), Gulf of Maine (GoM), and Mid-Atlantic Bight (MAB), during May/June of 2002 – 2012 ( Figure 1 ; Table 1 ). The samples are the same ones analyzed by metabarcoding using the V9 hypervariable region of 18S rRNA by Bucklin et al. (2019). There were a number of sampling gaps due to cancelled cruises, bad weather, and other causes: no samples were analyzed for MAB in 2003 or GoM in 2006; no samples were analysed for 2008; and only a single GB sample was analysed for 2012; COI metabarcoding data are missing for DE1105-25; EcoMon count data are missing for DE0305-38 and DE1105-127 ( Table 1 ). The collection site of one sample (AL0605-53 #13) is correctly shown on GB, although was listed within the EcoMon region of Southern New England (SNE).

Samples were quantitatively sub-divided using a box splitter (Motoda, 1959) to reduce zooplankton volume to ~25 mL. The sample was then washed with distilled water; inserted into a 50 mL Falcon tube above 35 µm Nitex mesh, which served to suspend the material and dry the pellet; and centrifuged at 3500 g for 4 min. The pellet was moved to a new 50 mL Falcon tube, and SDS buffer (Tris-HCl, 10 mM; EDTA, pH 8.0, 100mM; NaCl, 200mM; SDS 1%) 3 mL or equal to pellet volume, whichever was smaller) was added. The sample was homogenized using a hand-held homogenizer (D1000, Thomas Scientific) with saw tooth blade for 4 min at level 5. Proteinase K (MP Biomedicals) was added (0.2 mg/mL of sample) and tubes were incubated overnight in a water bath at 55-56°C. After centrifugation (3500 g for 15 min), 400 uL of the supernatant was transferred to individual sterile 2 mL Eppendorf tubes for storage as necessary at -20 or -80°C. Total genomic DNA was extracted using the E.Z.N.A Mollusc DNA kit (Omega Bio-tek) following manufacturer instructions. All samples yielded DNA of sufficient quantify and quality for metabarcoding analysis. Total genomic DNA was quantified on a Thermo-Fisher NanoDrop 2000 and normalized to a final concentration of 5 ng/µl.

Purified DNA was used to amplify a 313 base-pair (bp) region of mitochondrial cytochrome oxidase I (COI) using the primers: mlCOIintF and jgHCO2198 (Geller et al., 2013; Leray et al., 2013). Forward and reverse primers were altered for multiplexed sequencing by adding 5’ adapters (Illumina, Inc., San Diego, CA). The PCR reaction used 20ng of DNA, with Platinum Taq reagents, 4µL buffer, 2.4µL MgCl, 0.8 µL dNTPs, 0.2µL HiFi Taq Polymerase, and 0.8µL of each primer (10µM), with the following protocol: one denaturation cycle at 94°C for 60 sec; 38 cycles of 94°C for 30 sec, 46°C for 30 sec, 72°C for 90 sec; a single extension cycle of 72°C for 5 min; and an infinite hold at 4°C. COI amplicons were checked for successful amplification by running in a 2% agarose gel with a 50 bp marker.

Library preparation entailed adding index primers in a second PCR amplification of the purified amplicons using a master mix composed of (per sample): 5.0 µl purified PCR product; 5 µl Nextera XT Index 1 Primer; 5 µl Nextera XT Index 2 Primer; 25 µl 2x KAPA HiFi HotStart ReadyMix; 10 µl PCR-grade water; for a total volume of 50 µl. The PCR protocol was: 95°C for 3 min; 8 cycles of: 95°C for 30 sec, 55°C for 30 sec, 72°C for 30 sec; and 1 cycle of 72°C for 5 min. The indexed PCR product was purified using AMPure XP beads, with a final elution volume of 25 µL. Successful library attachment was verified using an TapeStation 4200 D1000 High Sensitivity assay (Agilent Technologies). Libraries were quantified using a Qubit 3.0 fluorometer, normalized according to amplicon size, pooled, and denatured with 0.2 N NaOH. Bi-directional sequencing was carried out at the University of Connecticut Center for Genomic Innovation (CGI; https://cgi.uconn.edu/) using an Illumina MiSeq sequencer using the MiSeq Reagent Nano Kit Ver.2 (500 cycles; 1 million clusters) spiked with a minimum of 20% PhiX (Illumina, Inc.). All samples were analyzed in a single MiSeq run, including negative controls and replicate samples, to allow full intercomparison of samples (e.g., sequencing depth).

Demultiplexed reads for the COI region were processed using a custom script for the Mothur pipeline (Ver. 1.44.3; Schloss et al., 2009) and run on the Xanadu computing cluster of the UConn Computational Biology Core (CBC; https://bioinformatics.uconn.edu/). Contiguous sequences (contigs) were assembled from forward and reverse Illumina MiSeq reads and trimmed to the overlapping section. Sequences were trimmed to a uniform length by removing the beginning and terminal ends of sequences that extended beyond the targeted COI gene region. Sequences containing ambiguous bases, quality Phred scores < 30, and with lengths shorter than 150 bp were removed from analysis. Concerns that PCR error may contribute to errors in biodiversity assessment (Kelly et al., 2019) were addressed by using the UNOISE method (Edgar and Flyvbjerg, 2015) within Mothur (Ver. 1.44.3) to de-noise aligned sequences. Sequences were screened for chimeras using the VSEARCH command (Rognes et al., 2016); sequences with chimeras were removed from analysis.

Among the 27 samples used for COI metabarcoding, 12 were selected at random for sequencing in the same MiSeq run, using a second aliquot from the PCR product for the target COI metabarcode region. The results from these samples were treated as technical replicates and examined using the Wilcoxon rank sum test carried out in MatLab (Ver. 2020B), which indicated no statistical differences among the replicates. In all cases, data from only one of each pair technical replicates were used for definitive statistical analysis.

Taxonomic identification of COI metabarcode sequences was determined using a naïve Bayesian classifier algorithm in Mothur (Ver. 1.44.3). Taxonomic assignments for species-level identifications used bootstrap values ≥ 97% after 100 iterations. Before performing zooplankton community analyses, sequences with abundances < 2 (i.e., global singletons) across the entire dataset were removed. Taxonomic classification and species identification were based upon the MetaZooGene Barcode Atlas and Database (MZGdb; https://metazoogene.org/MZGdb), which includes publicly available COI barcode sequences downloaded from GenBank and BOLD (Bucklin et al., 2021b). The results reported in this study used the North Atlantic Atlas and Database (https://metazoogene.org/MZGdb-NATL), (which includes barcodes for all holo- or mesozooplankton species reported from the region. To ensure completeness of the regional database, DNA barcodes are included in the North Atlantic MZGdb, even if the specimen used for DNA sequencing was obtained from a different ocean region.

The MZGdb allows targeted searches by taxonomic groups and geographic regions of interest and provides the capacity to map and visualize the geographical distribution of species observations and collection locations of specimens used for DNA sequencing on global to regional scales. The MZGdb is a collaboratively-developed product of MetaZoogene (SCOR WG157; https://scor-int.org/group/157/) and the Coastal & Oceanic Plankton Ecology, Production, & Observation Database (COPEPOD, https://www.st.nmfs.noaa.gov/copepod/). The MZGdb builds upon a taxonomically-arranged database of zooplankton observations, biometric traits, photographs, and DNA barcode data in COPEPEDIA (https://copepedia.org/), which stores and compiles information at multiple taxonomic levels, including species and taxonomic groups.

COI sequence numbers are reported for species of the same 17 taxonomic categories used to classify metabarcoding results for V9 18S rRNA by (Bucklin et al. 2019) ( Supplementary Table S1 ). Morphological (microscopic) abundance counts per 10m² were recorded for these same taxonomic groups ( Supplementary Table S2 ). Results were further analyzed for six taxonomic groups (Calanoida, Cyclopoida, Eucarida, Chaetognatha, Hydrozoa, and Gastropoda) for which both metabarcoding results and microscopic taxonomic counts were available for most samples; groups with zeros or missing data for many of the samples were not included in the analyses. All data were transformed (Log10+1) prior to analysis.

Numbers of COI sequences and abundance counts from morphological taxonomic analysis reported in the EcoMon database were statistically compared for the six groups using functional regression analysis (Ricker, 1973). Multivariate statistical analyses of COI sequence numbers for the six zooplankton groups were carried out to examine patterns of variation between regions and years using MatLab (Ver. 2020B). One distance measure used was Bray-Curtis dissimilarity coefficient (Bray and Curtis, 1957; McCune et al., 2002), with results displayed by cluster diagrams. Differentiation among the 3 regions was evaluated by Non-Metric Multidimensional Scaling (NMDS) and by Nonparametric (Permutation-based) MANOVA using the FATHOM Toolbox for MatLab (Jones, 2017; https://www.usf.edu/marine-science/research/matlab-resources/index.aspx/). The Shannon (H) and Simpson (D) Diversity Indices (Pielou, 1977) were calculated using COI sequence numbers for the 6 taxonomic groups for each sample. Regression analysis of values for the two indices showed no difference in all cases and results are reported for the Shannon Index (H).

COI metabarcodes identified a total of 181 species across 23 taxonomic groups of metazoan zooplankton listed in the summary file (Wang et al., 2007) generated by Mothur (Ver. 1.44.3; Schloss et al., 2009). Many of these species showed very low sequence numbers including many zeros ( Supplementary Table S3 ), which prevented statistical analysis. Multivariate statistical analysis was carried out for 23 species with total COI sequence numbers >1,000 across all samples ( Table 2 ). All data were transformed (Log10+1) for analysis. Patterns of variation in COI sequence numbers for the 23 species were statistically evaluated between regions and years in MatLab (Ver. 2020B), using the same tests as for the group comparisons, including NMDS and Nonparametric (Permutation-based) MANOVA (Jones, 2017; https://www.usf.edu/marine-science/research/matlab-resources/index.aspx/), Bray-Curtis dissimilarity coefficient (Bray and Curtis, 1957; McCune et al., 2002), and Shannon Diversity Index (H; Pielou, 1977).

Morphological counts are available in the EcoMon database for 7 of the 23 species with >1,000 COI sequences; the remaining 16 species were either grouped at a higher taxonomic level or were not listed among species to be identified (NMFS/NEFSC, 2019; https://www.ncei.noaa.gov/archive/accession/0187513). Numbers of COI sequences and miscroscopic counts were statistically compared for these species using functional regression analysis (Ricker, 1973).

Metabarcoding using a portion of the COI barcode region was carried out for 27 samples collected from three EcoMon Survey regions (GB, GoM, MAB) during 2002-2012 ( Table 1 ). The taxonomic groups selected were used in a previous metabarcoding study by Bucklin et al. (2019) that analyzed the V9 hypervariable region of 18S rRNA for the same EcoMon samples. Metabarcoding yielded a total of 4,992,468 COI sequences and 1,404,242 Amplified Sequence Variants (ASVs). Considering all samples together, numbers of COI sequences and ASVs for species with >50 sequences ( Supplementary Tables S3 , S4 ) were highly significantly correlated across all taxonomic groups (r = 0.978, p = 6.102 e^-11). The definitive analysis used COI sequence numbers.

Statistical analysis focused on 6 taxonomic groups (Calanoida, Cyclopoida, Eucarida, Gastropoda, Hydrozoa, Chaetognatha) for which sufficient numbers of observations (non-zero) were available for both molecular (metabarcoding) and morphological (counts) data. These groups were the same ones analyzed for V9 18S rRNA metabarcodes by Bucklin et al. (2019) for these same samples; a seventh group (Peracarida) could not be analyzed due to many zero values for COI sequence numbers ( Supplementary Table S1 ). Functional regression analysis (Ricker, 1973) of COI sequence numbers versus morphological counts was statistically significant for five of the six taxonomic groups; the exception was the Hydrozoa ( Figure 3 ).

Interannual and regional patterns of diversity of the 6 taxonomic groups based on COI sequence numbers and morphological counts analyzed by NMDS (Jones, 2017) revealed similar patterns, with considerable overlap among the 3 EcoMon regions, but some evidence of distinctive samples in some regions and years for both molecular and morphological analysis ( Figure 4 ). Based on Nonparametric (Permutation-based) MANOVA (Jones, 2017) analysis, the 6 groups showed different patterns of variation among regions and years (p < 0.001), but samples did not differ statistically significant among regions (p = 0.772). Cluster diagrams based on the Bray-Curtis dissimilarity coefficient (Bray and Curtis, 1957; McCune et al., 2002) based on both COI sequence numbers and morphological counts showed clear differentiation of MAB samples, with some overlap between GB and GoM ( Figure 5 ). There are two groups of GB stations based on both COI sequence numbers and EcoMon counts ( Figure 5 ), although the two GB clusters grouped together for EcoMon counts ( Figure 5B ), while GB #10-14 grouped with MAB and GB #15-19 grouped with GoM for COI sequence numbers ( Figure 5A ).

The Shannon Diversity Index (H; Pielou, 1977) showed lower levels for GoM during 2002-2005, based on both COI sequence numbers and morphological counts, with more variation among the 3 regions based on COI sequence numbers for 2007, 2009, 2010 and 2011 ( Figure 6 ). The Simpson Index was also calculated from the same data, with results that were statistically nearly identical to the Shannon Index based on regression analysis (COI sequences: r = 0.965, p = 2.138 e^-15; EcoMon counts: r = 0.982, p = 2.636 e-¹⁷).

COI metabarcodes identified 181 species across 23 groups at varying taxonomic levels; 67 species of 15 groups had > 50 COI sequences; 23 species had >1,000 COI sequences ( Supplementary Table S3 ). Classification and identification of species based on COI sequences used the North Atlantic regional MetaZooGene Atlas and Database (MZGdb), which includes 75,976 barcodes for 12,985 zooplankton species reported to occur in the region (https://metazoogene.org/mzgdb-natl, accessed March 13, 2022). Selection of species for the MZGdb is based upon the COPEPOD database (https://www.st.nmfs.noaa.gov/copepod/) and is designed to ensure an accurate, reliable, and taxonomically-complete reference sequence database with appropriate geographic coverage (Bucklin et al., 2021b).

Multivariate statistical analysis of metabarcoding results focused on 23 species with COI sequence numbers >1000 across all samples ( Table 2 ). Two-dimensional NMDS results for the 23 species showed grouping of samples for each EcoMon region, with a distinct cluster of GoM samples and overlap between MAB and GB samples, with the notable exception of the 2010 MAB sample ( Figure 7 ), which was dominated by the siphonophore, Nanomia cara ( Supplementary Table S3 ). The 23 species showed some variation among regions and years (p = 0.038) based on Nonparametric (Permutation-based) MANOVA (Jones, 2017) analysis, but did not differ statistically significant among regions (p = 0.359). The Bray-Curtis similarity index cluster diagram based on COI sequence numbers for 23 species also showed clear regional patterns, except for the 2002 GB sample, which clustered with MAB ( Figure 8 ).

The pattern of variation of the Shannon Diversity Index (H) based on COI sequence numbers for 23 species showed similar patterns of year-to-year variation in each of the 3 EcoMon regions, with marked differences between 2002-2005 versus 2007-2011 ( Figure 9 ), similar to the interannual variation of H based on COI sequence numbers for the 6 groups ( Figure 6 ).

Functional regression analysis of COI sequence numbers versus morphological counts from the EcoMon database were done for 7 species for which both types of data were available for most stations. Of these, 5 species showed significant regression relationships: Calanus finmarchicus, Centropages hamatus, C. typicus, Pseudocalanus spp., and Temora longicornis ( Figure 10 ); regressions were not significant for two species (not shown): Acartia longiremis (r = 0.607, p = 0.277) and Nannocalanus minor (r = 0.439, p = 0.711). In all cases, the species showing significant regressions had higher numbers of COI sequences, higher counts, and had data of both types for more stations; the species with insignificant regressions had more missing observations and recorded zeroes.

The NW Atlantic continental shelf was designated as a Large Marine Ecosystem (LME) based on the importance of the region for commercial harvesting and the need for conservation measures (Sherman et al., 2002). Despite many challenges over recent decades, including rapid warming from climate change (Friedland et al., 2020), the region has remained an important and productive region for commercial harvesting of numerous species. The importance of the pelagic community, and the zooplankton assemblage in particular, in ecosystem function and services has been acknowledged and examined for many decades (Sherman and Duda, 1999; Walsh et al., 2015; Friedland et al., 2019). Marked differences have been observed among the regions of the NW Atlantic continental shelf ecosystem ( Figure 1 ) in temporal patterns of variation in ecosystem dynamics, including zooplankton diversity and biomass (O’Brien et al., 2013; Figure 2 ).

The importance of biodiversity in the functioning of marine ecosystems is well established (Gamfeldt et al., 2015). A number of studies have examined the more specific question of the role of species diversity of zooplankton in sustaining ecosystem services, including commercial fisheries (Byron and Link, 2010; Bi et al., 2014; Morse et al., 2017). Analysis of the taxonomic composition, diversity, abundance, and biomass of the zooplankton assemblage can serve as an early indicator of climate impacts and regime shifts in the region (Johnson et al., 2011; Borja, 2014; Stern et al., 2018).

A number of previous studies have established the importance of identifying and discriminating zooplankton species, even closely-related and morphologically-cryptic species, to allow understanding of ecosystem function and prediction of impacts of environmental variation and climate change (Johnson et al., 2011; Hare and Kane, 2012; O’Brien et al., 2013), and also to guide fisheries assessment and management (Kelly, 2016; Goodwin et al., 2017; Aylagas et al., 2018). The increasing evidence that COI metabarcoding can provide accurate and reliable species-level identification across the zooplankton assemblage is especially relevant and important for these applications (Andújar et al., 2018).

This study reports the results of comparative molecular (COI metabarcode sequence numbers) and mophological (EcoMon database records for microscopic counts) analysis of six taxonomic groups of marine zooplankton for which both types of data are available ( Supplementary Tables S1 , S2 ). Five of the 6 groups showed significant correlations between COI sequence numbers and morphological counts: Calanoida, Cyclopoida, Eucarida, Gastropoda, and Chaetognatha, but not Hydrozoa ( Figure 3 ). These results provide further evidence of accurate quantitative measurements for some – but not all – taxonomic groups of zooplankton in some – but likely not all – circumstances.

The Shannon Diversity Index (H) based on COI metabarcoding revealed interannual variation for each region, including changes between 2002-2005 and 2007-2011 in the GoM for H index values for 6 taxonomic groups ( Figure 6 ) and 23 species ( Figure 9 ). The changes in some cases were consistent with variation of H index values based on EcoMon morphological counts ( Figure 6 ) and with time-series records of total zooplankton displacement volume in the GoM (O’Brien et al., 2013; Figure 2 ). This finding provides further support for the potential value of COI metabarcoding for revealing and analyzing time-series variation of the zooplankton assemblage and monitoring of ocean ecosystems.

Another positive result, in terms of potential applications of metabarcoding for monitoring of pelagic ecosystems, is the power of COI metabarcoding for detection of species across a number of diverse taxonomic groups of the marine zooplankton assemblage. In this study, a total of 181 species across 23 taxonomic groups were identified based on comparison with the MetaZooGene Database (https://metazoogene.org/MZGdb; Supplementary Table S3 ). These numbers are similar to the numbers of species in EcoMon Survey records, which list 186 zooplankton species across 14 taxonomic groups of metazoan zooplankton detected in NEFSC records since 1977 (NMFS/NEFSC, 2019). However, there is marked lack of overlap in the species detected: only 53 species were found in common between the lists for COI sequences and EcoMon counts. A total of 24 species were detected by both metabarcoding and morphology, considering only species with more frequent observations (including 67 species of 15 groups with > 50 COI sequences and 43 species of 9 groups with >1% frequency of occurrence in EcoMon Survey samples since 1977).

Additional analysis and intercomparison of results from metabarcoding and microscopic counts of the same samples are needed to evaluate and understand both the similarities in total numbers and limited overlap in species identified. One concern is the significant challenges of accurate identification of congeneric and closely-related species of zooplankton based on morphological characters. The discrimination and identification of morphologically cryptic species is a particular power of COI barcoding and metabarcoding (Bucklin et al., 2016; Leray and Knowlton, 2017). In this study, COI metabarcodes detected multiple species of several taxonomically-challenging copepod genera for which species are frequently over-looked or ignored – and rarely counted – in morphological taxonomic analysis, including EcoMon Survey data (NMFS/NEFSC, 2019). Four species of Calanus were detected: C. finmarchicus and C. hyperboreus predominated ( Table 2 ); a few sequences were identified to C. helgolandicus in the 2012 GB sample and to C. glacialis in the 2002 GB and 2005 GoM samples. Seven species of Clausocalanus were detected: C. pergens was most abundant, but C. furcatus, C. jobei, C. lividus, C. mastigophorus, C. parapergens, and C. paululus were detected with small numbers of sequences in several stations; also noteworthy was the detection of four species of Pseudocalanus, including P. acuspes and P. minutus, as well as the cryptic species, P. moultoni and P. newmani ( Supplementary Table S3 ). COI sequence totals across all four Pseudocalanus species were significantly correlated to morphological counts for Pseudocalanus spp. in EcoMon records ( Figure 10 ). Congeneric, sibling, and cryptic species of these copepod genera can be difficult to discriminate morphologically (Hill et al., 2001; Bucklin et al., 2003; Bucklin and Frost, 2009; Crouch et al., 2020), yet species-specific patterns of distribution and abundance are important indicators of seasonal-to-decadal patterns of environmental variation, climate change, and regime shifts (Johnson et al., 2011; Greene et al., 2013; Conversi et al., 2015; Morse et al., 2017).

Further evidence of the usefulness of COI metabarcoding for monitoring of ocean ecosystems is the finding of significant correlations between COI sequence numbers and EcoMon microscopic counts for abundant species for which both types of data are available ( Figure 10 ). The accuracy and reliability of metabarcoding for (semi)quantitative analysis, including abundance or biomass, of zooplankton has been evaluated in previous studies (Elbrecht and Leese, 2015; Thomas et al., 2016; Bucklin et al., 2019). A consistent finding has been that quantitative estimates are more accurate for highly abundant taxonomic groups or species (e.g., Matthews et al., 2021).

An acknowledged limitation of COI metabarcoding of zooplankton diversity is the uneven detection of species across the broad span of taxonomic groups in the pelagic assemblage (Deagle et al., 2014; Clarke et al., 2017; Hajibabaei et al., 2019). Various solutions have been proposed, including using multiple COI sub-regions, with specially-designed primers for target groups (Leray et al., 2013; Corell and Rodríguez-Ezpeleta, 2014; Elbrecht and Leese, 2017; Elbrecht et al., 2019) and integrative multi-region sequence analysis and bioinformatics (Antich et al., 2021; Creedy et al., 2021). Most importantly, classification and identification of species based on COI metabarcodes requires a taxonomically-complete and geographically-appropriate reference sequence database (Leray and Knowlton, 2017; Singh et al., 2021). Continued effort is needed to allow and ensure progress toward inclusion of COI barcode sequences for all zooplankton species, including sibling and cryptic species, recorded from regions throughout the global ocean (Bucklin et al., 2021b).

Marine zooplankton diversity and distribution have been examined using metabarcoding based on a number of different gene regions (Bucklin et al., 2016). A previous study (Bucklin et al., 2019) analyzed the V9 hypervariable region of 18S r RNA for the same set of samples from the NEFSC EcoMon Surveys, providing an opportunity to compare and contrast the results and conclusions based on the two marker gene regions. The 18S r RNA “tree of life” gene occurs in all living organisms on Earth; the gene sequence is both universal and highly conserved, ensuring detection across all major groups of metazoan zooplankton (Amaral-Zettler et al., 2009).

A total of 21 taxonomic groups of zooplankton, ranging from phylum to order, were detected and classified based on V9 metabarcoding of EcoMon samples (Bucklin et al., 2019; Supplementary Table S1 ). Statistical comparison of sequence numbers and abundance counts for these same groups revealed similar patterns of temporal (among years) and spatial (among regions) variation based on 27 samples collected during EcoMon Surveys from 2002-2012. Functional regression analysis for 7 taxonomic groups revealed positive correlations between V9 sequence numbers and abundance counts, with significant correlations (p < 0.05) for 3 groups: Calanoida, Chaetognatha, and Gastropoda. Comparison between the results from V9 and COI metabarcoding analysis of the same set of EcoMon samples demonstrates the power and accuracy of species-level identifications of marine zooplankton.