In total, 282 Ligia species specimens were collected from 16 sites in three countries, including 10 sites in South Korea, three in Japan, and three in Vietnam ( Figure 1 , Supplementary Table 1 ). All collected samples were directly fixed in 95% ethyl alcohol and stored at -20°C before DNA extraction. Genomic DNA was isolated using walking legs of individuals using a DNeasy Blood and Tissue Kit (Qiagen, United States). The concentration of extracted DNA was measured using NanoDrop 2000 (Thermo Fisher Scientific, United States) and the integrity was checked with a 1% agarose gel electrophoresis.

Species identification was performed done using mitochondrial 16S rRNA genes by DNA Sanger sequencing, with a 16Sa (5′-CGCCTGTTTATCAAAAACAT-3′) and 16Sb (5′-CTCCGGTTTGAACTCAGATCA-3) primer pair (Xiong and Kocher, 1991; Simon et al., 1994). PCR was conducted to amplify the partial fragment of the 16S rRNA gene with a 50 μL total volume of PCR mix comprising 25–100 ng genomic DNA, 10 mM dNTP, and 10 pM of forward and reverse primers with the addition of 0.25 units of Taq DNA polymerase (Solgent, South Korea). The following thermal cycling conditions were adapted as follows: denaturation at 94°C for 2 min, followed by 35 amplification cycles of 94 °C for 20 s, 48°C– 50°C for 40 s, and 72°C for 1 min, and a final extension at 72°C for 5 min. Subsequently, 1 μL of each PCR product was electrophoresed on a 1% agarose gel containing eco-dye and observed under UV light. The confirmed PCR bands were purified using a QIAquick PCR Purification Kit (QIAGEN, Germany) and directly sequenced with an ABI Prism 3730 DNA sequencer (PerkinElmer, United States) using a Big Dye Termination Sequencing Kit (PerkinElmer, United States).

In addition to 16S rRNA, NaK was amplified by PCR and sequenced from 11 Ligia individuals (9 L. exotica and 2 L. cinerascens) representing the 11 localities examined here ( Supplementary Table 2 ). Previously described primers Nak-for-b (5’-ATGACAGTTGCTCATATGTGGTT-3’) and Nak-rev2 (5’-ATAGGGTGATCTCCAGTRACCAT-3’) were used for PCR amplification of Nak (Tsang et al., 2008).

The obtained 16S rRNA and Nak nucleotide sequences of the two Ligia species were aligned using the Clustal X2 program by default setting, and ambiguous regions where homology could not be confidently identified at the 5’ and 3’ ends in sequences were removed by BioEdit 7.2.5 with a final length of 379bp of 16S rRNA and 663bp of Nak (Hall, 1999; Larkin et al., 2007). Haplotypes, singleton sites, and parsimoniously informative sites were identified using DnaSP 6.11 (Rozas et al., 2017). All novel 16S rRNA haplotype and Nak sequences discovered in the present study were deposited in NCBI GenBank with accession numbers OQ916835-OQ916875 of 16S rRNA and OR641242-OR641252 of Nak. Haplotype diversity (h) and nucleotide diversity (π) per population were calculated using DnaSP 6.11. To investigate the genetic differentiation between populations of each species, pairwise fixation index (F _ST) values were estimated using the Arlequin 3.5 (Excoffier and Lischer, 2010). Mantel test was performed using Alleles in Space 1.0 to investigate the patterns of isolation by distance between pairwise genetic and geographic distances (Miller, 2005).

For population genetic analyses, 197 previously published 16S rRNA sequences for Ligia species from the National Center for Biotechnology Information (NCBI) GenBank were included in the dataset. These included one sequence described by Wetzer (2001), one by Taiti et al. (2003), 39 by Jung et al. (2008), 24 by Yin et al. (2013), one by Santamaria et al. (2013), 42 by Hurtado et al. (2018), and 13 by Kang and Jung (2020), as well as 76 sequences found in unpublished data ( Supplementary Table 3 ). In the case of Nak, 11 previously reported sequences for 6 L. exotica and 5 L. cinerascens were retrieved from the NCBI GenBank for phylogenetic analyses ( Supplementary Table 2 ). In total, 22 Nak sequences, consisting of 11 obtained in the present study and 11 retrieved from the NCBI GenBank, were included for phylogenetic analyses. To examine phylogenetic relationships among the populations in the two Ligia species, unrooted phylogenetic trees were reconstructed using the maximum likelihood (ML) and Bayesian Inference (BI) methods. The ML tree was constructed using the IQ-tree web server (http://iqtree.cibiv.univie.ac.at/) (Trifinopoulos et al., 2016). The selection of the best evolutionary model was performed in the IQ-tree software package (http://www.iqtree.org). The TIM2 + F + I + G4 model was selected as a best-fit model based on Bayesian Information Criterion (BIC). The number of bootstrap iterations was 1000. For BI tree reconstruction, GTR + F + I + G4 model was selected as a best-fit substitution model using jModeltest 2.2 (Darriba et al., 2012), then used to reconstruct the BI tree in MrBayes 3.2 (Ronquist et al., 2012) with 10,000,000 generations, in which sampling was conducted every 1,000 iterations. The first 15% of the iterations were discarded after Markov Chain Monte Carlo (MCMC) runs reached a stationary level. Three additional phylogenetic analyses were conducted: an unrooted phylogenetic network based on the neighbor net algorithm was reconstructed using Splitstree (Huson and Bryant, 2006), whilst haplotype network based on the TCS approach were obtained using TCS 1.21 and popART aiming to investigate genealogical relationships between haplotypes (Clement et al., 2000; Leigh and Bryant, 2015). Principal Coordinate Analysis (PCoA) was used to visualize differences among sequences in Darwin 6.0.9 program (Perrier and Jacquemoud-Collet, 2006).

To evaluate species boundaries and number of species in the dataset, four different Molecular Species Delimitation Analyses (MSDAs) were performed: Automatic Barcode Gap Discovery (ABGD), Assemble Species by Automatic Partitioning (ASAP), Poisson Tree Process (PTP), and Bayesian implementation of PTP (bPTP) (Puillandre et al., 2012; Zhang et al., 2013; Puillandre et al., 2021). The two former methods are based on genetic distance, and the two latter methods are based on phylogeny. ABGD and ASAP analyses were performed based on Kimura two-parameter distance matrix (Kimura, 1980) using the online servers (ABGD: https://bioinfo.mnhn.fr/abi/public/abgd/; ASAP: https://bioinfo.mnhn.fr/abi/public/asap/). In ABGD, two more parameters were considered including P values (P _min value 0.001, P _max value 0.1) and a relative gap width of 1.5. In ASAP, default values of all parameters were selected. The PTP and bPTP analyses were carried out using an online server (https://species.h-its.org/). In the PTP and bPTP analyses, rooted phylogenies were reconstructed using the ML algorithm implemented in the IQ-tree web server (http://iqtree.cibiv.univie.ac.at/). Detailed parameters were as follows: 100,000 MCMC generations, a burn-in of 0.1 and a thinning value of 100. Genetic distance within or among the candidate species inferred from MSDAs patterns were calculated based on Kimura 2-parameter methods using MEGA 11 (Tamura et al., 2021).

To examine the historical demographic trends of L. cinerascens and two L. exotica lineages, we performed three different approaches based on 16S rRNA sequences: 1) three neutrality tests including Tajima’s D, Fu’s Fs, and R ₂ in DnaSP 6.11 and Arlequin 3.5, to investigate changes in population size in the past (Tajima, 1989; Fu, 1997; Ramos-Onsins and Rozas, 2002), 2) a mismatch distribution analysis (MMD) to estimate the frequency of pairwise differences using DnaSP 6.11 (Rogers and Harpending, 1992), and 3) a Bayesian skyline plot (BSP) to evaluate fluctuation of effective population size since the appearance time of the most recent common ancestor using BEAST 2.6.0 (Drummond and Rambaut, 2007; Bouckaert et al., 2014). In MMD, raggedness index (r) and sum of square deviations (SSD) were calculated to investigate if the dataset was suitable for the population growth-decline model using Arlequin 3.5 (Harpending, 1994). In BSP, mutation rate was directly calculated for each lineage assuming the exponential population growth of coalescence model, based on BEAST 2.6.0.: 1.30 × 10^-8 for Lineage N of L. exotica, 1.07 × 10^-8 for Lineage S of L. exotica, and 1.44 × 10^-8 for L. cinerascens. In addition, 10 million steps were run with sampling of every 1,000 generations in the MCMC method. The investigation of ESS value and BSP was conducted using the TRACER 1.6 program (Heled and Drummond, 2008; Rambaut and Drummond, 2014a). Estimation of divergence time was performed on the Ligia phylogeny inferred from 16S rRNA using the BEAST 2.6.0. program. The time calibration was based on fossil records, as estimated by Neraudeau (2008) (Neraudeau, 2008; Broly et al., 2013). A fossil of the genus Ligia was found in France and the age was estimated to be in the Albian stage of the Cretaceous (113.0-100.5 Ma). Based on the fossil record, we set a calibration point (110.0 Ma) where the genus Ligia and other genera had been diverged. After time calibration, divergence time estimation was carried out using a Yule speciation model with the strict clock method. For phylogeny reconstruction, four outgroup sequences including those of Idotea baltica, Ligia oceanica, L. perkinsi, and L. hawaiensis were selected and retrieved from NCBI GenBank. The GTR + F+ G4 model was selected as a best-fit substitution model, estimated from jmodeltest 2.02 (Darriba et al., 2012). A total of 10,000,000 MCMC generations were performed as a posterior distribution parameter. The first 25% of the phylogenies were removed as a burn-in, and the resultant 8,001 trees were combined into a maximum clade credibility tree using the TreeAnnotator 2.6.0 program (Rambaut and Drummond, 2014b). The consensus tree was visualized in FigTree 1.4.2 (Rambaut, 2012). To analyze the hypothetical distribution of a common ancestor, an S-DIVA analysis based on Bayesian Binary MCMC (BBM) methods was conducted using RASP 3.2 (Nylander et al., 2008; Yu et al., 2015), which was based on an ultrametric tree generated using BEAST 2.6.0. Seven localities were assumed and coded as follows: A) South Korea; B) Japan; C) China; D) Taiwan; E) Southeast Asia; F) Africa; G) Americas.

Here, 282 Ligia individuals were collected from ten populations from South Korea (DD, UL, YD, TA, PH, TY, YS, WD, JJ, and SG), three from Japan (SK, IW, and MUD), and three from Vietnam (CB, DN, and PQ) ( Figure 1 , Supplementary Table 1 ). Partial mitochondrial 16S rRNA gene with a length of 379 bp, were sequenced allowing the identification of 198 individuals as L. exotica and 84 individuals as L. cinerascens. Twelve populations were identified for L. exotica (six in South Korea, three in Japan, and three in Vietnam) and six for L. cinerascens (five in South Korea and one in Japan). In only two populations, including JJ of South Korea and MUD of Japan, samples of both species were collected. In total, 41 haplotypes were detected including 23 haplotypes of L. exotica and 18 of L. cinerascens ( Supplementary Table 4 ). No haplotypes were shared interspecifically. From a sequence alignment of the 23 16S rRNA haplotypes of L. exotica, 49 polymorphic sites were identified, which consisted of seven singleton variable sites, 42 parsimoniously informative sites, and one gap ( Supplementary Table 5 ). In the case of 18 16S rRNA haplotypes in L. cinerascens, 19 polymorphic sites were detected including six singleton variable sites, 13 parsimoniously informative sites, and one gap ( Supplementary Table 6 ).

In Table 1 , genetic diversity was examined based on the 16S rRNA sequences for each L. exotica and L. cinerascens population. Among L. exotica populations, the highest diversity was observed in IW, Japan (h = 0.857), and the lowest in YD and SG (h = 0.000), South Korea. When we compared the genetic diversities of the populations between countries, the highest haplotype diversity was observed in Japan (h = 0.896), and the lowest in South Korea (h = 0.251). Nucleotide diversity was the highest in SK (π = 0.019), Japan, and the lowest in YD and SG (π = 0.000), South Korea. Among the examined countries, the highest nucleotide diversity was observed in Vietnam (h = 0.059), and the lowest was observed in South Korea (h = 0.001). In the L. cinerascens populations, the highest haplotype diversity was observed in PH (h = 0.711), South Korea, and the lowest was observed in JJ, South Korea (h = 0.000; only two individuals). Nucleotide diversity was the highest in TA (π = 0.011), South Korea, and the lowest in JJ (π = 0.000; only two individuals), South Korea.

To investigate genetic differentiation among the populations, F _ST values between the L. exotica and L. cinerascens populations were calculated ( Table 2 ). Most of the F _ST values between L. exotica the populations were statistically significant excluding the following eight cases: UL vs. YS, YD vs. YS, YD vs. JJ, YD vs. SG, YS vs. JJ, YS vs. SG, YS vs. IW, and WD vs. MUD ( Table 2A ). The F _ST values ranged from 0.149 (UL vs. YD) to 0.997 (SG vs. PQ) ( Table 2A ). Among South Korea, Japan, and Vietnam, all F _ST values were significant, and they ranged from 0.547 (South Korea vs. Japan) to 0.980 (South Korea vs. Vietnam) ( Table 2B ). In L. cinerascens, all F _ST values were also significant, with ranges from 0.095 (MUD vs. TY) to 0.484(PH vs. DD) ( Table 2C ). The results of the Mantel tests based on the 16S rRNA sequences of the two Ligia species revealed a significant correlation between pairwise genetic and geographical distances (r = 0.548, p < 0.01 for L. exotica; r = 0.393, p < 0.01 for L. cinerascens; Supplementary Figure 1 ).

In addition to the 41 haplotypes of 16S rRNA, 197 16S rRNA sequences of the two Ligia species were retrieved from the NCBI GenBank ( Supplementary Table 3 ). From the retrieved sequences, 92 additional haplotypes were obtained. Based on the resultant 133 haplotypes, ML and BI trees were reconstructed and showed one genetic lineage of L. cinerascens and two genetic lineages of L. exotica with high nodal supports were observed ( Figure 2A ). The phylogeny based on Nak also showed one L. cinerascens and two L. exotica genetic lineages ( Supplementary Figure 2 ). Moreover, three 16S rRNA haplotypes, LCH80, LCH81, and LCH133 exhibited unexpected placements apart from the three genetic lineages of the two Ligia species. Similar patterns were observed in the different analyses including phylogenetic network, TCS network and PCoA ( Figures 2B-D in order) analysis. Especially in the TCS network ( Figure 2C ), the mutation steps between L. exotica and L. cinerascens were more than 40, indicating their genetically apparent differentiation. The mutation steps between the two genetic lineages within L. exotica were more than 30, and there were star-like topologies with the haplotypes such as LCH01, LCH03, LCH23, LCH05, LCH11, and LCH35 at the center, implying recent rapid expansion of the L. exotica population size. According to phylogenetic and population genetic analyses ( Figure 2 ), six distinct genetic lineages in the two Ligia species were observed, including the three major lineages and independent haplotypes of LCH80, LCH81, and LCH133. Their geographic distribution patterns are as follows ( Supplementary Table 7 ). L. cinerascens were mainly distributed in South Korea, Japan, and northern China along the coast of Bohai Sea and Northern Yellow Sea. In case of L. exotica, the two genetic lineages had different geographical distribution patterns Lineage N (North) for one and Lineage S (South) for the other. Lineage N of L. exotica was distributed in South Korea, Japan, and China (along the coastline of Yellow Sea). Conversely, Lineage S of L. exotica was distributed in Japan, China (along the coastline of East China Sea), Taiwan, and Southeast Asia including Vietnam, Cambodia, and India, Africa including Mozambique, and Republic of South Africa, and America including USA, Mexico, Republic of Trinidad and Tobago, and Brazil. Lineage S was not found in South Korea. In the case of LCH80, LCH81, and LCH133, all haplotypes were collected in Japan. Detailed information on the collection sites are as follows: Kanazawa (LCH80) and Kitadaito (LCH81) by Hurtado et al. (2018), and Okinawa (LCH133) (unpublished data).

To investigate whether each of the six genetic lineages in the two Ligia species could be considered independent, four different MSDAs, including ABGD, ASAP, PTP, and bPTP ( Figure 3 , Supplementary Figure 3 ) were conducted based on 16S rRNA. The result of all the analyses strongly indicated that possibility two Ligia species could be divided into six independent species namely, L. cinerascens, Lineages N and S of L. exotica, LCH80, LCH81, and LCH133 supporting the results of the phylogenetic and population genetic analyses ( Figure 2 ). Based on the results of the species delimitation analyses, LCH80, LCH81, and LCH133 were renamed as Ligia sp. 1, Ligia sp. 2 and Ligia sp. 3, respectively.

When genetic distances of 16S rRNA within each lineage of L. cinerascens, Lineage N of L. exotica and Lineage S of L. exotica were calculated using the Kimura-two-parameter (K2P) model, average genetic distance was the highest in the Lineage S of L. exotica (Mean: 0.018), and the lowest in the Lineage N of L. exotica (Mean: 0.008) ( Table 3 ). When comparing distance values across the three lineages, average genetic distance was the highest between Ligia sp. 2 and L. cinerascens (Mean: 0.142), and the lowest between Lineage N of L. exotica and Ligia sp. 3 (Mean: 0.062). The maximum genetic distance value within each lineage did not exceed the minimum genetic distance value among the lineages.

To investigate the fluctuation history of population size within each lineage, three neutrality tests including Tajima’s D, Fu’s Fs, and R ₂ test, were conducted for each lineage ( Table 4 ). As a result, all values obtained from the three neutrality tests were statistically significant in each lineage, suggesting that the population size of each lineage had undergone rapid expansion. When mismatch distribution analysis (MMD) was conducted, L. cinerascens and Lineage N of L. exotica exhibited a bi-modal distribution, respectively ( Figure 4A ), whereas Lineage S of L. exotica exhibited a multimodal distribution. When the raggedness index and sum of square deviation were calculated to examine the mismatch distribution pattern of each lineage, only Lineage N of L. exotica showed as a statistically significant value, implying the rapid expansion of its population size. In BSP analysis for each lineage, increases in effective population size were observed in L. cinerascens and Lineage S of L. exotica and not in Lineage N of L. exotica ( Figure 4B ). The increases in effective population sizes were estimated to have taken place between ca. 250-50 Ka and 350-150 Ka, respectively.

To examine divergence times among the six putative species observed from L. cinerascens and L. exotica, a molecular clock analysis was conducted using the BEAST 2.6.0 program. As shown in Figure 4C and Supplementary Figure 4 , the results suggested that a monoclade of L. cinerascens and Ligia sp.1 was first diverged from the remaining four species in the Oligocene (24. 12 Ma [95% HPD: 18.54 - 30.11 Ma]). Thereafter, Ligia sp.1 diverged from L. cinerascens in the early Miocene (18. 47 Ma [95% HPD: 12.78 - 24.08 Ma]). Not long after that event, Lineage S of L. exotica was diverged in about 18.26 Ma (95% HPD: 13.32-23.59 Ma). In the mid-late Miocene (13-10 Ma), Ligia sp. 2 (13.96 Ma [95% HPD: 9.57-18.35 Ma]) and Ligia sp. 3 (10.16 Ma [95% HPD: 6.68-14.06 Ma]) diverged from Lineage N of L. exotica. Among the three lineages of L. cinerascens, and Lineages N and S of L. exotica, the diversification of haplotypes could have first occurred in the Lineage S of L. exotica, in the late Miocene (8.10 Ma [95% HPD: 5.17-11.33 Ma]). The haplotype diversity of the other lineages emerged in the Pliocene, about 4.55 Ma (95% HPD: 3.14-6.33 Ma) in L. cinerascens and 3.11 Ma (95% HPD: 1.90-4.59 Ma) in Lineage N of L. exotica. The results of S-DIVA analysis under a Bayesian binary Markov chain Monte-Carlo (BBM) model indicated that a common ancestor of the six Ligia species originated from Japan.

The present study investigated the phylogenetic relationships, phylogeographical patterns, and speciation events of two Ligia species, L. exotica, and L. cinerascens, based on 16S rRNA at the population level. The results of the F _ST value and Mantel test showed the genetic differentiation by geographical distance among the populations of the two Ligia species, which suggests the existence of genetic structure among the Ligia populations ( Table 2 , Supplementary Figure 1 ). This possibility was directly shown in the unrooted phylogeny, TCS network, and PCoA ( Figure 2 ). Through the analyses, we found six different and independent lineages as follows: L. cinerascens, Lineages N and S of L. exotica, Ligia sp.1, Ligia sp.2, and Ligia sp.3. The six genetic lineages were also supported through the four MSDAs, suggesting each of the six as an independent species.

Among the 12 populations of L. exotica investigated in the present study, eight populations (UL, YD, YS, WD, JJ, and SG) in South Korea, and two populations (MUD and IW) in Japan belonged to a Lineage N of L. exotica. The remaining four populations, including one population (SG) in Japan and three populations (CB, DN, and PQ) in Vietnam, were included in a Lineage S of L. exotica. In Japan, there existed both lineages N and S of L. exotica. However, only Lineage N was found in South Korea and only Lineage S was found in Vietnam. These patterns of genetic lineages may help to explain the higher haplotype and nucleotide diversities of the total population of L. exotica in Japan (h = 0.896, π = 0.046), compared to those in South Korea (h = 0.251, π = 0.001) and Vietnam (h = 0.519, π = 0.005).

Moreover, we conducted phylogenetic analyses using Nak, which is one of the nuclear markers that has been studied for Ligia species. The results of Nak-based phylogenetic analyses also showed three distinct genetic lineages: L. cinerascens, Lineage N of L. exotica, and Lineage S of L. exotica. This means that the existence of three genetic lineages was also strongly supported based on analyses of both a nuclear and mitochondrial marker. However, the genetic marker Nak may have limitations in the study of species delimitation; only three polymorphic sites in Nak were detected between L. exotica and L. cinerascens, and only one was identified between the L. exotica two lineages ( Supplementary Table 8 ). Other nuclear genetic markers, such as internal transcribed regions (ITS), should be evaluated in future studies to better understand species delimitation.

Yin et al. (2013) discussed the detailed distribution ranges of L. exotica and L. cinerascens based on previous studies in the geographic range of East Asia. The present study reconstructed phylogeographic patterns of L. exotica and L. cinerascens with updated information after Yin et al. (2013) ( Figure 5 ). In Figure 5 , we speculate the northern and southern limits of their distribution ranges. The southern limit of distribution for L. cinerascens was observed to be about 33-35°N. In the case of L. exotica, there were apparent challenges in defining the northern limit in South Korea, because there was no or little information in North Korea. Nevertheless, we can propose a plausible northern limit (39−41°N latitude) of the two genetic lineages of L. exotica based on the geographic records of Lineage N of L. exotica. In sum, we found two different patterns of distribution for the two Ligia species in East Asia. First, L. cinerascens and L. exotica showed overlapping distribution in 33−41°N latitude range. Most of the range of overlapping distribution was in Japan (33−41° N) when compared with that in China (35−39° N). Such a pattern has been reported by Yin et al. (2013), and they discussed that the pattern might be due to sea surface temperatures caused by the Kuroshio current and deeper inshore seawaters around Japan (Yin et al., 2013; Hurtado et al., 2018). Regarding the causes of the overlapping distribution of the two Ligia species, human-mediated activities such as shipping and damming have been mainly cited (Yin et al., 2013). Other factors influencing overlapping distribution might be genetic exchanges by rapid population size expansion, oceanic rafting mediated by large algae, or secondary contact in the past (Jung et al., 2008; Yin et al., 2013).

Second, two genetic lineages of L. exotica exhibited general allopatric distribution, which appears to be more evident in China along the coastline than in Japan. Similarly, such allopatric distribution patterns have been reported in many diverse marine taxa such as Collichthys lucidus (Richardson, 1844), Cyclina sinensis (Gmelin, 1791), Eriocheir De Haan, 1835 , Monodonta labio (Linnaeus, 1758), Oratosquilla oratoria (De Haan, 1844), and Tubaca arcuate (De Haan, 1835) (Xu et al., 2009; Ni et al., 2012; Cheng and Sha, 2017; Zhao et al., 2017; Zhao et al., 2021; Hwang et al., 2022; Shih et al., 2022). Such studies suggested that the allopatric distribution observed in the coastline of China is caused by sea level fluctuations in the Pleistocene (Xu et al., 2009; Ni et al., 2014). Allopatric distribution of the two genetic lineages of L. exotica presented here could have been affected by sea-level change, causing the geographical isolation between the two genetic lineages.

In addition, an analysis of population dynamics was conducted to check population size stability using neutrality test, MMD, and BSP in L. cinerascens, as well as Lineages N and S of L. exotica. In each of the three genetic lineages, population size expansion had occurred in the course of their evolutionary history. Population expansion may also be recognized through the TCS network results, which show complex star-like topologies within each genetic lineage. In BSP, we could confirm the period of expansion of effective population size, except in Lineage N of L. exotica. When Marine Isotope Stage (MIS) based on oxygen isotope values reflecting the temperature change and sea level fluctuation was considered, the expansion events of effective population sizes of L. cinerascens and Lineage S of L. exotica were estimated to be MIS6 (ca. 191-130 Ka) and MIS8 (ca. 300-243 Ka), respectively. MIS6 and MIS8 were glacial periods where sea level decreased and land bridges were formed in the Yellow Sea and East China Sea (Waelbroeck et al., 2002; Shackleton et al., 2003; Jiménez-Sánchez et al., 2013). The land bridge formed during the glacial period in the Pleistocene period could have facilitated L. exotica and L. cinerascens colonization events through the extension of the coastline (Jung et al., 2008).

Based on the results of the investigations on population expansion for each genetic lineage of the two Ligia species, time divergence phylogeny was reconstructed, and S-DIVA analysis was conducted. According to the results, a common ancestor of the three originated from Japan. Moreover, the estimated species divergence time was confirmed to be consistent with the sequence of geological events in the Japanese islands. The present study strongly supports the viewpoint of Hurtado et al. (2018) in which geological events in the Japanese archipelago might be some of the factors causing genetic lineage divergence in L. exotica. Divergence times of Ligia sp. 1, L. cinerascens, and Lineage S of L. exotica were estimated from late Oligocene to early Miocene (about 30-15 Ma). During the period, proto-Japanese Islands located along the eastern margins of the Eurasian continent began to separate from that continent (Tojo et al., 2017). Through the event, the East Sea was formed, and the northeast and southwest portions of the Japan Arc were separated independently (Otofuji et al., 1991; Otofuji et al., 1994; Tojo et al., 2017). The event could have caused the isolation of the habitats of Ligia species and seemed to induce the speciation of Ligia sp.1, L. cinerascens, and Lineage S of L. exotica from all the Ligia lineages. After the split of the Japanese main island from the Eurasian continent, more submergence and upheaval occurred between continental plate and oceanic plate around to early Pliocene (about 5 Ma; Ota, 1998; Tojo et al., 2017). During the period, the other three lineages (Ligia sp. 2, sp. 3, and Lineage N of L. exotica) might have diverged, and the diversification of haplotypes in Lineage S of L. exotica might have already progressed. Afterward, the archipelago in the current shape of the main island of Japan was formed, which might have led to dramatic haplotype diversifications of other genetic lineages, such as L. cinerascens, and Lineage N of L. exotica, and yielded their current genetic lineages.

Numerous endemic species in the genus Ligia have been reported in Japan, such as L. boninensis Nunomura, 1979, L. daitoensis Nunomura, 2009, L. hachijoensis Nunomura, 1999, L. miyakensis Nunomura, 1999, L. ryukyuensis Nunomura, 1983, L. shinjiensis Tsuge, 2008, and L. yamanishii Nunomura, 1990 (Nunomura, 1979; Nunomura, 1983; Nunomura, 1990; Nunomura, 1999; Tsuge, 2008; Nunomura, 2009). The range of distribution of these species was limited only to islands and brackish lakes in Japan: Ryukyu Island (L. ryukyuensis), Izu Island (L. hachijoensis, L. miyakensis), Bonin Island (L. boninensis, L. yamanishii), Daito Island (L. daitoensis), and Shinji Lake (L. shinjiensis). These species have different morphological characteristics, such as flagellar segments of the second antenna, protuberance on the first pereopod in males and the apical part of the stylus (Nunomura, 1979; Nunomura, 1983; Nunomura, 1990; Nunomura, 1999; Tsuge, 2008; Nunomura, 2009). However, there is no genetic information yet on them. If genetic information on the species and investigations on morphological characteristics that can distinguish the two lineages in L. exotica is added through further studies, it may facilitate the understanding of species delimitation and discovery of cryptic species of Ligia in East Asia (Choi et al., 2020; Choi et al., 2021; Shin et al., 2021; Kim and Hwang, 2023).