Cryptic Zooplankton Diversity Revealed by a Metagenetic Approach to Monitoring Metazoan Communities in the Coastal Waters of the Okhotsk Sea, Northeastern Hokkaido

Monitoring zooplankton communities is important to understand dynamics in marine ecosystems. However, it is difficult to identify cryptic species and immature stages of zooplankton using morphological classification, which is time-consuming and requires high skill levels. Here, we conducted a metagenetic analysis of the 18S region in 101 zooplankton samples collected weekly throughout 2014 and 2015 at the Okhotsk Tower in Mombetsu, Hokkaido, Japan, and compared the results of this analysis with those provided by morphological analysis. The metagenetic analysis detected 561 molecular taxonomic units (MOTUs), whereas the morphological analysis detected 201 taxonomic groups. Zooplankton communities were dominated by copepods throughout the sampling period; however, non-copepod taxa, which comprised high proportions of both MOTUs (mean 51.1%) and sequence reads (mean 19.1%), were also important. Cryptic diversity detected by the metagenetic analysis was primarily driven by Copepoda and by the larvae of benthic taxa such as Bivalvia, Gastropoda, and Polychaeta. Community structure and diversity varied between periods of warm and cold water, indicating strong correlations with water temperature and thus seasonality. Furthermore, metagenetic analysis revealed detailed seasonal changes in dominant taxa, including larval stages of metazoans with high taxonomic resolutions; these included commercially important organisms such as Japanese scallops. The metagenetic analysis revealed that changes in both water mass and bentho-pelagic interactions sustain ecosystems rich in zooplankton diversity in this area. Metagenetic analysis provides novel insight into zooplankton diversity, and generates massive sequence data that may be used in future research; thus, it is considered an effective tool for monitoring zooplankton communities.


INTRODUCTION
Marine zooplankton comprise an abundant and diverse group including ∼7,000 described species in 15 phyla (Bucklin et al., 2010). Zooplankton contribute substantially to production in higher trophic levels, including in various taxa important to the fishing industry, and zooplankton communities respond sensitively to environmental and climate changes (Hays et al., 2005). For example, negative impacts of global warming on zooplankton abundance have already been reported in the California Current (Roemmich and McGowan, 1995), and changes in zooplankton community structure are predicted as water temperatures increase in the future (Richardson, 2008). Changes in the distribution of ecologically important copepod species have been observed in large-scale spatial and temporal surveys of zooplankton, which in turn affect the recruitment of Atlantic cod in the North Atlantic (Beaugrand et al., 2003). Such changes in distribution are difficult to predict, as zooplankton distribution patterns are species-specific, and different species have different environmental preferences. To mitigate this difficulty, and to track changes in marine ecosystems in response to rapid climate and environmental changes, techniques must be developed that facilitate the monitoring of zooplankton communities at high taxonomic resolutions.
Zooplankton classification is traditionally carried out using morphological observations of individual specimens, a process that is time-consuming and requires great expertise. Furthermore, such methods are imprecise in identifying immature stages and cryptic species, as revealed by the frequency with which undescribed and cryptic species are observed using molecular techniques (Bucklin et al., 2010). Among molecularbased methods, the metagenetic approach, a promising molecular technique using high-throughput sequencing, is particularly useful in revealing community structure based on massive sequence data recovered from bulk environmental samples. Although metagenetic techniques are still being developed in zooplankton, comparisons with morphological techniques already suggest that metagenetic techniques can detect cryptic diversity in zooplankton communities with high taxonomic coverage and at high resolution (Lindeque et al., 2013;Hirai et al., 2015a;Mohrbeck et al., 2015;Harvey et al., 2017). To date, metagenetic techniques have been used to examine zooplankton communities over large spatial scales, and have accurately detected shifts in community compositions corresponding with changes in environmental conditions (Pearman et al., 2014;de Vargas et al., 2015;Hirai and Tsuda, 2015;Chain et al., 2016). Metagenetic techniques have also been applied to revealing zooplankton community structure and diversity at different sampling depths (Pearman and Irigoien, 2015;Sommer et al., 2017) and in different seasons (Casas et al., 2017), as well as to dietary analysis of zooplanktivorous fish . As the utility of metagenetic techniques has been demonstrated for various marine organisms (Aylagas et al., 2016;Danovaro et al., 2016), and as the method is not dependent on subtle morphological characters or the skill of taxonomists, it is expected to become a widespread tool in the rapid and effective monitoring of zooplankton communities and diversity (reviewed in Bucklin et al., 2016).
The northeastern coast of Hokkaido, Japan, provides nursery grounds for juvenile pink and chum salmon, and opens to the Okhotsk Sea Nagata et al., 2007), an area undergoing rapid environmental change as sea ice melts, possibly as a result of air temperatures which have risen over the past 50 years at a rate higher than that of the global average (Ohshima, 2008;Radchenko et al., 2010). An environmental monitoring program has been carried out since the 1990s at the Okhotsk Tower, Mombetsu, where zooplankton communities have been investigated using morphological classifications (brief monthly reports of zooplankton communities are available with images of dominant taxa at http://mombetsu.jp/soshiki/kankou/2013-0902-1336-155.html). Water in the region is dominated by the Soya Warm Current during the summer and is replaced by surface water from the Okhotsk Sea during the winter (Aota, 1975(Aota, , 1985. Drift ice covers the sea surface during the winter, and a peak in chlorophyll a concentration is observed in the spring after the sea ice retreats (Hamasaki et al., 1998;Taguchi et al., 2000). These seasonal environmental changes correspond with shifts of dominant zooplankton taxa (Hamaoka et al., 2010). Given such high seasonality, the northeastern coastal waters off Mombetsu provide ideal conditions for investigating the responses of zooplankton communities to environmental changes. Recently, the metagenetic method was used to identify eukaryotic plankton communities in water samples taken from the region (Nagai et al., 2015) and to identify metazoan zooplankton communities collected by plankton nets (Hirai et al., 2015b). However, only limited numbers of zooplankton samples were analyzed in these previous studies, and the seasonality of zooplankton communities was not fully evaluated.
In this study, we analyzed zooplankton communities collected at the Okhotsk Tower from 2014 to 2015 using a metagenetic analysis of the 18S V9 region, and compared the results of this analysis with those generated by morphological analysis. Furthermore, we used the metagenetic analysis to detect cryptic zooplankton diversity, to investigate seasonal changes in community structure in relation to environmental changes, and thus to identify key factors influencing zooplankton communities in the region. Finally, we discuss the utility of metagenetic analysis as a tool for monitoring zooplankton communities in the coastal waters of northeastern Hokkaido.

Sample Collection
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).

High-Throughput Sequencing
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 4 (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.

Bioinformatics and Data Analysis
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 productmoment 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.

Morphological Analysis
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 3 ) 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.

Seasonal Changes in Environmental Conditions and Zooplankton Abundance
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 3 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 3 to 4.97 × 10 4 individuals/m 3 , peaking in early May after peaks in chl a in both 2014 and 2015, and returning to low levels during the winter.

Zooplankton Diversity across the Sampling Period
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.

Seasonal Changes in Zooplankton Diversity and Community Structure
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 FIGURE 1 | Seasonal changes in environmental conditions and overall zooplankton abundance at Okhotsk Tower, in the coastal waters of northeastern Hokkaido. Water temperature and salinity are provided as means throughout the water column. The concentrations of chlorophyll a, nitrate, and phosphate were measured at the surface. Gray lines in salinity (32.5 and 33.6) indicate the influence of the Okhotsk Sea surface water and Soya Warm Current, respectively.
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.

Seasonal Changes in Zooplankton Taxonomic Compositions
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 FIGURE 2 | The number of major taxonomic groups identified by morphological and metagenetic analyses. Taxa attributable to Copepoda are further classified to order. Taxonomic groups represented by fewer than five morphological groups and molecular taxonomic units (MOTUs) are grouped under "Others." Detailed information for all taxonomic groups is shown in Table S1.
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).

Seasonal Patterns of MOTU Abundance in Major Taxonomic Groups
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)

Cryptic Zooplankton Diversity
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.

Seasonal Changes in Zooplankton Diversity
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.

Seasonal Changes in Dominant Taxa
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, noncopepod 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.

The Metagenetic Approach as a Tool for Monitoring Zooplankton Communities
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(Katakura et al., , 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.

DATA ACCESSIBILITY
Raw Illumina MiSeq data are available from the DDBJ Sequence Read Archive under accession number DRA005786.

AUTHOR CONTRIBUTIONS
JH, SK, and SN: contributed to research design; SK: carried out field sampling; HK: contributed to the measurement of environmental variables; JH: contributed to laboratory work, data analysis, and writing the paper; all authors read and approved the final manuscript.

FUNDING
This study was supported by a research project grant from the Japan Fisheries Research and Education Agency (FRA) and a grant from the Japan Society for the Promotion of Science (grant number 247024).