Biodiversity and Community Structure of Mesozooplankton in the Marine and Coastal National Park Areas of Korea

: Zooplankton communities are useful bioindicators that can provide information on the changes occurring in marine ecosystems. Therefore, investigation of zooplankton communities in marine and coastal national parks is essential. However, the surveys of zooplankton communities using morphological identiﬁcation require considerable time and labor. Metabarcoding is a practical alternative that can detect various taxa simultaneously. In this study, metabarcoding was newly applied along with the traditional morphological identiﬁcation to establish a method for zooplankton community survey in the Marine and Coastal National Park areas of Korea. By comparing the results of these two identiﬁcation methods, the strengths and limitations of metabarcoding were veriﬁed with the zooplankton communities appearing in these areas. The sensitive detection capability of metabarcoding enabled the identiﬁcation of potential bioindicator taxa associated with external factors (e.g., water temperature, salinity, topography, and chlorophyll a concentration) in these national parks. We propose the use of metabarcoding for e ﬃ cient surveys of mesozooplankton communities in the Marine and Coastal National Parks to establish monitoring of bioindicator taxa. It is also necessary to continuously search for taxa with high research value in these national parks using metabarcoding. Establishing an ongoing monitoring system that employs this approach can provide an e ﬀ ective tool for managing marine ecosystems in the Marine and Coastal National Parks. Author Contributions: Conceptualization, H.K., C.-R.L. and W.K.; methodology, H.K., C.-R.L., S.-k.L. and S.-Y.O.; validation, H.K. and W.K.; formal analysis, H.K.; investigation, H.K, C.-R.L. and S.-k.L.; resources, C.-R.L. and W.K.; data curation, H.K.; writing—original draft preparation, H.K.; writing—review and editing, C.-R.L., S.-k.L., S.-Y.O. and W.K.; visualization, H.K.; supervision, W.K.; project administration, C.-R.L.; funding acquisition,


Introduction
Marine ecosystems are changing as a result of global climate change and industrialization in coastal areas. Unlike terrestrial ecosystems, marine ecosystems can be difficult to access and are influenced by a unique set of external factors, including the degree of light transmission, oxygen concentration, water masses, currents, and salinity, which complicate assessments and predictions. As marine ecosystems change, bioindicators respond by changing their morphological or cellular structure, metabolic processes, behaviors, and communities [1][2][3]. Due to these characteristics, studying bioindicators that can confirm and monitor the changes in the marine ecosystem is becoming important worldwide [4,5].
Zooplankton represent the primary and secondary consumers in the aquatic food chain and are some of the most abundant and ubiquitous taxa in aquatic ecosystems [6][7][8]. The spatial and temporal distribution of zooplankton communities fluctuate in response to environmental changes in marine ecosystems, such as variations in temperature and salinity [9,10]. Therefore, zooplankton are useful bioindicators for detecting environmental changes in the marine ecosystem [11][12][13][14][15][16]. However, the investigation of zooplankton communities using traditional morphological identification requires high taxonomic knowledge as well as considerable time and labor [17]. Additionally, it can be difficult to identify the lowest taxonomic ranks (i.e., genus and species) as some zooplankton have ambiguous morphological characteristics, particularly in the larval stages [18]. The development of next-generation sequencing (NGS) technology led to the advent of metabarcoding, a method that can quickly and simultaneously detect any taxa within the database [19,20]. Thus, large-scale marine ecological surveys were made possible with bulk-sample metabarcoding [21][22][23][24]. As these advantages were revealed, many researchers conducted comparative studies to confirm that metabarcoding was effective for ecological surveys when compared to traditional morphological identification. Most previous studies report that metabarcoding detects more taxa than morphological identification methods. Additionally, differences in communities can be distinguished and identified more efficiently. However, it is still difficult to achieve accurate biodiversity and species composition surveys with metabarcoding because of the potential for distortion of species abundance as a result of technical biases and false negatives [25][26][27][28][29].
The national parks of South Korea are designated as regions that represent the natural ecosystems or the natural and cultural landscapes of Korea. According to the Korea National Park Service website (http://www.knps.or.kr), a total of 22 national parks are designated in South Korea [30]. Among these, only four are marine and coastal national parks. Taeanhaean National Park and the Byeonsan-bando National Park are situated along the Yellow Sea coast. Dadohaehaesang National Park includes areas of both the Yellow Sea coast and Southern Sea coast of Korea, and Hallyeohaesang National Park is located on the Southern Sea coast of Korea. The Marine and Coastal National Parks aim to preserve the valuable and highly diverse ecosystems within them. As such, ecological study and efficient management of the Marine and Coastal National Parks are essential.
In this study, metabarcoding was newly applied along with traditional morphological identification to establish a method for zooplankton community surveys in the Marine and Coastal National Park areas of Korea. Mesozooplankton (>200 µm) were selected as the target organisms because there were many previous studies conducted using regular zooplankton surveys at the Marine and Coastal National Park areas, thus allowing for comparison of the identification results of metabarcoding with those of morphological identification. The mesozooplankton communities in the Marine and Coastal National Park areas were compared and analyzed according to sea area and location because the two areas (Yellow Sea and Southern Sea of Korea) and four locations (Taean, Byeonsan, Dadohae, and Hallyeo areas) included representative diverse marine environments with variations in depth, topography, effects of currents, and inflow of freshwater [31][32][33]. The main objective of this study was to perform a metabarcoding analysis of the biodiversity and community structure of mesozooplankton communities in the Marine and Coastal National Park areas. First, we verified the strengths and limitations of metabarcoding by comparing the results with those obtained by morphological identification. Second, bioindicator taxa associated with spatial and environmental characteristics were identified based on the metabarcoding analysis. Finally, we discussed the potential of metabarcoding analysis as an efficient method to monitor the zooplankton community.

Mesozooplankton Sample and Spatial and Environmental Data Collection
Mesozooplankton samples were obtained from a spring season survey during "A Survey on Marine Ecosystems of the Marine and Coastal National Park Areas of Korea" conducted by the Marine Research Center of the Korea National Park Service from May to June 2019.
Sampling was conducted at 58 sampling stations, including the sampling stations in the four Marine and Coastal National Parks and adjacent sea areas ( Figure 1). All sampling stations were Sampling was conducted at 58 sampling stations, including the sampling stations in the four Marine and Coastal National Parks and adjacent sea areas ( Figure 1). All sampling stations were designated categories according to the sea area and location. The study area was divided into two sea areas and four locations (Taean, Byeonsan, Dadohae, and Hallyeo areas), based on the standard line drawn at 225° from Haenamgak (34°17'33.09" N, 126°31'26.02" E) of the Korea Hydrographic and Oceanographic Agency and areas under the jurisdictions in the Marine and Coastal National Parks, respectively. A 200 μm mesh conical net with a 60 cm diameter mouth was lowered vertically to the bottom of each sampling station and then raised at a rate of 1 m/s. A flowmeter (Hydrobios, 438115) was attached to the entrance of the net to measure the amount of seawater filtered. Sampling was performed in duplicate at each sampling point; one of the obtained samples was stored in 4% formalin solution for morphological identification and the other in 99% ethanol for DNA extraction and molecular identification. Spatial and environmental data were also obtained to verify the relationships with the zooplankton community structure and distribution. At each sampling station, spatial data were obtained using longitude and latitude data from a global positioning system (GPS). Environmental parameters, such as water temperature and salinity and depth, were measured at each sampling station using a SBE 9plus conductivity-temperature-depth (CTD) instrument (Sea-Bird Electronics Inc., Bellevue, WA, USA). Chlorophyll samples were collected at each sampling station by filtering both the surface and benthic seawater through glass fiber filters (GF/F; Ø 25 mm, pore size 0.7 μm, Whatman, Maidstone, England) for chlorophyll a analysis. The filter papers were then placed in lightresistant containers with 90% acetone and frozen until the chlorophyll a was extracted. The extracted chlorophyll samples were transferred to test tubes, and chlorophyll a concentration measured using Spatial and environmental data were also obtained to verify the relationships with the zooplankton community structure and distribution. At each sampling station, spatial data were obtained using longitude and latitude data from a global positioning system (GPS). Environmental parameters, such as water temperature and salinity and depth, were measured at each sampling station using a SBE 9plus conductivity-temperature-depth (CTD) instrument (Sea-Bird Electronics Inc., Bellevue, WA, USA). Chlorophyll samples were collected at each sampling station by filtering both the surface and benthic seawater through glass fiber filters (GF/F; Ø 25 mm, pore size 0.7 µm, Whatman, Maidstone, England) for chlorophyll a analysis. The filter papers were then placed in light-resistant containers with 90% acetone and frozen until the chlorophyll a was extracted. The extracted chlorophyll samples were transferred to test tubes, and chlorophyll a concentration measured using a fluorophotometer (10AU, Turner Designs, Sunnyvale, CA, USA). The environmental variables at each sampling station were measured from the surface to the benthic layer and then averaged.

Morphological Identification and Metabarcoding Process
Formalin-fixed mesozooplankton samples were transported to the laboratory for species identification and counting. In the laboratory, mesozooplankton samples were divided into subsamples of 500-1000 individuals using a Folsom zooplankton splitter. Each subsample was counted and identified in a Bogorov counting chamber under a Leica M165C stereomicroscope (Leica Microsystems, Wetzlar, Germany). Taxonomy experts identified most of the copepods to the species level, but some individuals that were difficult to identify at the species level were classified to the lowest possible taxonomic level (Tables S4-S6).
The samples for DNA extraction were vortexed at maximum speed and then centrifuged for 5 min at 13,000× g. Subsequently, the supernatant was removed and incubated at room temperature until ethanol had completely evaporated. DNA was extracted from the pellets using a Qiagen DNeasy Blood & Tissue Kit (Qiagen, Valencia, CA, USA) following the manufacturer's instructions. In the last step, each eluted DNA sample was recombined according to the sampling station. Three PCR replicates were performed for DNA amplification and each DNA sample was diluted by 1:10. We chose a primer set to target the V9 region of the 18S ribosomal DNA, because it has the ability to detect the whole of zooplankton communities. Also, it is one of the most commonly used to investigate zooplankton using metabarcoding [21,[34][35][36][37][38][39][40]. The 18S ribosomal DNA (rDNA) V9 variable region was amplified using the 1391F (5 -GTACACACCGCCCGTC-3 ) and EukBr (5 -TGATCCTTCTGCAGGTTCACCTAC-3 ) primers [35]. The PCR amplification was performed as follows: 3 min at 94 • C, 35 cycles of 45 s at 94 • C, 45 s at 65 • C, 30 s at 57 • C, and a final extension step of 10 min at 72 • C. The amplified PCR products were confirmed via electrophoresis and pooled together for each sample. The amplified PCR products were then purified using the QIAquick PCR Purification Kit (Qiagen, Valencia, CA, USA) and paired-end sequencing was performed at Macrogen Inc. (Seoul, Korea) on the Illumina MiSeq platform.
The 18S rDNA sequencing data produced by Illumina MiSeq were analyzed using Quantitative Insights Into Microbial Ecology (QIIME) v 1.9.1 [41]. Forward and reverse sequences were concatenated into one read using PEAR with the default parameters [42]. Short (<200 bp) or low-quality assembled reads (Q < 30) were discarded and only the assembled reads were included in the bioinformatics analysis. Detection of chimeric sequences and operational taxonomic unit (OTU) clustering were performed using VSEARCH [43]. Chimeric sequences and singleton sequences were excluded from the analysis. All OTUs were clustered with 97% similarity and the most abundant sequence was selected in each OTU. These representative sequences were assigned taxonomic information by comparing the 18S rDNA eukaryotic database from the NCBI GenBank parsed using Biopython (http://www.biopython.org). In cases where the assigned taxonomic information of OTUs was unclear (e.g., uncultured/environmental sample sequences), it was inferred with the taxonomic information of the closest assigned species, considering lowest similarity thresholds for copepod taxonomic resolution (more than 96% for identification to family level; 85% or more to phlyum level) [44]. To revise the number of reads distorted by the technical bias problem, rarefaction for biodiversity analysis was conducted considering a minimum number of reads.

Bioinformatics and Statistical Analysis
All bioinformatics and statistical analyses were conducted and visualized with plots using ggplot2, Phyloseq, ggplot2, vegan, pairwise Adonis, dunn.test, rcompanion, and ade4 packages in R v 3.5.1 [45][46][47][48][49][50][51][52]. All p-value adjustments were applied as the false discovery rate (FDR) [53]. Taxonomic information and species counts (read counts) obtained using the morphological identification and metabarcoding were converted to Biological Observation Matrix (BIOM) format for the analysis of biodiversity and community structure, respectively. The indices of richness (observed species (OTUs) and Chao1), diversity (Shannon's diversity), and evenness (equitability) for each BIOM file were calculated using QIIME script (alpha_diversity.py). Statistical significances in the biodiversity indices for the sea area variables were determined by the Wilcoxon rank sum test. The Kruskal-Wallis test and pairwise comparisons were conducted to identify significant differences among the location variables, with the Dunn's test as a post hoc test.
To examine the differences between mesozooplankton community structures, the unweighted pair group method with arithmetic averages (UPGMA) was analyzed based on Bray-Curtis dissimilarities. To test the similarity in the zooplankton community structure identified by the two methods, two UPGMA cluster trees were compared with formed zooplankton communities. Procrustes analysis was conducted with 1000 permutations using the protest function. Constrained analysis of principal coordinates (CAP) was also performed to identify the relationships between zooplankton community structures and the following categories: sea area (Yellow Sea and Southern Sea of Korea), location (Taean, Byeonsan, Dadohae, and Hallyeo area), and spatial, environmental variables (latitude, longitude, water temperature, salinity, and chlorophyll a concentration). Statistical differences from the CAP analysis and among community structures were evaluated using ANOVA and the pairwise Adonis with the test of 999 random permutations, respectively. The taxonomic compositions of the mesozooplankton communities identified with the two methods were compared based on the phylum level and the most frequently detected family level.

Environmental Characteristics in the Marine and Coastal National Park Areas in Korea
During the survey period, the environmental data collected from of the Marine and Coastal National Park areas were compared according to the sea area and location (Table 1). Among the sampling stations, the average water temperature was higher in the Southern Sea of Korea than the Yellow Sea. The salinity of the Taean and Byeonsan areas was lower than that of Dadohae and Hallyeo areas. The average chlorophyll a concentration in the Yellow Sea was higher compared with that in the Southern Sea of Korea. The average depth was the greatest in the Hallyeo area and the lowest in the Byeonsan area. The deepest individual sampling point was N2 (76.05 m) in the Dadohae area and the shallowest was H2 (3.29 m) in the Dadohae area.

Mesozooplankton Biodiversity Analysis
We performed a comparison between the number of species identified by the morphological identification and the number of OTUs based on the similarities of sequences in metabarcoding ( Table 2). This is an indirect comparison because the 97% similarity distance measures used for the OTU clustering have insufficient resolution to distinguish between zooplankton species. With morphological identification, a total of 79 taxa were identified in mesozooplankton samples from the Marine and Coastal National Park areas. Fifty-five taxa were found in the Yellow Sea, 73 taxa were found in the Southern Sea of Korea, and 52 taxa were shared by both sea areas. The number of taxa identified in each location was as follows: 37 in the Taean area, 30 in the Byeonsan area, 57 in the Dadohae area, and 61 in the Hallyeo area. Using Illumina MiSeq sequencing, 18S rDNA sequencing data were produced from  (Table S4). The biodiversity indices were compared by sea area and location ( Figure 2, Tables S1 and S2). The results of morphological identification and metabarcoding showed similar patterns in biodiversity indices according to the sea area. Although the diversity indices calculated from the two methods were slightly different, the richness and evenness of the zooplankton communities were the same ( Figure 2, and Table S1). In contrast, the patterns of all biodiversity indices calculated among locations were completely different when using the morphological identification and metabarcoding ( Figure 2, and Table S2). The zooplankton richness of the Hallyeo area using the morphological identification was high compared to other areas, but when the metabarcoding approach was used, there were no statistical differences among locations. Comparing the diversity indices and evenness of zooplankton between the two methods, these biodiversity indices were distinctly lower in the Byeonsan area than the Dadohae area when calculated using the morphological identification results. However, these biodiversity indices calculated using metabarcoding were significantly higher in the Hallyeo area than the Dadohae area. Statistical differences in the biodiversity indices according to the sea area were calculated using the Wilcoxon rank sum test. As a post hoc analysis, all p-values were corrected using the Benjamini-Hochberg procedure (** p < 0.01, * p < 0.05, N.S. no significance). In the case of location, the significance of biodiversity indices was calculated using the Kruskal-Wallis test. As a post hoc analysis, pairwise comparisons were conducted using Dunn's test. The results for Dunn's test were marked using the same letter for values that were not significantly different from each other. OTUs are operational taxonomic units. Boxplots for the mesozooplankton biodiversity indices were calculated using morphological identification and metabarcoding results according to the sea area and location. Statistical differences in the biodiversity indices according to the sea area were calculated using the Wilcoxon rank sum test. As a post hoc analysis, all p-values were corrected using the Benjamini-Hochberg procedure (** p < 0.01, * p < 0.05, N.S. no significance). In the case of location, the significance of biodiversity indices was calculated using the Kruskal-Wallis test. As a post hoc analysis, pairwise comparisons were conducted using Dunn's test. The results for Dunn's test were marked using the same letter for values that were not significantly different from each other. OTUs are operational taxonomic units.

Mesozooplankton Community Analysis
In both methods, mesozooplankton communities in the Marine and Coastal National Park areas were grouped into three clusters (Figure 3). Although there were some differences in the mesozooplankton samples that belonged to the cluster, there were similar patterns: Cluster 1 mainly contained the zooplankton samples from the Dadohae area; Cluster 2 tended to consist of zooplankton samples from the Hallyeo area, in addition to samples from the eastern parts of the Dadohae area; and the zooplankton samples of the Yellow Sea (included in the Taean and Byeonsan areas) formed Cluster 3. Through the Procrustes analysis, we confirmed that there was a significant correlation between the mesozooplankton communities formed by the morphological identification and metabarcoding (m12 squared = 0.80; correlation value = 0.44; p-value = 0.001).

Mesozooplankton Community Analysis
In both methods, mesozooplankton communities in the Marine and Coastal National Park areas were grouped into three clusters (Figure 3). Although there were some differences in the mesozooplankton samples that belonged to the cluster, there were similar patterns: Cluster 1 mainly contained the zooplankton samples from the Dadohae area; Cluster 2 tended to consist of zooplankton samples from the Hallyeo area, in addition to samples from the eastern parts of the Dadohae area; and the zooplankton samples of the Yellow Sea (included in the Taean and Byeonsan areas) formed Cluster 3. Through the Procrustes analysis, we confirmed that there was a significant correlation between the mesozooplankton communities formed by the morphological identification and metabarcoding (m12 squared = 0.80; correlation value = 0.44; p-value = 0.001). Using CAP analysis, the mesozooplankton communities in the Marine and Coastal National Park areas detected using the two methods were significantly different according to the sea area (morphological identification: p-value = 0.001, explanatory power = 16.6%; metabarcoding: p-value = 0.001, explanatory power = 29.0%) and location (morphological identification: p-value = 0.001, explanatory power = 25.0%; metabarcoding: p-value = 0.001, explanatory power = 40.1%) (Figure 4ad). According to the pairwise Adonis test, all mesozooplankton communities formed by the CAP analysis were significantly different (Table S3). In the contrast, taxonomic compositions between mesozooplankton communities differed depending on the identification method. At the phylum level, the identification results of the morphological identification and metabarcoding confirmed that Arthropoda was the largest taxon in the Marine and Coastal National Park areas (Figure 5a,b). However, while more Myzozoa were identified using morphological identification than metabarcoding, Cnidaria were conspicuously detected using metabarcoding. Interestingly, Rotifera Using CAP analysis, the mesozooplankton communities in the Marine and Coastal National Park areas detected using the two methods were significantly different according to the sea area (morphological identification: p-value = 0.001, explanatory power = 16.6%; metabarcoding: p-value = 0.001, explanatory power = 29.0%) and location (morphological identification: p-value = 0.001, explanatory power = 25.0%; metabarcoding: p-value = 0.001, explanatory power = 40.1%) (Figure 4a-d). According to the pairwise Adonis test, all mesozooplankton communities formed by the CAP analysis were significantly different (Table S3). In the contrast, taxonomic compositions between mesozooplankton communities differed depending on the identification method. At the phylum level, the identification results of the morphological identification and metabarcoding confirmed that Arthropoda was the largest taxon in the Marine and Coastal National Park areas (Figure 5a,b). However, while more Myzozoa were identified using morphological identification than metabarcoding, Cnidaria were conspicuously detected using metabarcoding. Interestingly, Rotifera were detected only by the metabarcoding method. Myzozoa and Cnidaria were found more prominently in the Hallyeo area compared with other locations, while Rotifers were detected more in the Taean and Byeonsan areas. Differences in the taxonomic composition of taxa identified using the two methods were more apparent when compared at the major family level (Figure 5c,d). The proportions of Acartiidae, Corycaeidae, Noctilucaceae, Oikopleuridae, and Podonidae identified applying the morphological identification were higher than when applying metabarcoding. In contrast, more Calanidae, Centropagidae, Diphyidae, Euphausiidae, Mysidae, Paracalanidae, and Sagittidae were detected with metabarcoding. Based on the results of the two identification methods, the taxonomic compositions of mesozooplankton communities in the Marine and Coastal National Park areas were compared according to the sea area and location. In the Taean area, both Centropagidae and Podonidae were more dominant compared to the other areas, and in the Byeonsan area, Acartiidae was more abundant compared to other areas. Paracalanidae was often observed in samples from the Southern Sea of Korea (Dadohae and Hallyeo areas). Oithonidae was also more common in two areas of the Southern Sea of Korea compared to the other areas. Calanidae, Euphausiidae, and Mysidae were identified more in the Dadohae area than in other areas. In the Hallyeo area, Notilucaceae accounted for nearly half of the mesozooplankton community when identified using the morphological identification, while Diphyidae and Sagittidae were also detected using metabarcoding. were detected only by the metabarcoding method. Myzozoa and Cnidaria were found more prominently in the Hallyeo area compared with other locations, while Rotifers were detected more in the Taean and Byeonsan areas. Differences in the taxonomic composition of taxa identified using the two methods were more apparent when compared at the major family level (Figure 5c,d). The proportions of Acartiidae, Corycaeidae, Noctilucaceae, Oikopleuridae, and Podonidae identified applying the morphological identification were higher than when applying metabarcoding. In contrast, more Calanidae, Centropagidae, Diphyidae, Euphausiidae, Mysidae, Paracalanidae, and Sagittidae were detected with metabarcoding. Based on the results of the two identification methods, the taxonomic compositions of mesozooplankton communities in the Marine and Coastal National Park areas were compared according to the sea area and location. In the Taean area, both Centropagidae and Podonidae were more dominant compared to the other areas, and in the Byeonsan area, Acartiidae was more abundant compared to other areas. Paracalanidae was often observed in samples from the Southern Sea of Korea (Dadohae and Hallyeo areas). Oithonidae was also more common in two areas of the Southern Sea of Korea compared to the other areas. Calanidae, Euphausiidae, and Mysidae were identified more in the Dadohae area than in other areas. In the Hallyeo area, Notilucaceae accounted for nearly half of the mesozooplankton community when identified using the morphological identification, while Diphyidae and Sagittidae were also detected using metabarcoding.

Potential Bioindicator Taxa Detection Using Metabarcoding
Morphological identification and metabarcoding were compared to identify potential bioindicator taxa reflecting spatial and environmental characteristics in the Marine and Coastal National Park areas. A CAP analysis revealed that the associations between mesozooplankton communities and all variables produced similar results using both methods. The mesozooplankton communities identified by both methods were significantly affected by all spatial and environmental variables (morphological identification: p-value = 0.001, explanatory power = 36.7%; metabarcoding: p-value = 0.001, explanatory power = 49.8%) (Figure 6a,b). Each mesozooplankton community cluster exhibited significant differences when using both methods (p-values = 0.001 for all clusters). Of the three community clusters formed, Cluster 1 exhibited no correlation between the external variables we obtained and the mesozooplankton community. In contrast, Clusters 2 and 3 were related to spatial and environmental variables. Cluster 2 was correlated with longitude, water temperature, and salinity; latitude and chlorophyll a concentration were correlated with Cluster 3. The taxonomic compositions between mesozooplankton community clusters formed by constraining spatial and environmental variables was shown in CAP analysis (Figure 6c,d). Paracalanidae, which was dominant in the Southern Sea of Korea, was more abundant in Cluster 1 and Cluster 2 compared to Cluster 3. A larger number of Calanidae were identified in Cluster 1 by both methods. Oikopleuridae and Oithonidae were frequently observed through morphological identification in Cluster 1, while Euphausiidae and Mysidae were detected more in Cluster 1 with metabarcoding. Notilucaceae, Diphyidae, and Sagittidae, which were associated with the Hallyeo area, were more common in Cluster 2 than other mesozooplankton clusters using both methods. The phylum Rotifera included in other with Acartiidae, Podonidae, and Centropagidae, which were more dominant in the Yellow Sea, were identified more in Cluster 3 by metabarcoding. Morphological identification and metabarcoding were compared to identify potential bioindicator taxa reflecting spatial and environmental characteristics in the Marine and Coastal National Park areas. A CAP analysis revealed that the associations between mesozooplankton communities and all variables produced similar results using both methods. The mesozooplankton communities identified by both methods were significantly affected by all spatial and environmental variables (morphological identification: p-value = 0.001, explanatory power = 36.7%; metabarcoding: p-value = 0.001, explanatory power = 49.8%) (Figure 6a,b). Each mesozooplankton community cluster exhibited significant differences when using both methods (p-values = 0.001 for all clusters). Of the three community clusters formed, Cluster 1 exhibited no correlation between the external variables we obtained and the mesozooplankton community. In contrast, Clusters 2 and 3 were related to spatial and environmental variables. Cluster 2 was correlated with longitude, water temperature, and salinity; latitude and chlorophyll a concentration were correlated with Cluster 3. The taxonomic compositions between mesozooplankton community clusters formed by constraining spatial and environmental variables was shown in CAP analysis (Figure 6c,d). Paracalanidae, which was dominant in the Southern Sea of Korea, was more abundant in Cluster 1 and Cluster 2 compared to Cluster 3. A larger number of Calanidae were identified in Cluster 1 by both methods. Oikopleuridae and Oithonidae were frequently observed through morphological identification in Cluster 1, while Euphausiidae and Mysidae were detected more in Cluster 1 with metabarcoding. Notilucaceae, Diphyidae, and Sagittidae, which were associated with the Hallyeo area, were more common in Cluster 2 than other mesozooplankton clusters using both methods. The phylum Rotifera included in other with Acartiidae, Podonidae, and Centropagidae, which were more dominant in the Yellow Sea, were identified more in Cluster 3 by metabarcoding. Figure 6. The association between spatial, environmental characteristics and mesozooplankton communities. CAP analysis for zooplankton communities based on Bray-Curtis dissimilarities according to each category and identification method: (a) spatial and environmental variables, and morphological identification; (b) spatial and environmental variables, and metabarcoding. The arrows on the CAP plots in (a) and (b) indicate the patterns in response to the spatial and environmental variables for the zooplankton community clusters. Bar plots between zooplankton community clusters formed using (c) morphological identification and (d) metabarcoding according to spatial and environmental variables in CAP analysis.
Depending on the metabarcoding results, the dominant or uniquely identified taxa were considered as potential bioindicator taxa that characterize the mesozooplankton cluster (Figure 6d). Paracalanidae, Diphyidae, and Sagittidae, and Noctilucaceae, which were common in Cluster 2, could be associated with high water temperature, salinity, and topography. Acartiidae, Podonidae, Rotifera, and Centropagidae, which were more dominant in Cluster 3, could be bioindicators for inflow of freshwater and chlorophyll a concentration.

Comparison between the Morphological Identification and Metabarcoding Results
In this study, the efficiency of metabarcoding was verified by comparing it with the results of morphological identification. Therefore, we could validate the use of metabarcoding for investigation of the mesozooplankton community of the Marine and Coastal National Park areas.
Consistent with the results of previous studies comparing the efficiency of morphological identification and metabarcoding, our results demonstrated that metabarcoding was able to detect much more zooplankton taxa than morphological identification. Additionally, mesozooplankton community structures clustered in similar patterns when the results of both methods were compared. The morphological identification method may overlook small-sized zooplankton species and premature or cryptic species that are difficult to distinguish morphologically. In contrast, the sensitive detection capability of metabarcoding is likely to detect small, immature, and cryptic individuals, which cannot be detected by the naked eye. In our study, many individuals of the phylum Rotifera, that were not morphologically identified, were detected by metabarcoding. There is less interest in Rotifera compared to other taxa and domestic taxonomic experts of Rotifera are rare.
In addition, through species identification using DNA barcoding, it was confirmed that there are many cryptic species in this phylum. As such, ecological studies of Rotifera have limitations [54]. However, they are important for understanding the aquatic ecosystem, as this taxon represents an important food source for large aquatic organisms such as crustaceans and fish [55]. These results reveal that metabarcoding may be more useful than morphological identification for the detection of Rotifera. Additionally, metabarcoding brought the presence of the taxon to our attention, so we will be more aware of Rotifera when morphologically examining zooplankton communities in the Marine and Coastal National Parks.
Our study also revealed some limitations of the metabarcoding method that were previously reported. Consistent with the results of previous studies, we found that the biodiversity and taxonomic composition of mesozooplankton communities were different between the morphological identification and metabarcoding methods. In particular, the abundance of Calanidae, which was relatively large compared to other taxa, tended to be overestimated by metabarcoding. Among the copepods collected from our studies, Calanidae individuals generally have a larger body size (up to 3 mm) than Acartiidae, Centropagidae, and Paracalanidae. The large body size of these organisms may contribute to the amount of DNA extracted from a sample, resulting in an overestimate [56][57][58]. The underestimated abundance of the dinoflagellate Noctiluca scintillans in metabarcoding appears to be due to the low efficiency of DNA extraction compared to other zooplankton taxa. The DNA extraction efficiency for dinoflagellates varies according to the protocol [59]. It is inferred that a relatively small amount of DNA was extracted from Noctiluca scintillans due to the use of a zooplankton-focused method of DNA extraction. These technical biases, including DNA extraction and PCR biases, distort the actual number of sequences [8,[60][61][62][63][64]. Additionally, Oikopleuridae, which was among the most frequently detected family levels, was not detected with metabarcoding. Considering the results of previous studies, which detected Oikopleuridae in the stomach of fish using the same primer [34,65], it is expected that the lack of detection of Oikopleuridae may have been caused by the technical biases generated during the sampling or experimental processes.

Potential Bioindicator Taxa in the Marine and Coastal National Park Areas of Korea in Spring
Zooplankton taxa can provide early detections of global climate change due to their sensitivity to environmental changes. Metabarcoding has a sensitive detection capability, which can identify potential indicator taxa in bulk samples or communities [66]. Using metabarcoding, we identified the characteristics of three clusters divided according to spatial and environmental variables. In addition, using the results from both identification methods, we found potential bioindicator taxa that were related to the characteristics of Cluster 2 and 3.
Cluster 1 was unable to identify any correlations between the cluster of mesozooplankton communities and the spatial and environmental variables. This indistinctness may be attributed to the diverse geographical characteristics and extensive range of the Dadohae area. The sampling stations in the Dadohae area are distributed in both the Yellow Sea and the Southern Sea of Korea; therefore, these sampling stations are affected by the environmental characteristics of both sea areas. In addition, due to the seasonal changes of the Kuroshio currents and the southward movement of freshwater in the Yellow Sea by wind, the zooplankton habitat here changes more frequently than in other areas [67,68]. Using metabarcoding, Calanidae, Euphausiidae, and Mysidae were found more in this cluster compared with others. However, we were able to identify the common characteristics of these taxa that reflect the characteristics of Cluster 1 in this study.
Cluster 2 was associated with longitude, water temperature, and salinity. This distinct clustering could be a result of the environmental characteristics of the Kuroshio Current and topographical characteristics of the Southern Sea of Korea. Paracalanidae, Diphyidae, and Sagittidae, detected in high abundance by metabarcoding, appear to be associated with high temperature and salinity, which are characteristics of the Kuroshio Current. The Kuroshio Current has relatively high temperature and salinity compared with other currents affecting the Korean Peninsula [69]. The Genus Paracalanus belonging to Paracalanidae is one of the common copepods on the coast of Korea, which are reportedly correlated with high water temperature or salinity [70][71][72][73]. Diphyidae can be easily moved through ocean surface currents and thrive explosively upon encountering a preferred environment [74,75]. Most jellyfish are known to prefer high water temperature and salinity in marine environments [76]. In addition, Chaetognatha, a phylum that includes Sagittidae, is moved by the Kuroshio Current and its distribution is closely related to the physical and environmental characteristics (e.g., high water temperature and salinity) of these currents [77][78][79]. Noctilucaceae was also detected more in Cluster 2 than others when metabarcoding was used, although not as much as the result of morphological confirmation. The hydrographical characteristics of the Hallyeo area and high salinity of the Kuroshio Current may also contribute to this result, as the most widely known of Noctilucaceae species, Noctiluca scintillans, is widely distributed globally and is one of the red tide forming species [80,81]. The distribution of Noctiluca scintillans in Cluster 2 appears to be affected by unique hydrographical characteristics (e.g., topography) in the Hallyeo area. The Hallyeo and part of the Dadohae areas in Cluster 2 are well developed, partially enclosed bays. This topography has the characteristic of accumulating buoyant cells of Noctiluca scintillans, causing large blooms [80]. In addition, previous studies reported that salinity is positively correlated with the number of Noctiluca scintillans individuals in Gwangyang Bay, a nearby sea area of Hallyeohaesang National Park. Thus, Noctiluca scintillans are likely well adapted to high salinity conditions [82,83].
The metabarcoding identification results revealed that the proportions of Acatiidae, Podonidae, Rotifera, and Centropagidae were found to be higher in Cluster 3 than in other clusters. This cluster consisted mostly of samples from the Taean and Byeonsan areas in the Yellow Sea, which is associated with the inflow of freshwater and high concentrations of chlorophyll a. The Taean area and Byeonsan area, in the Yellow Sea, have freshwater inflows from the Geum River, Mankyung River, and Dongjin River. In addition, these areas have constructed artificial seawalls to prevent the inflow of seawater to the land due to large tidal differences. To improve the water quality of the lake created by the artificial seawall, a large quantity of freshwater is released into the sea through floodgates. This inflow of freshwater appears to have created a habitat for coastal species of zooplankton that are adapted to the low level of salinity. This release of freshwater can cause a change in the zooplankton assemblage [84][85][86][87][88]. For example, as salinity decreases in the surrounding marine environment, high-salinity tolerant species are replaced by low-salinity tolerant species with similar functions in the marine ecosystems [89]. In our results, a high proportion of Acartiidae was found in Cluster 3 with the metabarcoding method. Using morphological identification, Acartiidae were identified to the species level as Acartia hongi, Acartia hudsonica, Acartia ohtsukai and Acartia omorii. Podonidae, which were abundant in the Taean area, were identified morphologically as Pleopis polyphemoides. This species has the characteristic of preferring brackish water and river estuary areas and is known as being highly resistant to low salinity [90][91][92][93][94]. The phylum Rotifera also consist of freshwater invertebrates that play a pivotal role in freshwater and marine ecosystems, as mentioned above [95]. With the inflow of freshwater, it can be inferred that the proportions of Acartiidae, Podonidae, and Rotifera, which prefer low salinity, were higher in Cluster 3 than in other mesozooplankton community clusters. The occurrence of a highly detected Centropagidae species appears to be closely related to the chlorophyll a concentration. As mentioned above, the average chlorophyll a concentration was higher in the Yellow Sea compared with that of the Southern Sea of Korea. A Centropagidae species detected using metabarcoding was identified as Centropages abdominalis and verified by morphological identification. Similar to our results, Kang and Kim [71] also found that the occurrence of Centropages abdominalis is positively related to the concentration of chlorophyll a, and the amount of phytoplankton greatly affects its growth and development.

Conclusions
In this study, a metabarcoding method was investigated and compared to the traditional morphological identification method to establish an efficient survey and research approach for mesozooplankton community analysis in the Marine and Coastal National Parks.
Consistent with previous studies, the use of metabarcoding resulted in the detection of more taxa than morphological identification and distinguished differences between mesozooplankton communities with less time and labor. In addition, it was possible to successfully detect taxa that are morphologically difficult to identify, such as the phylum Rotifera. However, the inconsistencies between the results of morphological identification and metabarcoding in biodiversity and taxonomic composition of the zooplankton community due to technological biases must be taken into consideration. The most important results were the discovery of potential bioindicator taxa that could represent the spatial and environmental characteristics of the Marine and Coastal National Parks with the sensitive detection capability of metabarcoding. It is expected that further investigation and research on these taxa will enable us to effectively detect environmental changes in these national parks.
For efficient surveys of the mesozooplankton community in the Marine and Coastal National Parks, we propose the use of metabarcoding as a tool for finding bioindicator taxa that represent environmental changes. It is necessary to continuously search for taxa in these national parks that are worth researching, such as Rotifera, through metabarcoding. In addition, the efficiency of metabarcoding should be validated for various conditions (e.g., smaller planktons, additional primer sets, and various seasons). Establishing a monitoring system using this approach will help to identify the mid-to long-term patterns of change in the zooplankton community and the changes in bioindicator taxa according to environmental changes, which will be an effective tool for managing the marine ecosystems in the Marine and Coastal National Parks.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/12/6/233/s1, Table S1: Statistical differences in mesozooplankton biodiversity for the sea area using morphological identification and metabarcoding, Table S2: Statistical differences in mesozooplankton biodiversity for the location using morphological identification and metabarcoding, Table S3: Statistical differences in mesozooplankton communities for the sea area and location between morphological identification and metabarcoding, Table S4: List of zooplankton species identified by morphological identification and metabarcoding, Table S5: Results of morphological identification, Table S6: Results of metabarcoding.