Liolophura species discrimination with geographical distribution patterns and their divergence and expansion history on the northwestern Pacific coast

The chiton Liolophura japonica (Lischke 1873) is distributed in intertidal areas of the northwestern Pacific. Using COI and 16S rRNA, we found three genetic lineages, suggesting separation into three different species. Population genetic analyses, the two distinct COI barcoding gaps albeit one barcoding gap in the 16S rRNA, and phylogenetic relationships with a congeneric species supported this finding. We described L. koreana, sp. nov. over ca. 33°24′ N (JJ), and L. sinensis, sp. nov. around ca. 27°02′–28°00′ N (ZJ). We confirmed that these can be morphologically distinguished by lateral and dorsal black spots on the tegmentum and the shape of spicules on the perinotum. We also discuss species divergence during the Plio-Pleistocene, demographic expansions following the last interglacial age in the Pleistocene, and augmentation of COI haplotype diversity during the Pleistocene. Our study sheds light on the potential for COI in examining marine invertebrate species discrimination and distribution in the northwestern Pacific.

The sampling localities of one southern Chinese (ZJ) and two Japanese (EH and MY) previously catalogued haplotype sequencing studies retrieved from NCBI are also depicted. Table S1 and S2 contain more accurate information on the populations and individuals. The basic map is from a free map providing site (https://d-maps. com), which is modified with Adobe Illustrator v. 25 Phylogenetic and population genetic analyses based on 16S rRNA. We constructed a nucleotide sequence alignment set with 34 16S rRNA haplotypes of L. japonica (Data S2), and identified 35 polymorphic sites (6.9%; Table S8) and 24 parsimoniously informative sites (4.7%). Phylogenetic analyses, including ML, BI, and NJ analyses, were conducted with the outgroup Acanthopleura echinata (Table S6). The resultant phylogenetic trees (Fig. S2) and unrooted phylogenetic network (Fig. 3a) consistently supported the three distinct genetic lineages of L. japonica, with the phylogenetic relationship between Lineages S1 and S2 being much closer than those inferred from the results of COI (Fig. 2a,b; Fig. S1). The TCS network (Fig. 3b) revealed that Lineages S1 and S2 were closely connected with only with 4-5 mutation steps between them, while Lineages N and S1 or Lineages N and S2 were remotely distanced by 18 mutation steps. Also, the overwhelming dominance of the RA1 haplotype implied a recent and rapid population expansion of Lineage N. In addition to RA1, haplotypes RB1 for Lineage S1 and RC1 and RC2 for Lineage S2 were dominant ( Fig. 3b; Table S7). Consistent with this, in the PCoA plot (Fig. 3c), the three genetic groups of L. japonica were spatially separated. Lineage N was distantly located apart from Lineages S1 and S2, while Lineages S1 and S2 were spatially much closer.
Examination of species discrimination of L. japonica based on COI and 16S rRNA. Using the Automatic Barcode Gap Discovery (ABGD), we performed distribution of pairwise genetic divergences, ranked pairwise difference, and automatic partition analyses based on COI and 16S rRNA of L. japonica, respectively ( Fig. 4a-c), which confirmed that there were distinct barcoding gaps between intraspecific and interspecific variations, strongly supporting the possibility of species discrimination of L. japonica. the COI-based analysis yielded two different barcoding gaps, while the 16S rRNA-based analysis revealed only a single barcoding gap ( Fig. 4a-c). The results of automatic partition at each value of the prior intraspecific divergence (P) divided L. japonica into three groups by COI and two groups by 16S rRNA, respectively (Fig. 4a-c). We also implemented two DNA taxonomy approaches to evaluate the possibility of species discrimination based on COI: the general mixed Yule coalescent (GMYC) approach (Fig. S3) and a Bayesian implementation of a Poisson Tree Processes model (bPTP) (Fig. S4). The results consistently and robustly supported the possibility that L. japonica can be divided into three different species, as shown in the results of ABDG (Fig. 4a-c). The molecular variance analyses using analysis of molecular variance (AMOVA), based on COI and 16S rRNA, were conducted to evaluate the degree of genetic differentiation among Lineages N, S1, and S2 (Tables S9, S10). According to the results, supposing that there are three genetic lineages (N, S1, and S2) or two genetic lineages (N and S1/S2), almost all variation in both cases is attributed to variation among groups (= among lineages), whereas variations within populations (within lineages) exhibit negative values in common. We confirmed that there was a high degree of genetic differentiation among Lineages N, S1, and S2, which supports the results of the COI barcoding gap analysis shown in Fig. 4a-c, although this was not statistically significant (P > 0.05; Table S9). When we assumed only two genetic groups, Lineages N and S1/S2, the genetic differentiation between the two www.nature.com/scientificreports/ groups was statistically significant (P < 0.001) in COI but not in 16S rRNA (P > 0.05) (Tables S9 and S10). The discrepancy between the number of barcoding gaps inferred from COI and 16S rRNA may have been affected by different gene evolutionary rates of the molecular markers 11 ; nucleotide substitution rate of 16S rRNA is known to be generally slower than that of COI (which is especially fast in the third codon positions: 105 out of 127 polymorphic sites). When an ML tree was constructed based on 22 polymorphic sites, which are found only in the first and second codon positions of COI that are much more conserved than the third codon position, the three genetic lineages were retained in the resultant tree ( Fig. S5), but Lineage S2 was nested within Lineage S1, as in the trees inferred from 16S rRNA (Fig. S2). Reflecting the powerful resolution of the COI barcoding marker well known from animals 12 and the high degree of variation among the three genetic lineages (Fig. 4, Figs. S3, S4), we suggested that L. japonica could be categorized into three different species: L. koreana, Yeo and Hwang, sp. nov. for Lineage N, L. japonica for Lineage S1, and L. sinensis Choi, Park, and Hwang, sp. nov. for Lineage S2. To examine whether it is reasonable to give these a species-level taxonomic status, as shown in Fig. 4d, we reconstructed a COI-based NJ tree with one congeneric species L. tenuispinosa 13 , which was originally described as a subspecies-level taxon of L. japonica 14,15 and was then revised as an independent species closely related to L. japonica by Saito & Yoshioka 16 in 1993. The resultant tree ( Fig. S6 and Data S3) showed that L. tenuispinosa forms a sister group with L. japonica (Lineage S1). This likely indicates that L. koreana and L. sinensis have taxonomic status as independent species.
Morphological comparison and geographical distribution of the three Liolophura species. We compared morphological characteristics among Liolophura koreana, sp. nov. (Lineage N), L. japonica (Lineage S1), and L. sinensis, sp. nov. (Lineage S2). Their morphological appearances are shown in Fig. 5a-c, which indicated that black spots on the tegmentum (Fig. 5d-e) and shapes of spicules on the perinotum (Figs. 5f-k, 6e-f) represent key morphological characteristics to distinguish them from each other. Although black dots in pleural areas, which are between the middle and lateral areas of the tegmentum on valves II-VII (or VIII), are commonly shared in all three lineages ( Fig. 5a-c), other black spots on the valves exhibit a high degree of variation in morphology ( Fig. 5a-c, Fig. S7). Herein, we described a new species of genus Liolophura, that is, L. koreana Yeo and Hwang from South Korea and Japan, with detailed descriptions of morphological characteristics observed by light microscopy (M205, Leica Camera AG, Germany) and FE-SEM (SU8220, Hitachi, Japan). In addition, we suggested the divergence of a new species, L. sinensis Choi, Park, and Hwang from southern China, with simple remarks based on distinct genetic difference (mainly COI barcoding gaps), with possible unique morphological characteristics as follows. Description. Body small-sized and broad oval-to oval-shaped ( Fig. 5a; Fig. S7  showing the three different genetic lineages for L. japonica: Lineage N members are most likely from the populations inhabiting a wide range of South Korea and Japan, Lineage S1 members from the populations inhabiting southern coastal areas of South Korea and Japan only, and Lineage S2 members from the southern Chinese population. As shown in Fig. S1, Acanthopleura spinosa was used as an outgroup. Numbers on branches indicate node confidence values: BP in ML, BPP in BI, and BP in NJ in order. (b) A phylogenetic network reconstructed using the neighbor net algorithm without an outgroup, showing three different genetic lineages for L. japonica inhabiting the northwestern Pacific coast: Lineages N, S1, and S2. The COI sequence alignment set used is shown in Data S1. Detailed information of the 106 COI haplotypes used in this phylogenetic analysis is summarized in Table S1 and S2. (c) An unrooted TCS network showing three distinct genetic clusters, corresponding to Lineages N, S1, and S2. Three different genetic groups correspond to the three genetic lineages shown in the phylogenetic tree (a), respectively. The haplotype frequency is displayed by the circle size. (d) A two-dimensional PCoA plot showing the three distinct genetic groups corresponding to Lineages N, S1, and S2 shown in the phylogenetic tree (a). The score on the first two axes (Axis 1 = 79.05% and Axis 2 = 15.32%) from the matrix of genetic distances estimated with the 106 COI haplotypes are indicated.  www.nature.com/scientificreports/ line anteriorly or white line laterally), with black dots on the pleural areas of valves II-VII (or VIII) (Fig. 5a,d, Fig. S7); articulamentum entirely black (dark brown); whitish and blackish spicules on the perinotum scattered irregularly, sometimes forming a band besides each valve (Fig. 5a, Fig. S7). Surface of the tegmentum in middle and lateral areas as in Fig. 6a,b and Fig. S8; posterior margin of the head valve nearly straight; dorsal shape of intermediate valves round-backed and side slopes slightly convex; the posterior valve margin with a distinct central apex, its shape subtriangular to triangular (rounded or linear), particularly valve II, mainly with a strong projection (Fig. 6a). Perinotum covered with large, solid, slightly curved, and obtusely pointed spicules (rarely with smooth and radial ribbed spicules apically), its density relatively lower than that of L. japonica (Figs. 5f,h,i, 6e).
Habitat. This new species appears to be attached to rocks in coastal areas with strong waves, or a calm inner shore in the northwestern Pacific Ocean (Fig. S9).
Etymology. The species is named per the locality of the new species.
Remarks. We found that Liolophura koreana, sp. nov. has no black spots on lateral areas of the tegmentum (Fig. 5d), and large, slightly curved, and obtusely pointed spicules on the perinotum compared to those of L. japonica (Figs. 5f,h,i, 6e). On the other hand, L. japonica from southern South Korea and southern Japan has two black spots on the lateral areas of valves II-VII (or VIII) (Fig. 5e), and small, almost straight, and cylindrical spicules compared to those of L. koreana (Figs. 5g,j,k, 6f). As shown in Fig. 7, L. koreana (Lineage N) was observed in all the South Korean and Japanese populations examined here, except for the EH population in Japan (refer to Tables S1, S2), which were found from JJ at 33°24′ N to MY at 38°32′ N. On the other hand, L. japonica (Lineage S1) was found from the southern coastal areas of South Korea and Japan, which were found only between JJ at 33°24′ N and TT at 35°53′ N. Interestingly, we found only L. koreana north from the latitude of 35°10′ N (BS) such as UL/DD (37°24′ N) and PH (36°02′ N) in South Korea. In Japan, there was found only L. koreana at MY (38°32′ N) too, but it remains to be explored to clarify its distribution range in Japan with much more sample collections covering northern Japanese coastal areas through further study. It was also confirmed that L. koreana and L. japonica show a sympatric distribution pattern between JJ at 33°24′ N and TT at 35°53′ N on the southern coastal area of South Korea and Japan.
Habitat. This new species appears to be attached to rocks in coastal areas of the northwestern Pacific Ocean.
Etymology. The species is named after its locality.
Remarks. L. sinensis, sp. nov. was examined and established mainly by molecular data such as the COI barcoding gap (Fig. 4a-c) presented in this study and photos provided from Prof. Yong-Pu Zhang (Wenzhou University, Zhejiang Province, China) without any direct real sample observation. The morphology of this new species is Phylogenetic network reconstructed using the neighbor net algorithm, showing three different genetic lineages for L. japonica: Lineage N, Lineage S1, and Lineage S2. The 16S rRNA sequence alignment set used is shown in Data S2. Detailed information of 34 16S rRNA haplotypes used in these analyses is summarized in Table S5 and S6. (b) An unrooted TCS network. There are distinctly observed three different genetic groups, corresponding to the three genetic lineages shown in the phylogenetic network (a). The haplotype frequency is displayed by the circle size. (c) A two-dimensional PCoA plot showing the three distinct genetic groups, corresponding to Lineage N, Lineage S1, and Lineage S2. The score on the first two axes (Axis 1 = 87.77% and Axis 2 = 4.4%) from the matrix of genetic distances estimated with the 34 16S rRNA haplotypes are indicated. www.nature.com/scientificreports/ very similar to that of a previously known species, L. japonica. L. sinensis from southern China has black spots on the lateral areas of valves II-VII similar to L. japonica, but black bowtie-shaped spots anterior to valves II-VII (Fig. 5c) are a unique characteristic for L. sinensis. L. sinensis (Lineage S2) was found in Zhejiang Province (ZJ) in southern China, around ca. 27°02′-28°00′ N ( Fig. 7; Tables S1, S5).

Demographic history and divergence time estimation analyses. Mismatch distribution analyses
(MDA) based on COI were performed for L. koreana, L. japonica, and L. sinensis, respectively. The MDA results (Fig. 8a) showed a unimodal curve for each of the three lineages. In addition, when neutrality tests were performed with COI and 16S rRNA (Table S11) A7  A4  A3  A8  A10  A14  A17  A2  A43  A6  A34  A36  A21  A37  A9  A5  A47  A42  A33  A15  A35  A40  A32  A22  A31  A11  A20  A38  A1  A46  A25  A28  A44  A45  A12  A26  A39  A30  A29    Intriguingly, we also found that there are ecologically subtle differences in the microhabitats of L. koreana and L. japonica: L. japonica appears to prefer rock surfaces found only inside the shoreline, where there are few waves, whereas L. koreana appears to be attached onto rock surfaces regardless of whether it is in an area with strong ocean waves or a calm inner shore (Fig. S9). On the eastern coastline of the Korean Peninsula (PH) and volcanic islands in East Sea (UL and DD), rocks may be directly exposed to incoming strong ocean waves without any blocking or buffering due to simple coastline. In these sites, only L. koreana were found. On the contrary,  www.nature.com/scientificreports/ the southern coastline of the Korean Peninsula (TY and GJ) and the islands located in South Sea of the Korean Peninsula (JJ and TS) have topographically complex coastline since it is ria coast, possibly providing various types of habitat conditions (i.e., strong sea waves can be weaker and seawater can be warmer). In these areas, both L. koreana and L. japonica may be found. The difference in geographic distribution of the two Liolophura species can be interpreted, at least in part, based on these physical differences in shorelines of northwestern Pacific Ocean. In addition, the southern coasts of the Korean Peninsula and the Japanese Archipelago are constantly under the direct influence of these warm currents (Fig. 7). In contrast, the seawater temperature on the east coast of the Korean Peninsula does not rise as high as the south ria coast, due to the influence of the cold Liman current that flows southward from Vladivostok, Russia, even though some of these warm currents travel northward through the East Sea (Japanese Sea) between the Korean Peninsula and the Japanese Archipelago. The difference in the seawater temperature may be one of the reasons why L. japonica is not observed on the east coast of the Korean Peninsula.
Recently, some studies have reported the delimited distribution of northern and southern types of marine invertebrates inhabiting coastal areas of the northwestern Pacific Ocean 17,18 . For instance, Xu et al. 18 reported that the chiton Acanthochitona rubrolineata (Lischke, 1873) can be divided into the northern type found north of Lianyungang (LYG; ca. 34°59′ N) and the southern type found south of Zhousha (ZS; ca. 29°98′ N) on the East Chinese coast based on COI, 16S rRNA, and 28S rRNA sequences. In addition, Cheng and Sha 17 Table S1-S3 contain the full names of localities and detailed haplotype information. The question mark indicates that collection of Liolophura samples from such coastal areas in Japan is required to clarify distribution patterns of L. koreana and L. japonica in the East Sea (= Sea of Japan). The basic map was obtained from a free map-providing site (https://d-maps. com), which was modified using Adobe Illustrator v. 25  www.nature.com/scientificreports/ sympatric distribution zone. On the Japanese coast, this species is divided into the northern type found north of the Araike Sea (AS; ca. 32°98′ N), and the southern type found south of Aomori (AO; ca. 40°82′ N), exhibiting a sympatric distribution zone. These recent studies, including our work, show that marine invertebrates exhibit distinct geographical distribution patterns (i.e. northern and southern types), which can be used for monitoring marine invertebrates whose migration north is fostered by global warming 19 . Further exploration of the detailed distribution patterns of the Liolophura species in the coastal areas of the northwestern Pacific uncovered in this study is required, with additional information on how long and how far the planktonic embryos and ciliated trochophore larvae float in the sea [20][21][22][23][24][25] .
Presumably, population expansions of the Liolophura species (Fig. 8b) might arise as a result of the last interglacial stage (129-116 Ka), called the Eemian interglacial period, named after the Eem River in the Netherlands, which was the last time when sea level was as high as or even higher than present-day sea level, and a time when the earth was slightly warmer than the present 26 . In addition, the molecular clock analysis (Figs. 8c, S10) estimated that L. japonica and L. koreana had diverged around the mid-Pliocene warm period 27 . During this period, the global average temperature was 2-3 °C higher 28 Fig. 8c) period. This approximation corresponds to the historical isolation of the East Sea (Sea of Japan) from the South and East China Seas, which was caused by a decline in the sea level during glacial cycles. Contemporary distributions of L. koreana (northern lineage) and L. japonica (southern lineage) are likely affected by two different surface water temperature zones in the northwestern Pacific 31 . In particular, after divergence, L. koreana and L. japonica may have adapted to the two aforementioned different surface water temperature zones. Subsequent species delineation of L. sinensis may have then occurred from the southern lineage (L. japonica) during the Pliocene to the late Pleistocene (3.19-0.96 Ma) period. Taken all together, L. koreana, L. japonica, and L. sinensis can be robustly separated with strong multi-pronged phylogenetic, population genetic, demographic, morphological, and ecological evidences. If several more samples are obtained from a wider range of coastal areas of the northwestern Pacific Ocean (particularly southern China and northern Japan), supporting the separation of L. japonica into three different species with additional evidence will be possible. The findings of the present study can be a stepping stone to developing a powerful indicator, such as a Rosetta Stone, for monitoring the northward migration of marine invertebrates as a result of global warming in the northwestern Pacific Ocean.

Materials and methods
Sample collection. A total of 342 L. japonica samples were collected from 12 locations in the intertidal coasts on the Korean Peninsula and Japanese Archipelago in the northwestern Pacific Ocean ( Fig. 1; Table S1). After collection, the entire chiton body was immediately fixed in 100% ethanol and stored at -20 °C until total cellular DNA extraction. PCR amplification and sequencing. Genomic DNA was extracted from the foot muscle tissue using an Exgene Tissue SV kit (GeneAll, South Korea) following the manufacturer's protocol. A 635-bp fragment of mitochondrial COI was amplified using the previously reported primer sets LCO1490 (5′-GGT CAA CAA ATC ATA AAG ATA TTG G-3′) and HCO2198 (5′-TAA ACT TCA GGG TGA CCA AAA AAT CA-3′) 32 . These primers, which were developed for use in invertebrates, can be used for PCR amplification of a partial mitochondrial COI gene fragment. A 506-bp sequence of the mitochondrial 16S rRNA was also amplified using the previously reported primer sets 16Sa (5′-CGC CTG TTT ATC AAA AAC AT-3′) and 16Sb (5′-CTC CGG TTT GAA CTC AGA TCA-3′) 33 . PCR amplifications were carried out in a 50-μl volume reaction containing 3 μl genomic DNA, 5 μl 10 × PCR buffer, 4 μl dNTPs mixture (2.5 mM), 1 μl of each primer (20 μM), 1 μl Taq DNA polymerase (5 U/μl) and 35 μl double-distilled water. Thermocycling consisted of an initial denaturing step of 2 min at 94 °C, 30 cycles of denaturation for 30 s at 94 °C, annealing for 30 s at 55 °C to 60 °C, and extension for 30 s at 72 °C, and a final extension step of 5 min at 72 °C. The amplified products were purified using HiYield™ GEL/PCR DNA Extraction Kit (RBC Co., South Korea), and then sequenced by the commercial sequencing service company Genotech Co. (Daejeon, South Korea) using an ABI PRISM BigDye terminator system and an ABI 3700  Genetic diversity and population genetic analyses. After trimming the obtained sequences using BioEdit 34 , to perform genetic diversity and population genetic analyses, 106 COI and 34 16S rRNA haplotype sequences were aligned with Clustal X 35 (Data S1, S2). The numbers of polymorphic sites and haplotypes, haplotype diversity, and nucleotide diversity (Tables S1, S4, S5, S8) were estimated for each population using the program DnaSP v.6.0 36 . AMOVA was conducted using ARLEQUIN v.3.5 37 to partition the genetic variance within and among populations. Furthermore, a statistical parsimony haplotype network was constructed with a 95% connection limit using TCS v.1.2.1 38 , and was used to assess the genealogical relationship by constructing the COI and 16S rRNA haplotype networks. To evaluate and visualize the geographic genetic structure among populations, PCoA based on the COI and 16S rRNA haplotypes were conducted using DARwin v.6.0.9 39 , which ordinated genetic distance estimates calculated using the haplotype data from this study.
Phylogenetic analyses. Phylogenetic analyses using the COI and 16S rRNA sequences based on the sequence alignment sets of 106 COI haplotypes (Data S1) and 34 16S rRNA haplotypes (Data S2) obtained in this study were performed using three different tree reconstruction algorithms: ML, BI, and NJ methods. In the ML tree, model selection in the IQ-Tree software package (http:// www. iqtree. org) tested and selected the substitution model HKY + F + I as the best-fit model under Bayesian information criterion. The tree was computed from 1000 ultrafast bootstrap replicates using the IQ-Tree website (http:// iqtree. cibiv. univie. ac. at) 40 . For BI, the models of sequence evolution for each gene were selected with jModeltest Ver. 2 41 , and were then trained in MrBayes (HKY + F + I). Each dataset was run for 10,000,000 iterations with a sample frequency of 1,000 iterations. After determining that the Markov Chain Monte Carlo (MCMC) runs reached a stationary level, the first 25% of the iterations ("burn-in" fraction) were discarded. The NJ tree was reconstructed using MEGA X 42 , and the confidence value for each node was calculated from 1000 bootstrap replicates. The haplotype sequences used for the phylogenetic analyses are listed in Tables S2 and S6, which include the sequences collected from the Zhejiang Province (ZJ) samples in southern China 8 and Miyagi (MY) and Ehime (EH) samples from Japan retrieved from the NCBI GenBank. The outgroup species employed for the phylogenetic analyses were A. spinosa for the COI dataset and A. echinata for the 16S rRNA dataset. Additionally, a phylogenetic network for each sequence alignment dataset was generated via the neighbor net algorithm 43 .
DNA barcoding gap analyses. The analyses of barcoding gaps based on COI and 16S rRNA were conducted using the online version of ABGD 44 (https:// bioin fo. mnhn. fr/ abi/ public/ abgd/ abgdw eb. html) to generate distance histograms, distance ranks, and automatic partitions. These analyses were conducted using the Kimura 2-P distance matrix 45 , and two different parameters: the range of prior intraspecific divergence from P min (0.001) to P max (0.1), and relative gap width (X = 1.5). Furthermore, we implemented two DNA taxonomy approaches to evaluate for the presence of cryptic species on the basis of COI: the GMYC approach 46 and a bPTP 47 . The GMYC approach was applied to an ultrametric tree produced by BEAST 2.6.0 48 with the Splits package (http:// splits. r-forge.r-proje ct. org). It is a process-based approach to detect the threshold at which processes within a species (i.e., coalescence) shift to the processes between species (i.e., speciation and extinction) in a gene tree. The bPTP was also performed to infer putative species boundaries on a given phylogenetic input tree. Unlike the GMYC model, the bPTP model requires a bifurcated phylogenetic tree rather than an ultrametric tree, and it can model speciation or branching events in terms of the number of substitutions. We used the following parameters for bPTP: MCMC 500,000 generations, 100 thinning, 10% initial iterations burn-in, and assessed convergence in each case to ensure the reliability of the results. www.nature.com/scientificreports/ common ancestor using BEAST 2.6.0 47 . The HKY model was selected as the best-fit nucleotide substitution model by jModelTest 41 and mutation rates of 2.0 × 10 −853,54 under a strict molecular clock 55 was used and MCMC was run for 30 million steps. Then, 10% of the iterations that had not reached convergence were burned, and a tree was sampled every 1,000 iterations. TRACER 1.7 program 56 was used to construct the BSP 57 . Divergence time estimation for the three phylogenetic lineages of L. koreana, sp. nov. (Lineage N), L. japonica (Lineage S1), and L. sinensis, sp. nov. (Lineage S2), as done in Choi et al. 58 , was conducted based on COI haplotype sequences using BEAST 2.6.0. 48 . The divergence time was estimated using the strict molecular clock 55 algorithm under the calibrated-Yule tree, which was calibrated using the earliest fossil record data of Mopalia (ca. 17.2 − 15 Ma) 59 for the basal node (17 Ma; normal distribution); a "monophyly" option was chosen in the BEAUti 2 program (Mopalia, COI haplotypes of the Liolophura species). Posterior distributions of the parameter were estimated using 1,000,000 MCMC iterations and sampled every 1,000 iteration, after discarding the initial 20% of iterations as burn-in. We used HKY as the best-fit substitution model selected by jModelTest ver. 2 41 . Additionally, an effective population size was determined using Tracer 1.7 56 . TreeAnnotator 2.6.0. 60 was used to produce a tree with maximum clade credibility and a median height after removing the initial 25% of iterations as burn-in. FigTree 1.4.2. 61 was used to visualize the topology of the resultant consensus tree.

Data availability
The sequence data have been deposited to the NCBI GenBank database under the accession numbers KT932836 − KT932912 for COI and KT932913 − KT932935 for 16S rRNA. The pre-processed datasets are available at Dryad Digital Repository, https:// doi. org/ 10. 5061/ dryad. rr4xg xd7h, which are named as follows: Data S1, Nucleotide sequence alignment of 106 COI haplotypes of L. japonica with the outgroup A. spinose; Data S2, Nucleotide sequence alignment of 34 16S rRNA haplotypes of L. japonica with the outgroup A. echinata; Data S3, Nucleotide sequence alignment of 106 COI haplotypes of L. japonica, L. koreana sp. nov., and L. sinensis sp. nov., one COI haplotype of L. tenuispinosa, and 14 COI haplotypes of eight Acanthopluera congeneric species with an outgroup Tonicia forbesii.