Integrative Description of Cryptic Tigriopus Species From Korea Using MALDI-TOF MS and DNA Barcoding

MALDI Time-of-Flight Mass Spectrometry (MALDI-TOF MS) provides a fast and reliable alternative method for species-level identification of pathogens and various metazoans. Compared to the commonly used mitochondrial cytochrome c oxidase subunit I (mtCOI) barcoding, advantages of MALDI-TOF MS are rapid species identifications and low costs. In this study, we used MALDI-TOF MS to determine whether spectra patterns of different species can be used for species identification. We obtained a total of 138 spectra from individual specimens of Tigriopus, which were subsequently used for various cluster analyses. Our findings revealed these spectra form three clear clusters with high AU value support. This study validates the viability of MALDI-TOF MS as a methodology for higher-resolution species identification, allowing detection of cryptic species of harpacticoida. In addition, we propose a new species, Tigriopus koreanus sp. nov. by utilizing integrative methods such as morphological comparison, mtCOI barcoding, and MALDI-TOF MS.


INTRODUCTION
Intertidal copepods of the genus Tigriopus Norman, 1869 are widely distributed in high-shore intertidal rockpools across several continents (Tsuboko-Ishii and Burton, 2018). The genus consists of 17 species and is known to be easily bred under laboratory conditions, as it tolerates abrupt changes in environmental factors such as temperature and salinity (Kwok and Leung, 2005). Tigriopus californicus (Baker, 1912), a well-known representative of the genus, shows extensive population divergence and is a model for understanding allopatric speciation (Barreto et al., 2018). Reproductive isolation and molecular information have confirmed cryptic speciation in T. californicus (Peterson et al., 2013), but it remains unresolved.
Similarly, Tigriopus, which can be found also in the coast of Korea, is represented by T. japonicus Mori, 1938 and has been the subject study of several contributions. However, it has long been recognized that the Korean populations of T. japonicus are molecularly different from those of Japan (e.g., Jung et al., 2006;Ki et al., 2009).
In Karanovic et al. (2018), it was shown that the populations of T. japonicus inhabiting Korean coasts are actually two different species: Tigriopus west Karanovic et al., 2018 and Tigriopus east Karanovic et al., 2018. These three species are morphologically similar, making it nearly impossible to separate them through optical microscopy. Their cryptic nature was unraveled by Karanovic et al. (2018) through molecular and morphological comparisons. For morphological comparison, they selected landmarks on somitic sensilla and segments of different appendages, which were used for two-dimensional geometric morphometrics. In addition, average Kimura 2parameter (K2P) pairwise distance of mtCOI partial sequences between the two Korean species were revealed to be 0.226, confirming that Tigriopus in the west-south coasts and east coasts of Korea are two different species.
MALDI-TOF MS (Matrix-Assisted Laser Desorption/ Ionization Time of Flight Mass Spectrometry) initially hits the sample with a pulsed laser to ionize it. This procedure is known as laser desorption ionization. Using a matrix-which is a small material that can absorb laser wavelengths-the sample collapses when the temperature of the sample-matrix rises rapidly by absorption of the laser. Fundamentally, MALDI is an ionization technique that uses a laser energy absorbing matrix to emit ions from large molecules with minimal fragmentation (Hillenkamp et al., 1991). The protein mass is analyzed with a TOF analyzer that measures the flight time of the ion emitted.
The representative protein spectra obtained by MALDI-TOF MS are used as a kind of barcode that can distinguish different species. This identification method was initially used for the identification of species of bacteria, fungi, and viruses (e.g., Claydon et al., 1996;Fenselau and Demirev, 2001), and is still widely used (e.g., Seng et al., 2010;Chalupová et al., 2014;Singhal et al., 2015). In addition, MALDI-TOF MS has been used for species identification of microalgae (Baumeister et al., 2020), microorganisms (Maier et al., 2006) and various metazoan taxa including insects (Feltens et al., 2010), fish (Volta et al., 2012;Shao and Bi, 2020), and calanoid copepods Laakmann et al., 2013;Bode et al., 2017). Previous studies with several taxa as study subjects confirmed that MALDI-TOF MS is a reliable and fast alternative method for species identification.
Recently, Rossel and Martínez Arbizu (2019) suggested that the identification of species through protein mass measurement is a faster and cheaper method compared to mtCOI barcoding. In general, barcoding takes at least a day or more to obtain the sequencing result following DNA extraction of the target organism, and additional hours to complete PCR, electrophoresis, and purification. Most importantly, barcoding method varies highly in success rates depending on the primer sets used and requires prior knowledge on the target organism to find a suitable primer set. On the other hand, the protein spectrum used for MALDI-TOF MS can be obtained immediately by incubating the individual in a matrix for several minutes after drying the sample and measuring it with a MALDI-TOF MS machine. Rossel and Martínez Arbizu (2019) argued that the results could be obtained relatively quickly, at a lower cost. They (Rossel and Martínez Arbizu, 2019) used MALDI-TOF for the assessment of 75 species and confirmed that cryptic species Leptastacus can be separated using this approach . In this study, we used MALDI-TOF MS to obtain spectra of Tigriopus from Korea and determine whether their varying peaks can be used to detect and differentiate cryptic species. In addition, a new species of Tigriopus is proposed based on approaches such as mtCOI barcoding, clustering of protein spectra, and morphological comparative analysis.

Collecting and Morphological Comparison
Specimens of Tigriopus inhabiting intertidal rock pools were collected with a spuit at seven different coastal areas of Korea (Figure 1 and Table 1) and brought to the laboratory in plastic container with seawater. In the absence of other animals, specimens of Tigriopus were kept alive in a 20 • C incubator under a photoperiod of 12 h light: 12 h dark, for 2-3 months after being brought to the laboratory. As food source for copepods, a small amount of Tetra Bits was provided.
For morphological observation, slides were prepared using the sandwich method of 20 adult females of each species (total 60 individuals). Figure works for morphological comparison of each appendage (Mxp, P1 enp-1, P1 enp-2, P4 exp-3, P5 exp) were completed through observation at 400x magnification with an Olympus BX51 differential interference contrast microscope equipped with a drawing tube. First, length of each segment was measured to calculate the length/width ratio. Based on the ratio, analyses were conducted to compare the differences between the mean of the three species in SPSS. Normal distribution was tested, and while some Mxp's data were below significance, because the number of samples was sufficient, the data were kept for analysis. In addition, Bonferroni correction was performed, and ANOVA test determined which features showed significant morphological differences between species (p < 0.05).

Scanning Electron Microscopy (SEM)
In total, 18 specimens from BS were prepared for Scanning electron microscopy (SEM). For SEM preparation, specimens were transferred into hexamethyldisilazane (HMDS) for drying (Shively and Miller, 2009), mounted on stubs, and finally coated with gold in an ion coater, SPT-20 (COXEM, Korea). Coated stub was observed under SEM on the in-lens detector at an accelerating voltage of 15.0 Kv. The specimens were photographed using a COXEM EM-30 SEM at the Biodiversity laboratory in Hanyang University. Digital photographs were processed and combined into plates using Adobe Photoshop CS6.

Pretreatment for MALDI-TOF MS
Total of 143 individual harpacticoids (138 Tigriopus and 5 Schizopera) were transferred to a 6-well plate a day before and fasted. On the day of the experiment, copepods were fixed in 99% ethyl alcohol for approximately 30 min, and then transferred to a 1.5 ml tube containing small amount of ethyl alcohol (0.5  µl). Once all ethyl alcohol in the tube evaporated, leaving the copepod in completely dried state, 4 µl of α-Cyano-4-hydroxycinnamic acid (HCCA) matrix [Acetonitrile 50%; U.P. Water 47.5%; Trifluoroacetic acid 2.5%; supersaturated HCCA (30 mg for total 1 ml matrix)] was added to each tube and incubated at room temperature for at least 20 min. Two microliter of this solution was placed on the target plate of the MALDI-TOF MS equipment (AXIMA Confidence MALDI TOF-Mass Spectrometer; Shimadzu) and once the solution dried completely to crystalize, the plate was placed in the instrument to be analyzed.
Protein mass spectra were measured between 2k to 20k Dalton on a MALDI-MS Application Launchpad 2.9.2 (Shimadzu Biotech) software using linear mode with laser power 80. To create a sum spectrum, 150 profiles were summed up. Each profile was measured 10 times, and the obtained individual protein spectra were exported to ASCII format files after range editing in Data Explorer (TM) software version 4.3.

Workflow in R
The ASCII files were imported to RStudio using the R package MALDIquantForeign (ver. 0.12) (Gibb, 2015) and the subsequent data processing and analysis of the data was achieved using MALDIquant (ver. 1.19.3) (Gibb and Strimmer, 2012) and MALDIrppa (ver. 1.0.5) (Palarea-Albaladejo et al., 2018). For general methods to obtain the featureMatrix mentioned below, we referred to Species Identification using MALDIquant manual 1 .
The range of all protein mass spectra were trimmed to 2,000-15,000 m/z. We conducted a basic quality control with raw data and tested whether all spectra contain the same number of data points and are not empty before overall analysis. Spectra were square root transformed and smoothed with the Savitzky-Golay method. The baseline was removed based on SNIP baseline estimation method and spectra were normalized using the Total-Ion-Currentcalibration (TIC) method implemented in MALDIquant. Peak detection was applied with a signal-to-noise ratio (SNR) of 8 and a half window size of 18. Peaks were repeatedly binned with the "binpeaks" command from MALDIquant with a tolerance of 0.002. Through this process a feature matrix was created, which was then Hellinger transformed for further analysis. Results of the Hellinger transformed matrix were visualized through clustering analysis, diagonal discriminant analysis (DDA), and non-metric multi-dimensional scaling analysis.
Cluster analysis was performed using 138 spectra of data obtained from adult individuals of Tigriopus. Dendrogram was produced with pvclust (ver. 2.2) R package (Suzuki and Shimodaira, 2006) using Ward's 2D clustering algorithm with Euclidean distances and 10,000 bootstrapping replications. Approximately unbiased (AU) p-value and bootstrap probability (BP) values were provided alongside the dendrogram, of which AU values of ≥ 85% were considered strong evidence of a cluster for the analysis.
Non-Metric Multidimensional Scaling Plot (NMDS) was generated with the vegan (ver. 2.5.6) R package (Oksanen et al., 2018) based on Bray-Curtis Dissimilarity distance with k = 2. To test multivariate homogeneity of group dispersion, ANOVA test was performed utilizing betadisper function provided with the vegan (ver. 2.5.6) package. PERMANOVA test was also performed using Adonis tool to test the fit of the data with 999 permutations.
Diagonal discriminant analysis (DDA) was performed using Hellinger transformed featurematrix with sda (Ahdesmäki and Strimmer, 2010) function to find peaks with the highest variation among three different species. Total of 196 variables, 138 observations and 3 classes (species) were used to compute t-score for feature ranking.

DNA Extraction and Amplification
For DNA extraction and amplification, specimens were transferred to a well of ultrapure water for 2 h to remove ethanol or seawater. Non-destructive DNA extraction was then carried out by utilizing worm lysis buffer (Williams et al., 1992). Specimens were transferred to microcentrifuge tubes containing 25 µl lysis buffer and placed in a Takara thermocycler (Takara, Otsu, Shiga, Japan) with the following settings: 65 • C for 15 min, 95 • C for 20 min and 15 • C for 2 min. After this, specimens were kept for morphological identification. Unpurified total DNA was kept at −20 • C for long-term storage.
MtCOI genes were amplified using PCR premix (LaboPass, COSMOgenetech, Korea) and 3 µl of DNA template with Cop-COI-2189 and jgLCO or LCO1490 primers (Bucklin et al., 2010). The amplification protocol consisted of an initial denaturation at 94 • C for 5 min followed by 40 cycles of denaturation at 94 • C for 1 min, annealing at 45 • C for 2 min, and extension at 72 • C for 3 min; this was followed by a final extension step at 72 • C for 10 min (Bucklin et al., 2010). Successful amplification was confirmed by electrophoresis on a 1% agarose gel. PCR products were sent to Bionics (Seoul, Korea) for purification and DNA sequencing. For sequencing, an ABI automatic capillary sequencer was used with the same set of primers used for amplification. All obtained sequences were visualized using Finch TV (ver. 1.4.0) 2 (Geospiza Inc., United States). The quality of each sequence was evaluated, and low-resolution peaks were checked by comparing forward and reverse strands. BLAST search confirmed the obtained sequences as copepods without contaminants. Sequence information from this study was deposited in the NCBI database (MW429840-MW429844).

Phylogenetic Analyses
Thirty two additional sequences were downloaded from GenBank and included in our analyses (Supplementary Table 1). Sequences were aligned with the ClustalW algorithm (Thompson et al., 1994) in MEGA version 7.0 (Kumar et al., 2016). Average pairwise distances were also computed in MEGA version 7.0 using the K2P model. Phylogenetic analyses were performed using Neighbor-Joining (NJ), Maximum parsimony (MP), Maximum Likelihood (ML), and Bayesian Inference (BI) approaches.

Protein Mass Spectrum Pattern
Spectra data obtained from a total of 138 Tigriopus and 5 Miraciidae specimens were used for the analysis ( Table 1). In the protein spectra obtained from raw data (Figure 2), the x-axis represents protein mass, and the y-axis is intensity, showing how much protein is present for each mass value corresponding to the x-axis. Similar spectra peak patterns were observed among species collected from the same region. Prominent peak values were also confirmed to be similar by region. Protein spectra pattern did not seem to differ by body parts (prosome/urosome), or by sex (Figure 2). Presence of eggs did not seem to effect ability to obtain spectra. Even specimens fixed in 99% ethyl alcohol stored at 4 • C for several months yielded spectra without a problem (Table 1, marked with "-fix"), but when only female egg sacs were used for the analysis, accurate protein spectra data could not be obtained.

Cluster Dendrogram With p-values (%)
Cluster analysis of 138 adult individuals of Tigriopus collected from six different regions in Korea are shown in Figure 3. Analysis using Hellinger transformed matrix confirm that species are divided into three clusters (Busan, T. east and T. west) with fairly strong support (AU p-values = 91, 84, and 92%, respectively). Of the two p-values, AU is considered to be superior to BP in terms of bias by using multiscale bootstrap resampling as opposed to ordinary bootstrap resampling (Suzuki and Shimodaira, 2006).
It is worth mentioning that clustering analysis was attempted with both raw spectra feature matrix as well as Hellinger transformed matrix. Results from the Hellinger transformed matrix yielded a dendrogram with higher AU and BP values. While optimal combination of data transformation has been documented by Rossel and Martínez Arbizu (2018), they recommended Hellinger transformed data in most cases. Other MALDI-TOF MS studies conducted on copepods also utilized Hellinger transformed datasets (Rossel and Martínez Arbizu, 2018;.
The indication of a third new species of the cryptic species Tigriopus spp. is further supported below via NMDS analysis.

Non-metric Multi-Dimensional Scaling (NMDS)
Non-Metric Multi-Dimensional Scaling (NMDS) analysis of all 138 spectra of Tigriopus used in the experiment are shown in Figure 4. This is the result of visualizing a two-dimensional space by measuring similarities between spectra.
In the West Sea cluster, spectra from Jeju, Boryeong, and Baengnyeong island are mixed. Based on this, all of them are T. west, although there is a considerable geographical difference between Baengnyeong island and Jeju. In addition, it was confirmed that only the Busan cluster was gathered separately above. Therefore, it was determined that Tigriopus in six regions in Korea were divided into three groups.
It is quite clear from the NMDS plot that three groups are visibly distinguishable from one another (homogeneity test, P < 0.05, PERMANOVA test, P = 0.001). Stress value of this plot was 0.1518 which is in usable range according to Clarke (1993), who suggested stress value < 0.20 to be usable. Considering that our dataset contains 138 samples of different variables, and since stress value increases with the number of samples and variables, the stress value of 0.15 is quite acceptable.

Diagonal Discriminant Analysis
In Figure 5, the numbers on the left show top 40 protein masses that differed greatly among species. Simply put, in the case of T. koreanus sp. nov., they have a relatively large amount of 2,986−2,988 m/z compared to the other species, followed by T. west with a relatively large amount of 2,958 m/z. Such inferring peaks with high variation among different species can thus be considered when delimitating Korean species of Tigriopus. Figure 6 illustrates variation in important peaks shown in Figure 5. It is evident individuals exhibit similar peak patterns for each species. The parts not marked in blue mean that they are the peaks below the SNR criteria and have lower intensity compared to other species.

MtCOI Distance
Partial mtCOI sequences were successfully obtained from five specimens of the three putative different species (JMJ; BN, BR, SGP; BS). As a result of calculating the average pairwise distance with the mtCOI sequences obtained from each cluster, the mtCOI distances within each species are shown in Table 2. In the case of T. californicus, the average intra-species distance was 0.222, which is considered to reflect the fact that it is a cryptic species complex.
Remarks -From measuring and comparing the morphological features of 20 females of each species, subtle morphological differences in the length:width ratio (l/w) of the Mxp basis, l/w ratio of the P1 enp-1, P4 exp-3, and P5 exp (Table 3 and Figures 9, 10) among the Korean species of Tigriopus were revealed.
T. koreanus sp. nov. and T. west showed significant differences in Mxp basis l/w (P < 0.001), P1 enp-1 l/w (P < 0.001). The difference between T. koreanus sp. nov. and T. east showed significant differences in all four morphological features (P < 0.001). The Mxp basis l/w of T. koreanus sp. nov. was about 1.8, which was stouter than the other two species (l/w = 2). For P1 enp-1 l/w, T. west is about 2.8, the longest, T. east is 2.6, and T. koreanus sp. nov. is 2.4, which is different in all three species. On the contrary, P4 exp-3 l/w showed the narrowest trend at T. east at about 2.9, followed by T. koreanus sp. nov. with 2.7 and T. west with 2.6. P5 exp l/w was similar in T. west and T. koreanus sp. nov., but both species showed significant differences from T. east.
According to DDA, compared to the other two species, T. koreanus sp. nov. showed strong intensity at 2,986−2,988, 3,066−3,067, and 3,026 m/z, and these peaks can be regarded as representative peaks of this species (Figure 5). In the case of T. west, ten representative peaks were identified, including 2,958, 3,036−3,037, 3,617−3,618, and 4,242. Tigriopus east showed distinct differences from the other two species in 22 mass, including 2,950, 3,310−3,312, and 3,287−3,288.
Etymology -Species name refers to the type locality (i.e., Republic of Korea).

MALDI-TOF MS
There was no difference in raw spectra (Figure 2A) when an individual specimen was cut in half by prosome and urosome, processed in separate tubes by body parts. The protein spectra of Schizopera sp. (Miraciidae) were also obtained in order to observe spectra variation from species not belonging to Tigriopus (Figure 2F). The protein spectra of Schizopera sp. were significantly different from Tigriopus spp. (Figures 2A-E). Figure 5 shows relative difference between peaks when the three species are compared together, and it is sufficient to find the difference between the three Korean species of Tigriopus. Until now, various factors which could affect difference in spectra such as fixative solution and storage period after fixation have been investigated (Rossel and Martínez Arbizu, 2018). Usefulness of MALDI-TOF MS in cryptic species identification such as bacteria, mosquitoes, fish and sandflies have also been documented (Müller et al., 2013;Dieme et al., 2014;Paauw et al., 2015;Chavy et al., 2019;Maasz et al., 2020). This study further expands on the usefulness of MALDI-TOF MS, confirming that higher-resolution species identification (identification of cryptic species) is possible for copepods as well. Despite its ability to identify cryptic species, there are several limitations. First, closely related species have similar protein composition, so the aspects of mass spectra patterns may appear similar. In such case, using only the overall pattern of mass spectra may lead to insignificant results. We expect that increasing the SNR to focus specifically on peaks with strong signal difference help detect cryptic species better. Just as DNA barcoding varies in resolution by different gene markers, and markers are chosen depending on the specific study purpose, parameters such as SNR, halfwindow, minfrequency for MALDI-TOF methodology must be adjusted accordingly.
Secondly, MADLI-TOF MS does not provide insight on the phylogenetic relationships between species. Because this method simply compares the overall protein pattern for species identification, protein molecules, composition, or function is not considered. Granted, if there is a difference in genetic distance, it is reflected on protein composition, resulting in mass spectra difference. However, this does not reflect proportional difference in the genetic distance. As such, to gain insight on phylogenetic relationship whilst discerning species complex, MALDI-TOF MS should be used complementally with DNA barcoding method.
Thirdly, it can be misleading to delimit species solely using MALDI-TOF without any prior knowledge about the species being investigated. In such case, cluster result obtained from the dendrogram may be difficult to interpret, since the number of clusters (or in this case, species) differ depending on the criterion setting of the y-axis. However, since this study presented three pieces of evidence (morphology, molecule, and protein) together, other evidences allowed us to present species classification criterion in MALDI-TOF. What this study can suggest is that Korean species of Tigriopus can be classified based on the height FIGURE 6 | PeakPatterns plot of 138 Tigriopus spectra. X-axis; peaks with the highest variation among three Tigriopus species, Y-axis; individual numbers were displayed in 4 intervals.
value of 3 in the dendrogram through protein mass comparison. Since criteria vary depending on the taxa, in order to use MALDI-TOF for species identification, a library must be established through further studies dealing with many taxa. To achieve this, role of taxonomists to confirm morphology while building on curated spectra, or molecular data would be essential. Finally, MALDI-TOF MS can be affected by other factors, such as environmental factors (Karger et al., 2019;Chac et al., 2020). Karger et al. (2019) suggest that subtle differences within mass spectra can occur depending on habitat type, season, and region. This means environmental factors can directly influence differences in clusters, resulting differentiation of species. To improve the aforementioned limitations, sufficient accumulation of data from various species and habitats is required.
As is the case with other methods, it is still difficult to accurately identify a species with MALDI-TOF MS alone, so it is necessary to accumulate data library by morphological taxonomists, and in the beginning, it will be necessary to verify the species identification through other methods.

MtCOI Distance
Based solely on this result (Supplementary Table 3), it could be confirmed that the Busan community differs from the West coast and Jeju communities by about 6%. According to a study FIGURE 7 | Molecular phylogenetic Maximum Likelihood (ML) and Bayesian Inference (BI) tree of Tigriopus species based on 36 mtCOI sequences. All positions containing gaps and missing data were eliminated. The first number along the branches represents ML, the second number represents BI bootstrap values, and * indicates sequence obtained from this study.
on copepod DNA barcoding (Baek et al., 2016), the mean interspecific distance of harpacticoida was 1.6%, which is similar to the results of the three Korean Tigriopus species obtained in this study (Table 2). Also, in other copepod orders, there were cases where the mtCOI sequences between congeners differed only by 0.26−1.1%. Therefore, the difference of 6% is considered to be one of the supporting evidences for the species classification. As such, it is evident from existing studies, that range of intraspecific and interspecific distance varies for each copepoda order. The fact that the range of distance always vary among taxa is an indication that classification of species using only fragment sequence of one molecular marker is with limitations. With DNA templates, researchers can obtain sequences of different lengths and locations depending on the primers used. Therefore, distinguishing cryptic species with just mtCOI barcode may prove to be an ambiguous task, and there is a possibility that they will result in the same species, although they are in fact different species. In accordance with this, in several studies (e.g., Huys et al., 2007;Thum and Harrison, 2009;Wyngaard et al., 2010;Blanco-Bercial et al., 2011;Marrone et al., 2013;Vakati et al., 2019) suggest that multigene analyses are required for more accurate results when identifying species through molecular information. Therefore, rather than using only one marker such as mtCOI for species identification, providing protein analysis alongside can strengthen the support for the overall results.
In a previous study, Rossel and Martínez Arbizu (2019) showed that Leptastacus cf. laticaudatus can be separated through MALDI TOF-MS experiments in three putative cryptic species. The mtCOI distances of the three cryptic species included in this L. cf. laticaudatus complex ranged from 15.32 to 31.15%. Although mtCOI distance between T. koreanus sp. nov. and T. west was less (6%) compared to distances reported by Rossel and Martínez Arbizu (2019), results from protein analysis showed that they were indeed different clusters.
Clusters observed from the protein data are different from clusters from phylogenetic trees based on mtCOI data. Therefore, it is difficult to extract evolutionary and systematic relationships based solely from these results. It is also challenging to interpret whether T. east and T. west clusters are closer to the T. koreanus sp. nov. clusters.
It is predicted that the reason that a Jeju individual had a high similarity to Busan community is because Busan and Jeju are relatively closer than other West Sea regions. In the same vein, when comparing the mtCOI data uploaded to the NCBI, there was only 4.5−6.3% difference from the data known as T. west on the South Sea coast (Yeosu, Sacheon, Tongyeong) or on the northeast coast of Jeju (Gimnyeong, Seongsan), which are adjacent to Busan.
Although the East sea community is relatively close to Busan in geographic distance, the molecular distance differs as much as the distance from the Japanese species. A probable explanation for this is that T. east in the East Sea region have been influenced by the North Korea Cold Current coming from the north, rather than the Kuroshio Current coming from the south. If samples of Tigriopus from South Sea coast and North Korea are obtained in the future, it would be possible to verify this hypothesis through protein spectra measurements. In the BI and ML tree (Figure 7), the topology supports this hypothesis. Because the origin of T. east is different, it appears as a clade different from the species in the region related to the Kuroshio Current.
On the other hand, according to the MP tree (Supplementary Figure 1) and NJ tree (Supplementary Figure 2), the Korean species form one clade (bootstrap; 23, 67%, respectively). This supports the possibility that individuals originating from the same species by the Kuroshio Current were geographically divided, resulting in speciation. Considering the geological aspect, according to a previous study (Knowlton and Weigt, 1998), it was estimated that millions of years would have passed per 1.5−2.6% of mtCOI mutation rates in crustaceans. According to this result, the estimated timing of speciation of T. east and T. west from their most recent common ancestor would be between 20 and 11.6 mya (Miocene).   The two trees obtained in this study were similar to the topologies of the trees suggested in previous studies. The topologies of MP tree and NJ tree in this study were similar to that of Karanovic et al. (2018), which also used mtCOI gene marker. In these cases, T. east and T. west formed the same clade. On the other hand, the topologies of the BI tree and ML tree in this study (Figure 7) were similar to those of the BI tree in Ki et al. (2009) and Karanovic et al. (2018). Although the markers used were different in Ki et al. (2009), based on ITS-5.8S rDNA, the topology was corresponding.
As can be seen from different topologies depending on the marker, method, and model, in fact, no one knows what evolutionary histories remain veiled. Naturally, there are only a few hypotheses that we can infer from these results.
However, the second hypothesis is difficult to support since Tigriopus is a fast-evolving genus, and Miocene may be too old as time for speciation. Therefore, in this study, topological results of the BI tree and ML tree were deemed more reliable, and we suggest that T. east had different origin from T. koreanus and T. west.
Nevertheless, it is evident that Tigriopus species from Busan is different from the previously known T. west and T. east. The mtCOI barcodes further support that the population of Tigriopus inhabiting Korea consists of at least three different species.

Morphology
When comparing the morphological differences, the microscopic morphological features that distinguish the new species and other two species of Tigriopus were found in the lengthwidth ratio of the P1 enp-1 and Mxp. Karanovic et al. (2018) found the differences between T. west and T. east by geometric morphometrics in P4 exp-3, female P5 exp, and male P2 enp. In this study, differences in P4 exp-3 l/w and P5 exp l/w of the two species were also confirmed, but when comparing the P5 exp l/w of the three species together, it was confirmed that it was difficult to clearly distinguish the three species by these morphological features as T. koreanus shows overlap, the mid-level characteristics of the other two species.
Geometric morphometrics is a good tool to find significant differences in morphological characteristics, but it requires a lot of time and effort. Therefore, based on the results of this study, we recommend MALDI-TOF MS as a simple and efficient method for species identification rather than these timeconsuming methods.
Furthermore, studies that are likely to be conducted in the future should include experiments to confirm the difference between nauplii and copepodites of the species of Tigriopus. It would also be meaningful to apply this method to other sibling species complexes such as T. californicus ( Table 2; mtCO1 distance within genus; 0.22, S.E.; 0.02) to clear its convoluted taxonomic status.

CONCLUSION
Various identification tools such as mtCOl barcoding and geometric morphometry were utilized alongside morphological analysis to test the usefulness of MALDI-TOF for the identification of species of Tigriopus. Based on our findings, we conclude that identification of cryptic species can be made simpler, faster, and more cost-effective with MALDI-TOF MS than with mtCOI barcoding.
Through MALDI-TOF MS experiments on harpacticoid copepods of Tigriopus, we showed that the separation of specimens, fixation, presence of egg sac, and gender did not have a fatal effect on the clustering results. In addition, it was confirmed that the populations of Tigriopus from Baengnyeong island, Boryeong, and Jeju are in fact Tigriopus west, although there is a considerable geographical difference between Baengnyeong island and Jeju.
Although it has been known as one species for a long time due to the difficulty of morphological distinction even with optical microscopes, results from protein analysis in this study supported that Tigriopus, which inhabits rock pools in Korea's coastal intertidal zone, actually consists of three species, and that the Busan population is a third distinct species different from the previously reported two species.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: NCBI (accession: MW429840−MW429844) and Supplementary Materials.