1.1. Zooplankton community and approach
Zooplankton and ichthyoplankton species are distributed unevenly across marine and aquatic ecosystems, in which their relative species’ compositions, numbers, presence, and survival often are highly sensitive and responsive to environmental conditions [
1,
2,
3]. These species frequently interact together, accentuating or dampening the effects of environmental stressors on one another. Notably, coastal zooplankton and ichthyoplankton communities are regarded as some of the most vulnerable to ocean acidification (OA), hypoxia (H), and warming [
4,
5,
6] – and are subject to other environmental, chemical, physical, and biological fluctuating conditions (e.g., [
7,
8]). Moreover, their larval stages often lack species-level diagnostic morphological characters, precluding their identifications using microscopy and other morphology techniques [
9,
10,
11,
12,
13]. Thus, DNA sequences often provide important means to identify the species in these communities [
11,
12,
13,
14].
In many coastal areas around the world, including the Salish Sea in the temperate Northeastern (NE) Pacific region (
Figure 1), zooplankton and ichthyoplankton diversity levels and the ranges of environmental conditions they encounter are high [
6,
15,
16]. The Salish Sea (which includes the Strait of Georgia, the Strait of Juan de Fuca, and Puget Sound) is one of the most biologically rich and diverse inland seas, housing about 260 fish species and more than 3000 invertebrate species [
15]. These species are adapted to a broad range of physical and chemical conditions, including salinity, temperature, pH, and dissolved oxygen levels that can fluctuate extensively with season, rainfall, tidal influences, land run-off, and other variables [
6,
15,
16]. Because of its relatively low freshwater inputs and deep (>100 m) water column in many areas, the Salish Sea supports a diverse zooplankton and ichthyoplankton community assemblage that is more representative of coastal and oceanic regions than characteristic of most estuaries [
6]. This high zooplankton and ichthyoplankton species diversity distinguishes the Salish Sea from other estuarine ecosystems that have been the focus of most hypoxia and OA studies, making it challenging to fully assess the community diversity using traditional morphological methods. These species diversity levels are relevant because a wide variety of life history strategies, trophic roles, behaviors, and body sizes are represented in these plankton communities, allowing for complex interactions among them, their environment, and their predators. Given that the effects of OAH on trophic interactions are species-specific {5,6], these stressors likely exert a wide range of effects on zooplankton and ichthyoplankton communities.
Moreover, OA conditions in the Salish Sea are commonly at or beyond those regarded as deleterious by the United Nations’ Intergovernmental Panel on Climate Change (IPCC) [
17,
18,
19,
20,
21]. For example, pH levels <7.6 have become common during summer in some areas of Puget Sound, with extremely low values (<7.3) observed in its Hood Canal [
20,
21] (
Figure 1, sites F and G). These extreme conditions, together with their highly regional patterning [
6,
16,
20,
21], provide a unique opportunity to learn about the effects of changing water chemistry and environmental interactions on planktonic organisms. To help address this, the present study places particular focus on copepods and other crustacean, mollusk, and ichthyoplankton species, since field and laboratory investigations have indicated that these taxa often are sensitive to OAH and other perturbations [
16,
17,
18,
19,
20,
21,
22,
23,
24].
The various species comprising zooplankton and ichthyoplankton communities are not ecologically interchangeable, often respond differentially to environmental stressors, and may occur at various times and places [
7,
8,
10,
11,
12,
13]. In other words, these are dynamic communities, upon which other organisms depend for food [
25,
26,
27], and that are critical to ecosystem services and seafood resources [
12,
13,
16]. However, the sheer numbers of samples and the diverse taxonomic coverage needed to identify their component species and monitor such complex community changes are prohibitive using traditional morphological analyses alone (see [
12,
28,
29,
30]).
1.2. Metabarcoding markers evaluated for the zooplankton and ichthyoplankton communities
Our research aims to contribute genetic identification methodology, ground-truthed with zooplankton morphological identifications, towards understanding and monitoring the dynamic responses of zooplankton and ichthyoplankton species and their communities in future studies. Here we employ a high-throughput sequence (HTS) metabarcoding assay approach, which prior studies have indicated can capture a broad diversity of taxa that often are difficult to identify using microscopy and/or other morphological techniques; these include many larval bivalves, gastropods, polychaetes, copepods, and fishes [
28,
29,
30,
31,
32]. Metabarcoding also can detect species that are low in abundance and/or cryptic [
3,
9,
29,
30,
31,
32,
33,
34,
35,
36].
Effective metabarcoding gene regions (markers) require an appropriate level of sequence variability to distinguish among species/taxa of interest, and the primers that flank these diagnostic sequences must anneal to conserved sequences on either end [
31,
32,
33,
34,
35,
36]. When the appropriate gene region/marker and primers are selected, the internal sequence region contains nucleotide variants that can be used to distinguish among the species/taxa in a sample in comparison to reference sequence databases [
28,
29,
30,
31,
32,
33,
34,
35,
36].
Researchers thus aim to select the metabarcoding markers and corresponding primers that are most likely to identify the target taxa/species to the desired taxonomic level [
13,
28,
29,
30,
31,
32,
33,
34,
35,
36]. However, some markers are limited to higher systematic levels because they lack distinctive sequence variation (i.e., the gene region either is too conserved or the marker length is too short to include enough variant nucleotides), and/or available DNA sequence reference repositories are incomplete [
13,
31,
32,
33,
34,
35,
36,
37]. Mitochondrial (mt)DNA gene regions long have been used to identify invertebrate and fish species because their evolutionary lineage tracing is simpler than that of nuclear genes, since mtDNA usually is uniparentally inherited without recombination [
38,
39,
40,
41]. Consequently, many mtDNA gene regions are well-represented in sequence repositories [
3,
9,
11,
13,
28,
29,
30,
31,
32,
33,
34,
35,
36,
41]. Also, the overall abundance of mtDNA in multicellular animals often greatly outnumbers that of their nuclear DNA (summarized by [
42]), which can facilitate detection of rare species and/or smaller-sized individuals in field samples (see [
41]).
One of the more popular mtDNA gene regions for identifying invertebrate and fish species is
cytochrome c oxidase subunit I (COI), which was advocated as the “barcode of life” in the 2000s [
43]. Consequently, many invertebrate and fish species have been sequenced for
COI, yielding high representation in sequence repositories (e.g., [
43,
44,
45,
46,
47,
48,
49,
50]). Our study evaluates its species identifications in the Salish Sea using the metabarcoding
COI marker “LrCOI” designed by Leray
et al. (2013) [
49] as modified by Geller
et al. (2013) [
50]
, which averages 313 base pairs (bp) in length (
Supplementary Table S1), has been widely employed for invertebrates from marine and freshwater habitats [
13,
49,
50,
51,
52], and additionally resolves some fish taxa [
49].
Other frequently used mtDNA sequence regions for species/taxon diagnostics of invertebrates and fishes include the two ribosomal (rRNA) subunit genes,
16S and
12S RNA [
14,
32,
36,
39,
40]. These ribosomal sequence regions overall are more conserved than other mtDNA gene regions, including
COI (above) and
cyt b (below) [
38,
39,
40,
41,
53], since many of their nucleotides are critical to structure and function of the mt ribosomes [
41,
53,
54]. Thus, congeneric species and other relatives can share the same sequences, especially in the short marker regions used in many high-throughput sequencing (HTS) metabarcoding platforms [
36,
51,
52,
53,
54]. This can be alleviated by using longer-read markers; however, the degraded quality of many environmental (e)DNA water samples often precludes the use of longer markers [
37].
Our study evaluates the mt
16S RNA “Cop16S” metabarcoding marker that was designed to identify marine Calanoida Copepoda (Crustacea, Arthropoda) in Tasmania, Australia by Clarke
et al. (2017) and averages about 301 NTs in length [
55] (
Supplementary Table 1). The original study found that Cop16S also resolved other crustacean species, including Cladocera (water fleas), Euphausiacea (krills), and Decapoda (crabs, shrimps, etc.) [
55]. Background data for the present investigation revealed that Cop16S identified copepod and euphausiid species from NE Pacific marine coastal eDNA water samples [
56]. In the freshwater Great Lakes, Cop16S identified copepod species and other (but not all) crustaceans [
31], illustrating the marker’s broad salinity range utility. The present investigation compares species identifications of Salish Sea invertebrates among the LrCOI, Cop16S, and Mol16S (below) metabarcoding markers.
We additionally evaluate invertebrate species identified by a mtDNA
16S RNA metabarcoding marker called “Mol16S”, which was designed to identify freshwater Mollusca to species by Klymus
et al. (2017) [
36] for use in zooplankton and eDNA water samples. Mol16S further resolved other invertebrate species belonging to diverse taxa, including: Annelida (segmented worms), Arthropoda (including Crustacea), Bryozoa (=Ectoprocta; moss animals), Entoprocta (=Kamptozoa), Porifera (sponges), and Rotifera [
31,
36]. Mol16S ranges from 125-185 bp in length (
Supplementary Table S1) and additionally was found to identify some fish species (for which an additional fish “blocker” primer can be added to reduce the fish sequence read signal) (see [
36]). Background data for the present study conducted in marine NE Pacific waters showed that Mol16S identified marine Mollusca, Crustacea, Echinodermata, and some fishes (Teleostei, Actinopterygii, Chordata) species from coastal and deep-sea hydrothermal vent habitats with eDNA water samples [
56]. This metabarcoding marker’s versatility reflects that the shared evolutionary histories of some invertebrate taxa have encompassed biogeographic movements among marine, estuarine, and fresh water habitats [
57], with some species being euryhaline today (see [
58]). Our investigation explores the utility of Mol16S to identify mollusks and other taxa in the Salish Sea.
We also analyze the widely-employed “MiFish12S” mt
12S RNA metabarcoding marker that was designed to identify marine Actinopterygii fish species by Miya
et al. (2015) [
59]. MiFish12S averages about 172 bp in length [
59] and has been broadly used across marine to freshwater ecosystems mostly for eDNA water samples [
59,
60,
61,
62,
63]; here we evaluate it for companion ichthyoplankton studies (
Supplementary Table S1). Our study compares fish species identified by MiFish12S and FishCytb (described next).
Like
COI, the mt
cytochrome (cyt)
b gene evolves relatively quickly, has considerable nucleotide diversity (particularly in the third codon “wobble” position), and a long history of sequencing popularity for fishes and invertebrates, with relatively high representation in repositories (e.g., [
39,
40,
41,
62,
63]). Our study evaluates a “FishCytb” marker that originally was designed to identify freshwater fishes by Snyder and Stepien (2020a) and Stepien et al. (2020) [
62,
63]. FishCytb was targeted to differentiate among all native and nonindigenous fish species in the Great Lakes with eDNA water samples, and averages 154 bp in length) [
62,
63] (
Supplementary Table S1). Preliminary background results also showed its utility to identify marine fish species in NE Pacific coastal eDNA water samples [background data, CAS]. This marker’s versatility traces to the freshwater evolutionary origin of Actinopterygii fishes (summarized by [
64]). Moreover, (a) many fish taxa (including families, some genera, and some species) have moved among marine, estuarine, and/or freshwater habitats during their evolutionary and biogeographic histories [
64], (b) some species are euryhaline, and (c) catadromous and anadromous species transverse salinity realms during their life stages (e.g., some salmonid species, including in the Salish Sea) [
15]. Our investigation compares fish species identifications among metabarcoding markers.