Metabarcoding Analysis of Ichthyoplankton in the East/Japan Sea Using the Novel Fish-Specific Universal Primer Set

The spatiotemporal distribution of fish larvae and eggs is fundamental for their reproduction and recruitment in aquatic ecosystems. Here, a metabarcoding strategy was employed as an alternative to a conventional ichthyoplankton survey, which requires a considerable amount of time, labor, and cost. First, a piscine-specific universal primer set (FishU) was designed to amplify the region, flanking the highly conserved mitochondrial 12S and 16S ribosomal genes, and it was optimized for the MiSeq platform. Based on both in silico and in vitro analyses, the newly designed FishU primers outperformed the two previously reported fish-specific universal primer sets (ecoPrimer and MiFish) in taxon coverage, specificity, and accuracy in species identification. The metabarcoding results by FishU primers successfully presented the diversity of ichthyoplankton directly from the zooplankton net samples in the East/Japan Sea, presenting more accurate and plentiful species numbers than those by MiFish primers. Thus, the metabarcoding analysis of ichtyoplankton using the newly designed FishU primers is a promising tool for obtaining useful data to understand the reproduction of fish, such as spawning sites, reproductive periods, population structures, feeding ecology, and diet.


INTRODUCTION
DNA metabarcoding is a recently established technique that makes it possible to conduct taxonomic identification of entire assemblages in environmental samples via high-throughput sequencing (Yu et al., 2012). Because of its reliable results with relatively low cost and labor, this technique has become one of the most widely used methods in ecological studies, including biodiversity assessment (Valentini et al., 2016;Bylemans et al., 2018), the detection of invasive species (Borrell et al., 2017;Klymus et al., 2017), and feeding ecology (Iwanowicz et al., 2016;Yoon et al., 2017;Hawlitschek et al., 2018). The environmental DNA (eDNA) derived from organisms in water samples can also be directly analyzed by eDNA metabarcoding, another promising method that has little impact on the ecosystem during sample collection (Thomsen et al., 2012;Sigsgaard et al., 2017;Yamamoto et al., 2017;Bylemans et al., 2019;McDevitt et al., 2019). As a result, metabarcoding is now being considered by many researchers as a reliable alternative to replace time-consuming and laborious traditional survey methods (Berry et al., 2015;Aylagas et al., 2016;Shaw et al., 2016;Watts et al., 2019).
Most metabarcoding analyses are currently based on the assemblage of PCR products amplified by a universal primer for specific taxa. Several fish-specific universal primers have been developed targeting 12S rRNA regions, such as EcoPrimers (Riaz et al., 2011) and MiFish (Miya et al., 2015), 16S rRNA region primers (Evans et al., 2016;Shaw et al., 2016), cytochrome b (Minamoto et al., 2012;Thomsen et al., 2012) and cytochrome c oxidase I (Balasingham et al., 2018;Collins et al., 2019). Among them, the primers targeting the 12S rRNA and 16S rRNA regions generally showed higher performance in the broader taxonomic range and species-level assignment than the primers targeting the cytochrome b region (Zhang et al., 2020). The MiFish primer set demonstrated reliability for fish biodiversity analysis of eDNA samples for both seawater (Ushio et al., 2017;Yamamoto et al., 2017) and freshwater (Sato et al., 2018). The bioinformatics tool, web-based MiFish pipeline 1 also provides a standardized and convenient bioinformatic platform, which has helped researchers obtain reliable metabarcoding data without considering the complicated bioinformatics processes from the raw reads generated by the sequencer (Sato et al., 2018). Therefore, the MiFish platform is now being widely used for the study of fish biodiversity using eDNA metabarcoding analysis (Andruszkiewicz et al., 2017;Yamamoto et al., 2017;Bylemans et al., 2019). In particular, owing to its small amplicon size (∼170 bp), the MiFish platform is currently considered one of most reliable tools for the eDNA metabarcoding analysis of fish taxa. However, the small amplicon size of MiFish is a double-faceted property. Although it may increase the chance of detection rate, its small amplicon size may not have enough sequence to discriminate a species from the closely related ones (Bylemans et al., 2018;Zhang et al., 2020). In fact, MiFish sequences could not discriminate species in several genera, including Sebastes and Takifugu (Yamamoto et al., 2017). In addition, the short MiFish sequence may require sufficient reference sequences in the database for the local fish assemblage, which would be another barrier to adopting the MiFish platform directly for local analysis. Even a single nucleotide difference in the region may result in the assignment of a wrong species in the MiFish platform (Yamamoto et al., 2017). Therefore, it was necessary to design a primer set covering larger barcode sequences for higher discrimination power.
Ichthyoplankton refers to the eggs and larvae of fish, which is usually found at depths of less than 200 m, in what is known as the epipelagic zone. Their surveys have been conducted for a long time in the East/Japan Sea to obtain information about fish resources, ranging from the spawning areas and seasons to the estimated numbers of spawning stocks and changes in the distribution of certain species (Keller et al., 1999;Shoji et al., 2011). Due to the simple coastlines, fisheries in the East/Japan 1 http://mitofish.aori.u-tokyo.ac.jp/MiFish/ Sea have been strongly affected by the collision of two main currents along with the Korean peninsula: the North Korea cold current (NKCC) from the North and East Korea warm current (EKWC), a branch of the Tsushima current (Kim et al., 2007b). In addition, the transport processes of eggs and larvae are largely dependent on the variability in these two currents, affecting the recruitment and productivity of each fish stock (Yatsu et al., 2005). Therefore, ichthyoplankton surveys have been one of the main methods to understand fisheries in the East/Japan Sea. However, traditional ichthyoplankton surveys depend mainly on microscopic morphological observations, which require a considerable amount of time and effort to sort and count the eggs and larval fish from the zooplankton net samples. In addition, morphological identification of several fish larvae and eggs remains a challenge even for experienced specialists. These drawbacks have been a major bottleneck in conducting largescale ichthyoplankton surveys.
Therefore, there is a need for the use of metabarcoding analysis for ichthyoplankton surveys. However, the use of short barcodes, such as MiFish, often results in the misidentification of exotic species due to ichthyoplankton being transported from distant locations by currents. To avoid potential issues with identification, an accurate species analysis using longer barcodes was needed. In this study, we designed novel fishspecific universal primers (FishU) to analyze ichthyoplankton from the zooplankton net samples. FishU yields longer amplicons that have a higher taxon specificity than MiFish. Their reliability was compared with previously designed fish-specific universal primers, including ecoPrimer and MiFish. We also compared the metabarcoding data generated by two primers, MiFish and FishU, using zooplankton samples obtained from the East/Japan Sea. These data show the FishU primer set provides a higher accuracy in species assignment and recovers more species than the MiFish primer set. Thus, the newly developed FishU primers represent a useful tool for ichthyoplankton surveys in countries with relatively limited reference sequence information, facilitating a molecular strategy in fishery surveys.

Design of Universal Primer Set for Fish Taxa
A pair of universal fish primer sets (FishU) was designed to meet three goals: (i) wide taxon coverage to present most of the fish species, (ii) a high degree of taxon specificity to amplify fish DNA from mixed zooplankton samples, and (iii) an optimized amplicon size for the Illumina MiSeq platform (300-500 bp). FishU primers were designed to amplify the region between the end of 12S and 16S based on the multiple alignment of 1,808 fish mitogenome sequences obtained from the MITOFISH database (Iwasaki et al., 2013) (Supplementary  data 1). The expected amplicon size by FishU was approximately 373 bp (Figure 1).
The FishU primers were validated by in silico analysis using 1,808 fish mitogenome sequences in the database. The OBITools software (Boyer et al., 2016) was used to calculate the amplified species numbers according to the mismatched nucleotide numbers of the three universal primers (ecoPrimer, MiFish, and FishU). Taxonomic resolutions amplified by three primers were calculated by in silico PCR using the ecotaxspecificity script of OBITools (Boyer et al., 2016).

Examination of the Reliability of the FishU Primer Set
The taxon specificity of the FishU primer set was examined by PCR analysis (Figure 2). The amplification by FishU was tested with 51 individual fish samples collected from Korean waters, which covered 50 genera, 37 families, and 14 orders.
Specificity for fish species was tested by PCR with the most abundant invertebrate taxa in the zooplankton net samples, including Metridia pacifica (copepod), Euphausia pacifica (krill), Sagitta elegans (arrow worm), and Penaeus monodon (shrimp). Genomic DNA was extracted using the AccuPrep R Genomic DNA Extraction Kit (Bioneer, Republic of Korea) following the manufacturer's instructions. PCR amplification was performed with extracted DNA as templates and three different primer sets (ecoPrimer, MiFish, and FishU). PCR amplification with the FishU primer set was conducted using the following cycling conditions: initial denaturation at 94 • C for 3 min, followed by 35 cycles of 94 • C for 30 s, 48 • C for 30 s, and 72 • C for 30 s with a final extension at 72 • C for 3 min.
PCR conditions with MiFish and ecoPrimer sets were the same as above except for the annealing step at 50 • C for 10 s. The PCR reaction mixture (20 µL) contained 100 ng template, 1 µL of each primer (10 pmol), 2 µL of dNTPs (10 mM), 0.2 µL exTaq HS (Takara, Japan), and 2 µL of 10 × buffer and distilled water to a final volume of 20 µL. The amplification products were separated by 1.5% agarose gel electrophoresis and stained with loading star dye (Dynebio, Republic of Korea).
The taxon coverage and specificity of the three primers was also compared through metabarcoding analysis. The library for MiSeq sequencing was constructed using the Nextera XT index kit (Illumina, United States). Genomic DNA (200 ng) in two random zooplankton samples was amplified using three primer sets with the adapter. The first PCR conditions of FishU involved initial denaturation at 94 • C for 3 min, followed by 30 cycles of 94 • C for 30 s, 48 • C for 30 s, and 72 • C for 30 s with a final extension at 72 • C for 3 min. The first PCR amplification was performed with a 20 µL reaction volume containing 200 ng of template, 1 µL of each FishU with adapter primer (20 pmol), 0.5 µL of dNTPs (10 mM), 0.2 µL of ex hot start Taq (TaKaRa, Japan), and 2 µL of 10 × exTaq buffer. The PCR reaction mixture and conditions with the MiFish and ecoPrimer set were identical to the FishU mixture except for the primer with adapter (each 5 pmol). The first amplicons were visualized by 1.5% agarose gel electrophoresis, and the fragment of expected size was collected. The collected fragments of each primer set were purified using the AccuPrep R Gel Purification Kit (Bioneer, Republic of Korea) and eluted with 20 µL of elution buffer. The purified products were used to construct a library using the Nextera XT index kit (Illumina, United States) according to the manufacturer's instructions. Nextera libraries were purified using the AccuPrep R Gel Purification Kit (Bioneer, Republic of Korea). After the quality and quantity of the library were measured using the 2100 Bioanalyzer (Agilent Technologies, United States) and qubit dsDNA HS Assay Kit (Invitrogen, United States), sequencing was performed using the Illumina MiSeq (2 × 300 bp pair ends).

Metabarcoding Analysis of Ichthyoplankton From Zooplankton Net Sample
Zooplankton net samples were collected from 12 sample stations from the coastal waters along the Korean Peninsula as part of a project funded by the Ministry of Oceans and Fisheries, Korea (Figure 3, Supplementary Table 1). Plankton samples were collected by oblique tows using a bongo net (60 cm in diameter, 330 µm in mesh size) monthly from January to June 2016. The collected zooplankton samples were immediately stored in five volumes of 99.5% ethanol (Daejung Chemicals & Metals Co., Ltd., Republic of Korea), transferred to the laboratory, and stored at −20 • C until further analysis.
After rinsing with distilled water using a sieve (pore size, 200 µm), each zooplankton net sample was divided into two halves. The wet weight of half of the sample was measured, and six volumes of lysis buffer (Biosesang, Republic of Korea) were added for genomic DNA extraction. Genomic DNA was extracted with 700 µL of homogenized sample using an AccuPrep R genomic DNA extraction kit (Bioneer, Republic of Korea) following the manufacturer's instructions. The isolated DNA was quantified using a NanoDrop TM 1,000 Spectrophotometer (Thermo Fisher Scientific, United States) and qualified by agarose gel electrophoresis. The other half of the sample was used for the measurement of dry weight, according to a previous study (Jacobs and Grant, 1978). The comparisons of biodiversity between FishU and MiFish for ichthyoplankton were determined by NGS analysis. The library was constructed as mentioned above, but the first PCR products were pooled monthly and purified. Nextera libraries were sequenced from Illumina using 2 × 300 bp sequencing on MiSeq.

NGS Data Analysis
Using the CLC Genomic Workbench v.8.0 (CLC Bio, United States), the short reads (<100 bp) and low-quality sequences (QV <20) were trimmed from the raw data. Pairedend reads of FishU sequences were assembled using Mothur software (Schloss et al., 2009). From the assembled reads, the reads that overlapped by >6 bp, no mismatches, and sizes between 350 and 500 bp were filtered. The reads that could not be matched with options were discarded. The sequences of the forward and reverse primer regions were trimmed with one nucleotide mismatch option. Using the UCHIME software package (Edgar et al., 2011), operational taxonomic units (OTUs) were assigned at 99.7% identity from the obtained paired-end contigs, and chimeric sequences were removed based on de novo chimera detection. After OTUs with fewer than 10 contigs were removed, the species name was assigned using the BLASTn algorithm of BLAST + 2.2.38 (Camacho et al., 2009) against the NCBI non-redundant nucleotide database (accession date: 09/04/2019). The top-scored species name was assigned for each OTU with >97% sequence identity to the database. OTUs between 95 and 97% identity to the database were described as "genus name with highest score" followed by "sp." OTUs FIGURE 3 | Sample collection stations. The 12 sample collection stations were the coastal waters of three cities including Daejin (D1, D2, D3, and D4), Ayajin (A1, A2, A3, and A4), and Susan (S1, S2, S3, and S4). Maps created using Ocean data view v. 4.7.8 (Schlitzer, 2016, Ocean Data View, https://odv.awi.de).
with an identity of less than 95% were classified as "unknown." MiFish raw reads were processed using the MiFish pipeline with default settings (Sato et al., 2018). The OTUs obtained from MiFish were assigned to the highest blast score species at the species level with at least 97% sequence identity, and those between 95 and less than 97% sequences were assigned to the same genus. The representative haplotypes of two primers were designated in OTUs with abundances of more than 10% in each species, and these OTUs were subjected to phylogenetic analysis by maximum-likelihood algorithm using MEGA-X (Kumar et al., 2018). To compare the species resolution between the two primers, we calculated the genetic distance using MEGA-X based on the pairwise distances (p-distance).

Design of Fish-Specific Universal Primer, FishU
After comparing the mitochondrial rRNA sequences from 1808 fish species in the database, a universal primer set, FishU, was designed to amplify the variable region between the highly conserved mitochondrial 12S and 16S ribosomal sequences (Figure 1). The predicted size of PCR products amplified by the FishU primer ranged from 300 to 400 bp with an average of 373 bp, which was optimized to the MiSeq platform. Unusually large amplicons were identified in only 14 species, including Hoplolatilus cuniculus (818 bp), Centropyge multicolor (558 bp), and Centropyge joculator (558 bp).
The amplification efficiency of the newly developed FishU primers was compared with two previously designed fishspecific universal primers, ecoPrimer and MiFish. The recovered species numbers by the different mismatched bases within the primer sequences were compared (Table 1 and Supplementary Tables 2-4). When PCR was conducted without any mismatch, the highest number of species was amplified by the FishU primer set, in which 1,529 species out of a total of 1,683 examined bony fish species (Actinopteri) were able to be recovered (90.85%). Although the ecoPrimer set also amplified 82.06% of them, only 7.13% of the examined species were amplified using MiFish, and only 7.13% of the examined species were amplified using MiFish primers ( Table 1). When one mismatch was allowed, the percentage recovered by the MiFish primer increased up to 81.28% while 94 and 98.10% of the examined species could be amplified by ecoPrimer and FishU primer, respectively. As the mismatch numbers increased to three, the amplification success rates became similar among all three compared primers (94.18% by ecoPrimer, 97.39% by MiFish, and 98.46% by FishU). In contrast to the other two universal primers, the FishU primer also showed a high recovery percentage in the Chondrichthyes, up to 98.99% (Table 1). However, only 0-13.13% of the Chondrichthyes species could be amplified by ecoPrimer and MiFish primer sets with up to one mismatch. The taxonomic resolutions of the region amplified by the three universal primers were compared ( Table 2). The sequence variations and the length of the barcode region by each universal primer set are key to the high degree of taxonomic resolution. The percentage of correctly assigned species numbers by different universal primers was calculated. Barcodes using ecoPrimer showed the highest misidentification percentage by one nucleotide mismatch with 42.37% misassigned rates ( Table 2). By contrast, only 10.15% were misidentified by one nucleotide mismatch in the FishU primer. Even in five nucleotide mismatches, 78.56% of the barcodes by FishU were correctly assigned, and only 27.12% of those by ecoPrimer were adequately assigned. Compared with those by ecoPrimer, barcodes by MiFish showed a much higher accuracy in the taxon assignment ranging from 81.15 to 60.50% according to the number of mismatches. However, its rates were lower than those by FishU ( Table 2).
The taxon specificity of the three fish universal primers, namely ecoPrimer, MiFish, and FishU, were compared with four individual zooplankton, which were among the most abundant taxa in the zooplankton net samples (Figure 2). FishU amplified only Engraulis japonicus (415 bp), and ecoPrimer and MiFish showed a high degree of cross-reactivity in the examined invertebrate species, including Sagitta elegans, Euphausia pacifica, Metridia pacifica, and Penaeus monodon (Figure 2). The metabarcoding analysis of the zooplankton net samples was performed with three fish-specific universal primer sets to analyze the ichthyoplankton assemblage ( Table 3). Only ecoPrimer showed cross-reactivity with low proportions of Arthropoda (1.2 and 5.2%), Mammalian (0.002% in S2), and Enteropneusta (0.1 and 0.03%). In contrast, both MiFish and FishU presented only fish taxa. These results indicate that FishU primers are useful for examining the assemblage of ichthyoplankton directly from zooplankton net samples with a long barcode size, a low cross-reactivity, and a high specificity for fish taxa.

Comparative Analysis of Ichthyoplankton Survey in East/Japan Sea by Two Universal Primers: MiFish and FishU
Metabarcoding analyses of ichthyoplankton from the zooplankton net samples were conducted using two fishspecific universal primer sets (MiFish and FishU). A total of 1,816,585 and 4,383,919 raw reads from 6 month samples were generated using the MiFish and FishU primer set, respectively ( Table 4). On average, 203,551 merged reads were obtained using the MiFish primer after eliminating the low-quality reads, accounting for approximately 67% of the raw reads. All processed merged reads were identified as fish taxa (Table 4). Compared with those by MiFish, only 35% (220,909) of raw reads by the FishU primer were successfully merged, representing a 2.01-fold decrease compared with those by MiFish. An average of 75% of merged reads were identified as fish barcodes in those by FishU primer (Table 4). A total of 180 representative haplotypes were obtained from the MiFish primer set, including 42 species in 22 families in 6 orders. Forty-seven haplotypes from 33 species in 21 families in 9 orders were obtained by the FishU primer set, which was fewer than those from the MiFish primer set (Figures 4-6). Among them, only 16 species were commonly identified by the two primers, which were 38.1 and 48.5% of those by MiFish and FishU, respectively. In the order Perciformes, 46 (8 families) and 18 haplotypes (7 families) were obtained using the MiFish and FishU primers, respectively (Figure 4). Three (Seriola quinqueradiata, Chirolophis japonicus, and Uranoscopus japonicus) and two (Suruga fundicola and Chirolophis saitone) species were detected only by the FishU and MiFish primers, respectively ( Table 5). MiFish generated a higher number of haplotypes in each species compared with those by FishU primers. In particular, haplotypes with a low identity to the reference database, such as those in the Stichaeidae and Trichodontidae, were more difficult to assign to their correct species (Figure 4B). In contrast, multiple haplotypes were identified in only a limited number of species in those by FishU primers. Some haplotypes were assigned different species names, presumably based on different reference data. For example, only a single nucleotide was different between Chirolophis saitone by MiFish and Chirolophis japonicus by FishU (Figure 4). The reference sequence of C. saitone for the FishU barcode has not yet been deposited in the database, and it should be supplemented.
Haplotypes in the Scorpaeniformes, especially those in the families Cottidae and Psychrolutidae, showed a low identity to the reference database in both primers, indicating that the reference data in the family should be supplemented in the East Sea ( Figure 5). In addition, haplotypes in the family Sebastidae by MiFish showed a 100% identity to multiple species, indicating the low resolution of the MiFish region to distinguish those taxa (Figure 5B), which was also identified in a previous study (Yamamoto et al., 2017). In contrast, two sebastids, Sebastes owstoni and Sebastiscus marmoratus, were clearly assigned by FishU with 99% identity to the database (Figure 5A), which was supported by a previous conventional survey (Sohn et al., 2014). In addition, the assignment of species in Pleuronectidae were highly different between MiFish and FishU. Two species, P. yokohamae and G. zachirus, were assigned by the MiFish pipeline, and P. herzensteini and G. stelleri were assigned by FishU (Figure 6). The habitats of these two pleuronectid species by MiFish did not match the information of FishBase 2 . For example, G. zachirus is distributed from the Kuril Island to the Bering Sea coasts of Russia, not the coastal water of the East/Japan Sea, suggesting an incorrect species assignment by the MiFish pipeline. In addition to those in the family Sebastidae, species in the Pleuronectidae assigned by the MiFish pipeline should be manually checked before their use in regional/local surveys. 2 www.fishbase.org  Based on the metabarcoding analysis of zooplankton net samples using the two primers, the spawning period of each species was estimated in the East/Japan Sea ( Table 5). The average species richness by MiFish was 8.3, ranging from 5 to 13, and 9.8 by FishU, ranging from 4 to 14. Four species, including Konosirus punctatus, Engraulis japonicus, Ammodytes personatus, and Scomber japonicus, were detected with a similar pattern by both primers. Although the spawning periods for two related species, S. lalandi and S. quinqueradiata, were similar to each other from April to May along the Korean waters (Shiraishi et al., 2010;Miura et al., 2020), ichthyoplankton for S. quinqueradiata were detected only by the FishU primer in June ( Table 5). Chirolophis japonicus and Liza haematocheila were also detected in January and June as reported in previous studies (Kim et al., 2015;Park et al., 2015). In the species belonging to Pleuronectiformes, two species were identified in the genus Glyptocephalus: Glyptocephalus stelleri by FishU and Glyptocephalus zachirus by MiFish (Figure 6). According to the previous study, the main spawning season for G. stelleri is from April to May (Cha et al., 2008). Two species belonging to the genus Pseudopleuronectes were also detected from March to June: Pseudopleuronectes herzensteini by FishU and Pseudopleuronectes yokohamae by MiFish. The previously reported reproduction periods of P. herzensteini eggs also support the correct species assignment by the metabarcoding of the FishU primer (Lee and Kim, 2016). Collectively, the FishU primers outperformed MiFish in icthyoplankton metabarcoding analysis because of the longer barcode, allowing for a clearer species assignment and a broader taxon coverage, representing more taxa, including Synodontidae, Ophidiidae, Uranoscopidae, Cynoglossidae, and Triglidae.
The ichthyoplankton metabarcoding results and previous records for the spawning season were compared ( Table 5).
Because most of the studies on the spawning seasons have been conducted in commercially important species, including Scomber japonicus, Ammodytes personatus, and Engraulis japonicus, our comparisons were limited to these species (Table 5; shown in blue). The detection of ichthyoplankton generally reflects the spawning season. For instance, Sebastiscus marmoratus and Maurolicus muelleri showed spawning periods similar to those reported in previous studies (Yuuki, 1982;Takano et al., 1991). Metabarcoding analysis successfully discriminated the spawning season of two relative species, Ammodytes peronatus and Ammodytes hexapterus (Table 5). However, the onset time of larval detection in other species was often identified some months after the previous records. For instance, the spawning period of Konosirus punctatus was known to be from March (Kim et al., 2007a), and its larvae were first identified in June (Table 5). Similarly, two species belonging to the genera Scomber, S. japonicas, and S. australasicus, were detected in June, 3 months after their recorded reproduction periods. The spawning period of G. herzensteini was known from December to February (Park et al., 2007), but was detected from February to April.
Our ichthyoplankton metabarcoding analyses also revealed unreported spawning seasons of some species. For instance, we were able to detect the larvae of Cynoglossus robustus in the East/Japan Sea from January to April, and it was previously found on the southeastern coast of Korean peninsula from June to July (Park et al., 2013). In addition, the spawning seasons of several species, including Suruga fundicola, Anisarchus medius, Uranoscopus japonicus, Lepidotrigla microptera, and Neobythites sivicola, were first estimated in this study. Collectively, the metabarcoding analysis of ichthyoplankton shows great potential for understanding the reproduction of fish species around the Korean peninsula. However, further study is needed to produce data comparable to the conventional ichthyoplankton surveys.

DISCUSSION
In the present study, we designed a novel fish-specific universal primer set, FishU, for the metabarcoding analysis of ichthyoplankton and compared its reliability with two currently used universal primers, ecoPrimer and MiFish. Because most of the current metabarcoding analyses are conducted based on PCR analysis amplified by universal primers, the choice of the universal primer is crucial. FishU outperformed the other two fish-specific universal primers for the ichthyoplankton survey in several aspects. First, the FishU primer showed a high degree of taxon specificity, amplifying almost exclusively fish taxa directly from the zooplankton net samples. In addition to the sequence specificity within the primer region, the barcode region amplified by FishU contained tRNA val flanked by the 12S and 16S regions, a characteristic structure of the fish mitogenome (Figure 2). Therefore, each primer can land on different genes, 12S and 16S, respectively, lowering the chance of amplifying different taxa. By contrast, ecoPrimer and MiFish primers have been designed within the 12S gene, increasing the chance of amplifying a similar sequence of other taxa (Figures 1, 2). Copepods, euphausiid, and Chaetognatha species are among the major zooplankton taxa (Iguchi, 2004), and even a low degree of cross-reactivity to those species would be problematic for the metabarcoding of ichthyoplankton directly from the net sample.
In addition to taxon specificity, the longer barcodes of the FishU primer contributed to a higher accuracy in species identification by both in silico and zooplankton net sample analyses. The average length of the barcodes generated by FishU was 373 bp, which represents a 3.5-and 2.2-fold increase compared with those generated by ecoPrimer (106 bp) and MiFish (172 bp), respectively. Although 200 bp or smaller sizes have been reported as the ideal length for metabarcoding to maximize the recovered species numbers (Coissac et al., 2012;Clarke et al., 2017), these small-sized barcodes may trade off the accuracy of species identification. For example, a single base substitution in the barcodes by ecoPrimers resulted in almost half of the barcode not being able to correctly assign species (Table 2). In addition, the short length of the MiFish barcode may not be variable enough to discriminate between closely related species, especially in Pleuronectiformes or Scorpaeniformes, as reported in previous studies (Yamamoto et al., 2017). By contrast, Pseudopleuronectes herzensteini and Sebastes owstoni in those orders were accurately assigned only by FishU, which was supported by previous studies (Sohn et al., 2014;Lee and Kim, 2016). In addition to these two orders, two related species in each genus, Alcichthys alcicornis and Alcichthys elongatus and Maurolicus japonicus and Maurolicus muelleri, should be reexamined.  • relative proportion of each month was 50-100%, relative proportion of each month was 1-50%, • relative proportion of each month was 0-1% (except 0%). The spawning periods from the previous studies are blue-shaded.
The application of the MiFish pipeline requires careful attention from those who would like to use the data directly from its bioinformatic analysis for the analysis of local fish assemblages. Although much higher haplotypes (180) and species numbers (42) were obtained by the MiFish pipeline in this study compared with those by the FishU platform (51 haplotypes and 31 species), most of the species assigned only by MiFish showed a low similarity to the database ( Table 5). Those assigned as genus names, such as Konosirus sp., Engraulis sp., or Ammodytes sp., may have come from the short barcode length. When merging each pair read generated by the MiSeq platform, it was especially challenging to eliminate chimeric sequences for multiple haplotypes, species with a high degree of genetic similarity, or lack of a reference database. For instance, three species in the genus Ammodytes could not be adequately assigned by the MiFish pipeline, generating a high number of haplotypes with low similarities, mainly due to the genetic similarity among those species and the lack of the local reference sequences in the database (Figure 4). The chimeric sequences generated by MiFish would be problematic, exaggerating the species numbers for countries without a sufficiently large reference database. In contrast, the barcodes by FishU showed no chimeric haplotypes with long read lengths. In fact, all the species exclusively identified by FishU primers showed a high sequence identity, with an identity of more than 99% to the database, which included Clupea pallasii, Neobythites sivicola, Seriola quinqueradiata, Uranoscopus japonicus, and Radulinopsis derjavini ( Table 5). This result indicates that FishU provides species information with a high degree of accuracy in ichthyoplankton metabarcoding analysis. The longer barcode size by FishU was not a problem, at least for ichthyoplankton metabarcoding. Instead, FishU outperformed the other two fish-specific universal primers in the ichthyoplankton survey with higher accuracy and reliability.
Compared with the other two fish-specific universal primers, the use of FishU for the ichthyoplankton metabarcoding has several advantages, including higher taxon coverage and specificity and accurate species identification. However, there are still several shortcomings in adopting the primer. The barcode region by FishU includes the tRNA-valine flanked by the partial 12S rRNA and 16S rRNA, whose reference sequences are fewer than those of the typically used barcodes, such as COI or MiFish region. Notably, there is much less information regarding those that have little commercial value. Because taxonomic identification in metabarcoding is entirely dependent on the reference database, the reference sequences for the FishU region, especially those for the local or indigenous species, should be supplemented. The recent fast growth of complete mitochondrial DNA information of various fish taxa is also helpful to compensate for these shortcomings (Iwasaki et al., 2013). Another factor to consider in choosing FishU is the low yield of merged reads from the raw data ( Table 5). The proportions of the average merged fish barcode numbers generated by the FishU primer were 35%, which was twofold lower than that by MiFish. These results are mainly due to the low merged rates for the longer barcodes by the FishU primer . These low merged-read rates can be compensated by the increased raw reads. To obtain processed merged reads similar to those of MiFish in this study, approximately twofold higher raw read numbers were required in the FishU platform (Table 4). Because longer barcodes have many advantages in ichthyoplankton metabarcoding, metabarcoding analysis by FishU would trade off the shortcoming in the cost. Storage after sample collection may also affect the recovery of longer barcodes by FishU, and we cannot rule out that the low merged-read numbers from January to April may have arisen due to the degradation of DNA after a long period of storage (Table 4) because it is widely known that DNA can be degraded over long storage times (Vamos et al., 2017). Therefore, higher merged reads could be achieved by metabarcoding immediately after the extraction of total genomic DNA.
Here, we attempted to determine the spawning season based on the metabarcoding analysis of ichthyoplankton from zooplankton net samples. The metabarcoding analysis of ichthyoplankton provides many useful data to understand the reproduction of fish, including details of their spawning sites, reproductive periods, and population structures, at relatively low cost and labor. The large amount of data obtained by metabarcoding analysis can be used for various purposes not previously possible using conventional ichthyoplankton surveys, including the effects of climate change during the spawning season in each species, an understanding of the reproduction of rare species, or the genetic populations of ichthyoplankton. In addition, the transport routes of fish eggs and larvae in the East/Japan Sea can also be understood by its larval metabarcoding analysis. Waters in the East/Japan Sea are the place where two main currents along the Korean peninsula, the NKCC from the north and EKWC, collide. Therefore, transportation and settlement of ichthyoplankton are largely dependent on the dynamics of these currents. In fact, we were able to identify a small number of tropical and cold-water fish larvae, which are further used as indicators of the spatiotemporal dynamics of these currents. For instance, Saurida microlepis and Chelidonichthys kumu have been conveyed from the South Pacific Ocean by the Kuroshio warm current and Tsushima current entering the East/Japan Sea through the Korean Strait. Because the Kuroshio current begins with the water around the Philippines and Taiwan, which is known to be the spawning site of many fish species, the abundance of these tropical fish would explain the effects of the warm current on fisheries in the East/Japan Sea. In contrast, some species detected during the winter season may have been conveyed from the NKCC. For example, Anisarchus medius is a demersal species, preferring water at temperatures below 0 • C, and its larvae were only found in February although A. hexapterus was originally distributed in cold waters from Arctic Alaska to the northern Pacific. Therefore, changes in the abundance of tropical or cold ichthyoplankton help explain the effects of the dynamics of the two main currents on fisheries in the East/Japan Sea.
Although it is a useful tool for the ichthyoplankton survey, metabarcoding analysis is not likely to replace currently used methods in the near future. In addition to diversity, several data, such as morphological or developmental parameters, can be obtained only using the conventional ichthyoplankton survey. For example, egg or larval fish abundances cannot be obtained by metabarcoding analysis using the current methodology. It is also unclear whether metabarcoding analysis can explain the behavioral characteristics, such as diurnal vertical migrations, of fish larvae or spawning ecology. Quantitative analysis is another issue to establish the correlation between metabarcoding analysis and conventional survey. Although it is generally accepted that there is a quantitative relationship between the sequence reads and biomass (Deagle et al., 2019;Lamb et al., 2019), a clear consensus has yet to be reached with regards to the quantification methodology and continues to vary between research groups. Furthermore, there are many chances for bias throughout the metabarcoding pipeline from the sample collection to the bioinformatic processes, making it more difficult to achieve an accurate quantification (Juen and Traugott, 2006;Plummer et al., 2015;Smith et al., 2017). We also identified the delayed detection of ichthyoplankton in many species compared with the previously known spawning season ( Table 5). The reproductive period of a species is usually estimated by the combination of gonad maturation and ichthyoplankton survey (Murua et al., 2003;Lowerre-Barbieri et al., 2011). The delayed detection can be explained by the lower chance of detection in eggs, which have lower copies of DNA, compared with those of larvae. As cell numbers increase from eggs to fish larvae during the developmental stages, the chance of detection is much higher for the larval stage than for eggs (Takeuchi et al., 2019). Therefore, further studies should be conducted to identify the spawning seasons more precisely via the metabarcoding analysis of ichthyoplankton.
Despite its current limitations, the metabarcoding analysis of ichthyoplankton remains a promising strategy to supplement conventional surveys. As the data accumulates by comparative analysis with conventional microscopic observation, this technique provides more reliable data. Once a universal and automated metabarcoding platform is established, cost-effective and long-term ichthyoplankton surveys with a higher degree of statistical reliability will be possible, which would provide useful information for the scientific management of local fish resources.

DATA AVAILABILITY STATEMENT
The sequencing data from this study has been deposited in BioProject (accession: PRJNA682864).

AUTHOR CONTRIBUTIONS
AK designed the experiments, performed the experiments, analyzed the data, prepared figures and tables, and wrote the manuscript. T-HY designed the experiments, analyzed the data, and prepared figures and tables. CL and C-KK conceived and designed the experiments, contributed reagents, materials, and analysis tools. H-WK conceived and designed the experiments, performed the experiments, analyzed the data, contributed reagents, materials, analysis tools, prepared figures and tables, authored and reviewed drafts of the manuscript, and approved the final draft. All authors contributed to the article and approved the submitted version.