Zooplankton samples (n = 101) were collected at the Okhotsk Tower, 1 km off the coast of Mombetsu in northeastern Hokkaido, Japan (44°20.2′N, 143°22.9′E), weekly throughout 2014 and 2015. Samples were collected during the daytime using Kitahara nets with a mesh size of 100 μm in a vertical tow from the bottom (~10 m deep) to the surface. Following collection, samples were preserved at −20°C for metagenetic analysis. Samples used in morphological classification (n = 95) were collected using the same method and were preserved in 5% buffered formalin. Vertical profiles of water temperature and salinity were measured using a RINKO-Profiler ASTD102 (JFE Advantech Co., Ltd.) on most days during the sampling period. Because zooplankton samples were collected from the bottom to the surface, mean water temperature and salinity throughout the water column were also measured. Mean values were used because they sufficiently represent seasonal changes in water temperature and salinity, which were the focus of our study. On zooplankton sampling days, surface water was collected using a bucket to measure concentrations of chlorophyll a (chl a), nitrate, and phosphate, according to methods described by Kasai et al. (2017).

In addition to the 101 environmental zooplankton samples described above, and to verify the accuracy of the metagenetic analysis, two mock community samples for which the community structure was known in advance were also prepared. Total genomic DNA from bulk zooplankton samples was extracted using a DNeasy Blood & Tissue Kit (Qiagen). The concentration of DNA in each sample was measured using a Qubit 3.0 Fluorometer (Thermo Fisher Scientific). The 18S V9 region was amplified using PCR with the eukaryotic universal primers 1389F and 1510R (Amaral-Zettler et al., 2009). Three PCRs were carried out using a KOD Plus version 2 (Toyobo). For the first PCR, a 15 μl reaction volume was created, containing 7 μl distilled water, 1.5 μl 10 × buffer, 1.5 μl dNTPs (2 mM), 0.9 μl MgSO₄ (25 mM), 0.9 μl each primer (5 μM), 0.3 μl KOD Plus polymerase, and 2 μl template DNA (1 ng/μl). The first PCR was carried out according to the methods described by Hirai et al. (2017a), with a modification of cycle numbers from 30 to 25. PCR products from the target region were confirmed by electrophoresis using a 2.0% agarose gel. In cases where no clear amplifications were obtained in the first PCR, DNA samples were purified using PowerClean DNA Clean-Up Kit (MO BIO), and the PCR was carried out again. The second and third PCRs, in which adaptors were attached for sequencing in an Illumina MiSeq and dual-index target sequences were used to discriminate samples, respectively, were carried out according to the methods described by Hirai et al. (2017a). The total reaction volume was adjusted to 15 μl containing 2 μl template DNA (1/20 diluted PCR products). The concentration of the final PCR products was measured using a Qubit 3.0 Fluorometer, and each sample was pooled and purified using a QIAquick PCR Purification Kit (Qiagen). Purified PCR products were transferred to Bioengineering Lab. Co., Ltd., where a single sequencing run was carried out using a MiSeq Reagent Kit v2 on an Illumina MiSeq to obtain 2 × 250 bp paired-end sequence reads.

Sequences were filtered for quality and were classified taxonomically according to the methods described by Hirai et al. (2017a). Bioinformatic procedures were carried out using MOTHUR (Schloss et al., 2009). Sequence reads were classified into broad taxonomic groups based on the V9_PR2 reference database (de Vargas et al., 2015) using a naïve Bayesian classifier (Wang et al., 2007). Because this study focused on metazoan zooplankton, only sequence reads identified as “Metazoa” with a threshold >70% were selected for further analysis. Metazoan sequences in the V9_PR2 reference database were aligned using MAFFT (Katoh and Standley, 2013), and were subsequently used as reference sequences for aligning massive sequence reads from the zooplankton samples collected at Okhotsk Tower. Following alignment, preclustering, and the removal of singleton reads, possible chimeras were removed using UCHIME (Edgar et al., 2011), by using both the reference-based mode, the V9_PR2 reference, and the de novo mode. Same numbers of sequence reads were randomly selected based on the minimum number of sequence reads within a single community sample, and massive sequence reads from these samples attributable to Metazoa were clustered into molecular taxonomic units (MOTUs) at a similarity threshold of 97%. In this study, we used MOTUs with an abundance threshold of ≥8 sequence reads. These similarity and abundance threshold values were determined based on the results of the mock community analysis, in order to avoid overestimating diversity by intra-specific sequence variations in 18S RNA genes and erroneous sequence reads produced during PCR and sequencing. MOTUs were grouped based on the results of a naïve Bayesian analysis at a threshold of >70%. Unclassified MOTUs were further investigated using a BLAST search of the NCBI database, and were assigned to taxonomic groups in cases where ≥95% sequence identity was observed. Additionally, a BLAST search was carried out for the 10 most abundant MOTUs in the major taxonomic groups of Bivalvia, Copepoda, Gastropoda, and Polychaeta. These abundant MOTUs were assigned either to the family (Bivalvia, Copepoda, and Polychaeta) or superfamily (Gastropoda) level, in cases where they were ≥99% similar to taxonomic groups in the public database or to higher taxonomic levels when the similarity threshold was ≥95%. Although there are possible methodological biases, including copy numbers of 18S rRNA genes and primer mismatches during PCR, sequence read abundance in metagenetic analyses can be used as a proxy for biomass in zooplankton (Lindeque et al., 2013; Hirai et al., 2015a). Therefore, we used proportions of sequence reads for quantitative analysis in the metagenetic analysis in this study.

Zooplankton diversity was measured using the total number of MOTUs in conjunction with Shannon and Simpson diversity indexes. Correlations between environmental variables and MOTU numbers were determined using Pearson's product-moment correlation coefficient (r). Bray–Curtis similarity was calculated between community samples for multidimensional scaling (MDS) analysis in PRIMER version 7 (Clarke and Gorley, 2015), based on the abundance of sequence reads and the presence/absence of MOTUs in each sample. In the abundance-based analyses, the proportions of sequence reads attributable to each MOTU were log_x+1-transformed. The strongest environmental determinants of zooplankton communities were analyzed in a redundancy analysis (RDA) in CANOCO v.5, using the forward selection method (ter Braak and Verdonschot, 1995). Seasonal changes in taxonomic composition were also investigated, based on the presence/absence of specific MOTUs or sequence reads.

Specimens were classified to the lowest-ranking taxon possible using morphological techniques, and the abundance of morphological groups in each sample was counted. These counts, as well as diversity indexes, MDS plots, and taxonomic compositions, were compared with those generated in the metagenetic analysis. Because this study did not collect biomass data; instead, the abundance of each taxonomic group (individuals/m³) was used in the quantitative analyses used to measure diversity indexes and community similarity. Correlations between diversity measured metagenetically and diversity measured morphologically were compared using Pearson's product-moment correlation coefficient (r). In the abundance-based community analysis, similarities between communities were calculated using log_x+1-transformed abundance data from each morphological group.

Mean water temperatures ranged from −1.6 to 20.2°C in 2014 and from −1.6 to 20.3°C in 2015 (Figure 1). Water temperature peaked from July to September, and temperatures < 0°C were primarily observed from January to March. Salinity levels were low in the winter and high in the summer, although occasional abrupt decreases in salinity were observed throughout both years. Low salinity levels (≤32.5) were observed primarily during November and May. Salinity levels ≥33.6 were primarily detected from July to October. Chl a concentration ranged from 0.1 to 13.3 mg/m³ during the sampling period, and peaked in spring in both 2014 and 2015; however, high chl a concentrations were occasionally observed in summer and autumn, particularly in 2015. Chl a concentrations remained low during winter. Concentrations of nitrate and phosphate were highest during winter, with peaks of 18.3 and 1.5 μM, respectively, and began decreasing as chl a concentrations increased in early in spring (March–April). Concentrations of these nutrients remained low during summer and autumn. Zooplankton abundance ranged from 2.63 × 10³ to 4.97 × 10⁴ individuals/m³, peaking in early May after peaks in chl a in both 2014 and 2015, and returning to low levels during the winter.

A minimum number of 15,939 sequence reads remained within a single community sample after quality-filtering of raw sequence reads. For standardization, 15,939 sequence reads were therefore randomly selected from each of the 101 community samples, and so a total of 1,609,839 reads from the samples were used in the following analyses. After clustering at the 97% similarity threshold and applying the abundance threshold, 561 MOTUs were identified. The number of MOTUs identified across the sampling period was higher than the number of morphological groups identified (201), and the number of taxa attributable to each taxonomic group was generally higher in the metagenetic analysis than in the morphological analysis (Figure 2).

Among the taxa identified, the highest diversity was observed in Copepoda, which accounted for 189 MOTUs and 101 morphological groups. Within Copepoda, the majority of diversity was driven by Calanoida (77 MOTUs and 33 morphological groups) and Cyclopoida-Poecilostomatoida (35 MOTUs and 25 morphological groups). Further diversity was contributed by Harpacticoida, the only group in which more morphological groups (42) were detected than MOTUs (32). A further 40 MOTUs attributable to Copepoda could not be classified to lower taxa. In addition to Copepoda, the high diversity detected in the metagenetic analysis was largely driven by meroplankton diversity. Among this group, Polychaeta in Annelida were particularly diverse, contributing 78 MOTUs, although morphological classification identified only 19 groups. The same pattern was observed in various other groups, where the numbers of MOTUs were generally higher (24 in Bivalvia, 9 in Cirripedia, 29 in Gastropoda, and 25 in Platyhelminthes) than the number of morphological groups identified (1 in Bivalvia, 2 in Cirripedia, 5 in Gastropoda, and 1 in Platyhelminthes). In some cases, MOTUs were identified in groups not detected by morphological analysis (7 in Brachiopoda, 7 in Bryozoa, 21 in Nematoda, and 5 in Polyplacophora). Gelatinous plankton attributable to Hydrozoa contributed 21 MOTUs and 12 morphological groups. Fifty MOTUs remained unclassified at any rank lower than Metazoa. Minor taxonomic groups (<5 morphological groups or MOTUs; these comprised a total of 13 morphological groups and 21 MOTUs) are grouped under “Others” in Figure 2, and detailed information about the diversity of all taxonomic groups is shown in Table S1.

The number of MOTUs identified per sample ranged from 29 to 133, and was always higher than the number of morphological groups identified, which ranged from 5 to 55 (Figure 3). Seasonal changes in the numbers of MOTUs and morphological groups were correlated (r = 0.64, p < 0.01). Water temperature was positively correlated with the number of MOTUs (r = 0.71, p < 0.01) and the number of morphological groups identified (r = 0.64, p < 0.01), resulting in large numbers of taxa from summer to autumn and low numbers of taxa during the winter. Similar seasonality was detected in the Shannon and Simpson diversity indexes, although the seasonality was not as pronounced in the Simpson index. Significant correlations between the metagenetic and morphological analyses were observed both in the Shannon and Simpson diversity indexes (p < 0.01); however, the correlation was stronger in the Shannon index (r = 0.61) than in the Simpson index (r = 0.32).

In addition to seasonal changes in diversity, MDS ordination revealed clear seasonal changes in community structure in both metagenetic and morphological analyses (Figure 4). Water temperature was the strongest environmental factor influencing both abundance-based (sequence reads or numbers of individuals) and non-abundance-based (presence/absence of MOTUs or morphological species) measurements (p < 0.01 for all analyses). Seasonal change was clockwise, characterized by gradual changes in zooplankton communities correlated with increases or decreases in water temperature. The same seasonal pattern was observed in both 2014 and 2015.

Copepod MOTUs contributed substantially to zooplankton diversity in every month throughout the sampling period (ranging from 27.4 to 79.5%; mean 48.9%; Figure 5A). Proportions were particularly high during the winter and low during the summer and autumn. Among copepods, the highest proportion of MOTUs was contributed by Calanoida (average 27.5%), followed by Cyclopoida-Poecilostomatoida (mean 8.9%) and Harpacticoida (mean 7.3%); the mean proportion of unclassified Copepoda was 5.2%. In the morphological analysis, copepods comprised 23.5–100% of samples (mean 58.0%; Figure 6A). Calanoida contributed most to copepod diversity in the morphological analysis (mean 28.5%), followed by Harpacticoida (mean 15.3%) and Cyclopoida-Poecilostomatoida (mean 14.1%).

Among non-copepod taxa, high proportions of MOTUs and morphological groups were attributable to meroplankton, including Bivalvia (MOTU maximum 14.3%, mean 4.7%; morphological maximum 11.1%, mean 2.4%), Cirripedia (MOTU maximum 7.5%, mean 3.8%; morphological maximum 20.0%, mean 3.9%), Gastropoda (MOTU maximum 12.7%, mean 6.0%; morphological maximum 15.4%, mean 4.4%), and Polychaeta (MOTU maximum 25.0%, mean 12.1%; morphological maximum 29.4%, mean 13.4%). Rare taxa were proportionally more common in the metagenetic analysis (mean 12.5%) than in the morphological analysis (mean 5.8%). Unclassified MOTUs comprised 7.0% of taxa in the metagenetic analysis.

Copepods were particularly abundant among the sequence reads, accounting for 17.7–99.6% of all reads per month (mean 80.9%; Figure 5B). Among copepods, Calanoida were the most abundant (14.3–96.0%; mean 55.6%), particularly in winter. Cyclopoida-Poecilostomatoida were also common (1.6–69.3%; mean 21.9%), and were sometimes more abundant than Calanoida. Harpacticoida and unclassified Copepoda sequence reads were comparatively uncommon (mean 2.9 and 0.5%, respectively). However, high proportions of Harpacticoida (up to 37.7%) were occasionally observed. Copepods were also abundant in the morphological analysis, accounting for 31.7–100% (mean 77.7%) of individuals (Figure 6B). However, high proportions of Copepoda individuals sampled were unclassified nauplii (mean 30.5%). Among the classified taxa, the most dominant taxa belonged to Calanoida (mean 27.7%), followed by Cyclopoida-Poecilostomatoida (mean 16.6%). Harpacticoida were comparatively uncommon (mean 2.9%).

Besides copepods, meroplankton were also comparatively abundant in sequence reads (Figure 5B). Cirripedia abundance peaked in early spring (maximum: 42.0%). Peaks in Bivalvia and Gastropoda abundance of 14.2 and 20.9%, respectively, were also observed, primarily in summer. Polychaeta abundance, which was observed primarily in spring and autumn, averaged 3.9% and peaked at 28.5%. In addition to meroplankton, Hydrozoa occasionally contributed high proportions of sequence reads (maximum 66.9%, during May and June 2014); however, proportions were generally low. High proportions of Euphausiacea were observed in April in both 2014 and 2015 (maximum 22.6%). In the morphological analysis, high proportions of meroplankton were observed coinciding with those detected by metagenetic analysis, although the maximum proportions differed (60.2% in Cirripedia, 32.1% in Bivalvia, 62.4% in Gastropoda, 53.1% in Polychaeta, 17.7% in Euphausiacea, and 4.2% in Hydrozoa; Figure 6B). Unclassified Metazoa and sequences classified as “Others” were not as abundant in the sequence read analysis (mean 2.7% for “Others,” 0.8% for unclassified Metazoa), and similar results were found in the morphological analysis, where the proportion of “Others” was 0.8%. Peaks in “Others” were observed both in sequence read abundance and in morphological abundance during the summer (peaking at 11.3 and 7.4%, respectively).

The monthly means of sequence reads and abundance showed similar peaks between metagenetic and morphological analyses of major taxa of Bivalvia, Copepoda, Gastropoda, and Polychaeta (Figure 7). In the morphological analysis, the major morphological groups were the larval stages of Bivalvia, Gastropoda, and Polychaeta, as well as Copepoda nauplii, that were not classified to low ranks. However, metagenetic analysis revealed seasonal changes in the abundance of MOTUs classifiable to various taxonomic groups (Figure 7), including shifts in dominant taxa on a monthly basis. In Bivalvia, the dominant MOTUs in August 2014 belonged to Mytilidae and Veneridae, whereas in October 2014 they belonged to Corbiculidae (Figure 7A). The dominant Bivalvia MOTUs in June 2015 and September 2015 were attributable to Mactridae and Mytilidae, respectively. One MOTU classified to Pectinidae (MOTU 109) was relatively abundant in June 2014 and May 2015. The morphological analysis classified adult Copepoda into species, and both MOTUs and morphological species showed seasonal shifts in dominant taxa (Figure 7B). Correspondence between metagenetic and morphological analyses was observed in Oithonidae (MOTU 1; including Oithona similis), Paracalanidae (MOTU 3; including Paracalanus parvus s. l.), Clausocalanidae (MOTU 6; including Pseudocalanus newmani or P. minutus), and Centropodoidea (MOTUs 2, 4, 9, and 11; including the genus Acartia and Eurytemora). There were, on the other hand, major MOTUs that were only detected in the metagenetic analysis, which were attributed to the families Calanidae (MOTU 5), Pontellidae (MOTU 8), and Metridinidae (MOTU 10). Gastropoda MOTUs were classified to superfamily-level, due to the relatively low sequence variations among taxonomic families. However, MOTU-based analysis successfully detected seasonal changes in dominant Gastropoda, compared to morphological analysis (Figure 7C). The dominant MOTUs Cerithioidea (MOTU 37), Buccinoidea (MOTU 38), and Cephalaspidea (MOTU 43) were replaced by Truncatelloidea (MOTU 29), Caenogastropoda (MOTU 41), and Bulloidea (MOTU 51) during May and August. Polychaeta larval stages were partially classified to family-level by morphological analysis, and metagenetic and morphological analyses showed the same peaks of some dominant taxa, including Spionidae and Sabellariidae (Figure 7D). Seasonal occurrences in dominant Polychaeta were clearly shown by the MOTU-based analysis. Although detailed monthly changes in most of major taxonomic groups were different when comparing 2014 and 2015, Polynoidae (MOTU 14) was dominant from March to July in both years. Spionidae (MOTU 19) showed a high proportion of sequence reads in November and December in 2014, and other dominant MOTUs also showed specific peaks of sequence reads, such as Nephtyidae (MOTU 35) peak in July in both 2014 and 2015.

Metagenetic analysis revealed a higher diversity of zooplankton taxa than those revealed by morphological analysis. The 18S V9 region used in the present study is highly variable; however, species-level data are not sufficiently high-resolution in some taxonomic groups of zooplankton, making accurate classification difficult in certain cases (Wu et al., 2015). Although a high similarity threshold (e.g., 99%) is expected to be used when analyzing metazoan zooplankton (Mohrbeck et al., 2015), we instead used a comparatively low similarity threshold (97%), a technique recommended to prevent overestimations of diversity (Kunin et al., 2010) that is often used in metagenetic analyses of zooplankton using the 18S V9 region (Pearman et al., 2014; Casas et al., 2017). Additionally, we used a strict abundance threshold (eight minimum sequence reads), because rare MOTUs are not always reliable and lead to overestimations of diversity (Bokulich et al., 2013; Hirai et al., 2017b; Leray and Knowlton, 2017; Sommer et al., 2017). Indeed, an even stricter minimum threshold of 16 sequence reads would have identified 429 MOTUs, a number much higher than that detected morphologically. Metagenetic analysis is not free from errors which sometimes lead to the overestimation of MOTU numbers; however, the strict criteria used in the present study nevertheless indicate a richer, often more cryptic diversity than previously recognized based on morphological observations in the study area.

Much of this cryptic diversity was contributed by Copepoda, particularly Calanoida, because molecular techniques are capable of classifying even nauplii and cryptic species (Bucklin et al., 2010). However, a substantial proportion was contributed by meroplankton taxa. Meroplankton taxa are primarily the larvae of benthic organisms, which are difficult to identify morphologically, leading to the increased diversity detected by metagenetic analysis. Although other metagenetic studies targeting the 18S region have reported cryptic diversity among meroplankton (Lindeque et al., 2013; Pearman et al., 2014; Sommer et al., 2017), the contribution of meroplankton to overall diversity was particularly high in the present study. For example, as many Polychaeta MOTUs were detected as calanoid MOTUs. This unexpectedly high diversity in meroplankton may be attributable to shallow water depth in the study area, suggesting that the metagenetic method is especially effective in assessing zooplankton diversity at this location. Metagenetic analysis also detected higher diversity in various other taxonomic groups, including gelatinous organisms such as Hydrozoa, providing further evidence for the existence of cryptic zooplankton diversity in the study area, and suggesting that morphological analysis alone is not sufficient to characterize plankton communities accurately.

The seasonal changes in temperature, salinity, chl a, and nutrient content measured in the present study conform with patterns measured from 1996 to 2014 in the study area (Kasai et al., 2017), and changes in zooplankton abundance followed patterns reported previously (Hamaoka et al., 2010). Although we observed abrupt changes in environmental variables, such as declines in salinity and increases in chl a, possibly caused by increased rainfall, influx of inland water, or turbulence, such changes have been reported as characteristic of study sites located in coastal regions (Hamasaki et al., 1998; Kasai et al., 2017). Thus, zooplankton samples used in this study were collected under typical environmental conditions and can be considered representative of common seasonal changes of zooplankton communities in the area.

The substantial variations in environmental conditions typical of the study area drove seasonal changes in zooplankton diversity. Patterns of presence and absence in different taxa differed clearly in winter and summer in both metagenetic and morphological analyses in correlation with water temperature. This result conforms to those reported in a previous morphological study that showed seasonal variability in species distribution was driven by differing water masses (Hamaoka et al., 2010). Because high zooplankton diversity, such as in copepods, is primarily observed in warm-water masses (Rombouts et al., 2009; Tittensor et al., 2010), water exchanges during summer and autumn may have generated the high diversity detected in the present study. Additionally, proportionally more non-copepod taxa, particularly among meroplankton, were detected by metagenetic analysis than by morphological analysis throughout the year, suggesting that meroplankton contribute substantially more to diversity than previously recognized, particularly from summer to autumn. Thus, metagenetic analysis revealed clear seasonality in zooplankton diversity at high taxonomic resolution, and environmental correlations suggest that these changes in diversity were driven by water mass exchanges and interactions between planktonic and benthic communities.

Although there are methodological limitations (e.g., PCR bias, gene copy numbers), sequence reads provide semi-quantitative information for zooplankton biomass (Lindeque et al., 2013; Hirai et al., 2015a; Albaina et al., 2016; Clarke et al., 2017), which is not directly comparable with the number of individuals in morphological analysis. We nevertheless observed similar seasonal patterns in diversity indexes, abundance-based community structures, and taxonomic compositions between the metagenetic and morphological analyses. These results support the utility of abundance-based analysis of metagenetic data in detecting dominant taxa. The relatively low correlation in Simpson index and some discrepancies in dominant taxa we measured may have been caused by the use of abundance rather than biomass data in the morphological analysis. Furthermore, fluctuations in diversity patterns, including significant dips in diversity, particularly in the Simpson index, may have been caused by occasional changes in environmental variables, or by the presence of large-bodied taxa (e.g., high proportions of Chaetognatha and Hydrozoa in 2014).

The metagenetic analysis provided detailed taxonomic data using sequence reads, whereas morphological classification was often unsuccessful. Our MOTU-level analysis of Bivalvia, Gastropoda, and Polychaeta revealed seasonal and annual changes in the occurrence and dominance of certain taxa that were not detected in the morphological analysis, due to the difficulty of identifying meroplankton. Additionally, copepod nauplii were particularly abundant in the morphological analysis; however, most sequence reads of copepods were successfully classified into detailed taxonomic groups throughout the sampling period in the metagenetic analysis. Although the majority of sequence reads were attributable to copepods, non-copepod taxa made non-trivial contributions, particularly from summer to autumn. Indeed, high diversity index values were observed in warm periods, indicating the co-existence of multiple species from various groups, which may promote zooplankton productivity in the study area. The proportional contribution of meroplankton to sequence read counts was greater than that detected in a metagenetic study in Tosa Bay, off Japan, (Hirai et al., 2017a) or in reefs in the Red Sea (Pearman et al., 2014). Furthermore, stable isotope analyses in the coastal waters of northeastern Hokkaido indicated an interaction between planktonic and benthic organisms (Hiwatari et al., 2011). As observed for zooplankton diversity, environmental factors such as water temperature, driven by water mass exchange, and a strong bentho-pelagic interaction in coastal waters, acted together to influence seasonal changes in dominant taxa.

Metagenetic approaches provide a sensitive method for monitoring zooplankton diversity, as demonstrated in studies of zooplankton communities in the English Channel (Lindeque et al., 2013), tropical and subtropical Pacific (Hirai and Tsuda, 2015; Hirai et al., 2015a), station ALOHA (Sommer et al., 2017), Red Sea (Pearman et al., 2014; Pearman and Irigoien, 2015; Casas et al., 2017), Monterey Bay (Harvey et al., 2017), Canadian coastal areas (Chain et al., 2016), and Storm Bay (Clarke et al., 2017). Given the alignment of results generated from metagenetic and morphological methods, we conclude that metagenetic analysis was successfully used to monitor changes in zooplankton diversity and communities in the study area, and that this in turn corresponded with environmental changes, particularly water temperature. Metagenetic analysis yielded higher taxonomic coverage and resolution, and the strong bentho-pelagic interactions and importance of meroplankton in this study meant that the metagenetic technique was particularly effective in evaluating zooplankton diversity in the study area. Furthermore, because the MOTU-level analysis of Bivalvia detected a Pectinidae MOTU, possibly derived from Japanese scallop (Mizuhopecten yessoensis), we also posit that the technique is useful in monitoring the dynamics of specific taxa, including commercially important species in the study area. In a previous study, immunostaining was used to detect M. yessoensis larvae in the study area (Katakura et al., 2013, 2015); however, metagenetic analysis functions as a taxonomically comprehensive method of classification that is not restricted to a specific species. Conversely, metagenetic analysis provides only semi-quantitative data (e.g., proportions of sequence reads), and thus, does not provide detailed information concerning the size, number, or developmental stage of individuals detected (Lindeque et al., 2013). Thus, a combined analysis, using both morphological and metagenetic techniques, is the optimal approach for monitoring zooplankton communities.

Relatively large numbers of sequences have been registered in 18S rRNA gene (Quast et al., 2013; de Vargas et al., 2015); however, we encountered many unclassified metazoan MOTUs. To further develop metagenetic analysis techniques, we recommend additional efforts to accumulate even more zooplankton sequences than exist in current databases. Furthermore, we used the 18S V9 region because of its capacity to identify various metazoan zooplankton; however, marker choice affects the results of the metagenetic analysis. For example, although mitochondrial genes such as COI have high species-level resolution, primer mismatches remain common, particularly in highly variable regions of genes (Deagle et al., 2014; Clarke et al., 2017). Thus, molecular markers should be chosen according to the purpose of each study.

The coastal waters of the Okhotsk Sea provide both a productive fishery ground and an important locality for monitoring marine ecosystems (Radchenko et al., 2010). Because climate and environmental changes are proceeding rapidly, it is important to understand existing zooplankton community structures when predicting future changes. However, the present study generated results over a span of only 2 years; accordingly, additional longer-term zooplankton monitoring efforts are recommended to further characterize marine ecosystems in the region and to identify the effects of long-term climate changes on these ecosystems. Nevertheless, because the data generated in this study do not depend on subtle morphological differences or on taxonomic expertise, they will provide a reliable and straightforward dataset for implementation in future studies using different spatial and temporal spans. The present study showed that metagenetic analysis successfully detected cryptic diversity in the study area, with implications for detailed taxonomic information concerning dominant zooplankton taxa. Thus, we consider metagenetic analysis is an especially useful tool for monitoring metazoan zooplankton in the coastal waters of the Okhotsk Sea, northeastern Hokkaido.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.