From November 2020 to November 2022, 164 specimens of N. yoldii were collected from nine intertidal sites along the NWP coastline and muscle tissues were preserved in anhydrous ethanol. Detailed information regarding the collection dates and locations is presented in Table 1 , and the distribution of the sampling stations is illustrated in Figure 1B . The sampling region encompassed three marginal seas within the NWP region. Namely, the Yellow Sea (NT), East China Sea (DYS, SS, NSD, LK, PT, XM), and South China Sea (SZ, ZJ) were included. Concurrently, the South China Sea experienced a decline in sea levels by approximately 100-120 meters, transforming it into a semi-enclosed sea, primarily linked to the Pacific Ocean via the Bashi Strait (Ni et al., 2014). The depth of the NWP increases significantly from north to south. The Yellow Sea has an average depth of only 44 m, while the East China Sea reaches around 350 m, and the South China Sea exceeds 1,200 m. The SS location marks the northern distributional limit of N. yoldii within the NWP since the 1960s (Tchang et al., 1963). SS situated at the mouth of the Changjiang, the third-largest river globally, with an average annual discharge ranging from 8-9×10¹¹ m³. Changjiang exerts significant influence on various ecological and environmental parameters within the East China Sea, separating the North Pacific Temperate Biotic Region (characterized by cold temperate fauna) from the Indo-West Pacific Warm-water Biotic Region (characterized by tropical/subtropical fauna) (Ni et al., 2017). The annual sea surface temperature north of the Changjiang Estuary ranges from 14°C to 19°C, while the sea surface temperature south of the Changjiang Estuary varies between 19°C and 25°C. Consequently, this biogeographic boundary along the NWP coastline has historically restricted the northward dispersion of N. yoldii, although recent decades have witnessed its extension beyond traditional confines (Wang et al., 2022). Notably, investigations in January 2016 (Wang et al., 2018) indicate that the NT location approaches its current northern boundary at approximately 32°51′N.

The samples were subjected to DNA extraction using a Tiangen DP304 blood/cell/tissue genomic DNA kit (Tiangen, Beijing, China). Polymerase chain reaction (PCR) amplification was used to obtain COI gene sequences for phylogenetic analysis. The mitochondrial DNA (mtDNA) COI region was amplified and sequenced using the primer pairs LCO1490-GTCAACAAATCATAAAGATATTGG and HCO2198-TAAACTTCAGGGTGACCAAAAAATCA (Folmer et al., 1994). PCR amplification conditions were as follows: an initial denaturation at 94°C for 3 minutes; 35 cycles of denaturation at 94°C for 45 seconds, annealing at 50°C for 45 seconds, and extension at 72°C for 1 minute; and a final extension at 72°C for 10 minutes. Subsequently, the unpurified PCR samples were subjected to electrophoresis on a 1.5% agarose gel at 120V for approximately 30 minutes. A BigDye Terminator v3.1 sequencing kit (ABI, U.S.) was used for PCR. The PCR products underwent purification via the ethanol precipitation method, followed by drying and storage at 4°C away from light. Sequencing was conducted using the Sanger method (double-stranded DNA chain termination technique), followed by capillary electrophoresis to acquire the DNA base sequence. The dried samples were dissolved in 10 µL of formamide, denatured in a PCR apparatus, and sequenced using the 3730XL platform (Thermo Fisher Scientific, U.S.). The initial sequence assembly was conducted using SeqMan software v.7.1.0, and alignment and trimming were performed using BioEdit v.7.0.9 (Hall, 1999). All sequences were deposited in GeBank, with accession number PP436912-PP437075.

The number of haplotypes for each locality was determined using DnaSP v.5.10. Additionally, the number of polymorphic sites, haplotype diversity, and nucleotide diversity were estimated using Arlequin v3.11 (Excoffier et al., 2005). Bayesian analysis was conducted on all N. yoldii specimens, alongside GenBank sequences of Neritina violacea (accession numbers MF688020), utilizing MrBayes v3.2.4 under the F81 model (Huelsenbeck and Ronquist, 2001), four Markov chains were applied for 1,000,000 generations, sampled every 1000 generations, and the first 25% discarded as the burn-in, respectively. A network depicting the relationships among the haplotypes was constructed using TCS (Clement et al., 2000).

A simulated annealing approach implemented in SAMOVA v. 2.0 was used to test varying numbers of groups (K) to identify phylogeographically homogeneous populations that were maximally differentiated from each other (Dupanloup et al., 2002). The analysis was run for 2 to 8 groups, and the grouping associated with the highest F_CT value was selected for further analysis. Analysis of molecular variance (AMOVA) was conducted for COI data using Arlequin v3.11, with 1,000 permutations to explore the partitioning of genetic diversity among populations (Excoffier et al., 1992). Arlequin v3.11 was also employed to compute the pairwise genetic differentiation (F_ST) among populations and groups identified by SAMOVA (Excoffier et al., 2010). To investigate the relationship between genetic and geographical distances, Mantel’s test with 10,000 randomizations for isolation-by-distance (IBD) (Slatkin, 1993) was performed between the linearized F_ST (F_ST/[1-F_ST]) and geographical distances (shortest marine connection measured in km with ArcMap v10.3) using IBD v1.52 (Bohonak, 2002) to ascertain whether this pattern conformed to the expectation that genetic similarity decreases with geographical distance (Novembre et al., 2008). Linearized genetic differentiation (F_ST/[1-F_ST]) values were plotted against geographical distance using Origin software to visualize the relationship between genetic and geographical distances.

Tajima’s D neutrality test and Fu’s F_S neutrality test (Fu and Li, 1993) were conducted using DnaSP 5, with the p-value serving as the criterion for significance. These tests are based on the assumption that mutations occur in each position of the alignment with equal probability consistently over time. Values greater than zero indicate bottleneck effects and/or non-directional selection, whereas values significantly lower than zero indicate population expansion and/or directional selection. Mismatch distribution analysis of the mtDNA sequence of each lineage was performed (Slatkin, 1995). The historical population sizes were estimated using extended Bayesian skyline plots approach employed with the HKY substitution site model using a strict molecular clock with rate 0.0023 (Frey and Vermeij, 2008) and running 1.0 × 10⁶ Markov Chain Monte Carlo (MCMC) simulations in the software BEAST 2.7 (Drummond and Rambaut, 2007; Heled and Drummond, 2008). The convergence of the MCMC runs was assessed using effective sample size (ESS) values for the parameters. The MCMC conditions were adjusted to achieve ESS values greater than 200 for all parameters ( Supplementary Table S1 ). The results were checked with Tracer 1.7 (Rambaut et al., 2018). The population expansion time (t) of N. yoldii was calculated using the formula t = τ/2u (Rogers, 1995), where u = μkg; μ is the nucleotide substitution rate, which is reported to be 0.0023 substitutions/site/myr for the COI gene in genus Nerita (Frey and Vermeij, 2008); k is the length of COI sequence; and g is the generation time of N. yoldii. A generation time of 3 years was set based on studies of the life history of its sister species, N. japonica (Nakano and Nagoshi, 1980). DIY-ABC v2.1.0 was utilized to gain further insights into the timing and divergence history of N. yoldii based on COI gene sequences (Cornuet et al., 2008; Cabrera and Palsbøll, 2017). We tested five plausible divergence scenarios ( Supplementary Figure 1 ) for the four lineages identified by the Bayesian tree. In all scenarios, lineage B and C were the first to diverge. The scenarios modeled the possible hypotheses regarding the origins of lineages A and D, which could either have originated from lineage C or resulted from an admixture between lineage B and C.

The alignment of the 557 bp COI gene fragment was analyzed in 164 individuals, revealing 384 variable sites and 35 unique haplotypes (Hap1–Hap35). The Bayesian tree for N. yoldii based on the COI gene fragment indicated four lineages: A, B, C, and D, and revealed a polytomy among the lineages ( Figure 1A ). The haplotype network diagram ( Figure 2 ) aligned closely with the Bayesian tree, indicating that lineages B, C, and D were interconnected, while lineage A was distinct from lineages B and D. Lineage C exhibited a star-like topology with Hap5 as the central haplotype. Populations NSD and ZJ exhibited the highest haplotype diversity, whereas populations SS and LK exhibited only four haplotypes. Hap5 was the most widespread haplotype being found in all populations. The average haplotype diversity (h) across all populations was 0.742, and nucleotide diversity (π) was 0.0018. The NSD population displayed the highest haplotype diversity (h = 0.874), whereas the SS population exhibited the lowest (h = 0.437). The highest nucleotide diversity was observed in LK (π = 0.0037), while PT and XM had the lowest (π = 0.0004).

According to the SAMOVA, F_CT was maximized (F_CT = 0.437, P= 0.036) when all populations were divided into three groups: northern (NT, DYS, SS, NSD, and LK), central (PT and XM), and southern (SZ and ZJ). AMOVA indicated that 91.71% of the haplotype variation occurred within populations, while 7.13% occurred among groups and 1.16% among populations within groups ( Table 2 ). Pairwise F_ST comparisons revealed significant differentiation, with SZ and ZJ exhibiting significant differentiation from the other populations ( Table 3 ). However, the Mantel test indicated no significant relationship between the linearized F_ST values and geographical distance across the nine populations ( Figure 3 ). The genetic differentiations among groups were all significant, with F_ST values between the northern and central group higher than those between the northern and southern group ( Table 4 ).

The mismatch distribution and neutral tests of the COI gene were used to detect population expansion. Both Tajima’s D and Fu’s F_S in the lineage C were negative and statistically significant, indicating demographic expansion. Only the mismatch distribution of lineage C was a unimodal Poisson distribution, confirming demographic expansion in its history. For the other three lineages, there were no explicit signals of population expansion, as neither neutral tests were not significant or the observed mismatch distribution did not identify unimodal ( Table 5 ; Supplementary Figure S4 ). Reconstruction of population sizes over time indicated a historical population expansion in lineage C, whereas the other three lineages remained relatively stable ( Figure 4 ). Demographic expansion of the lineageC occurred during the late Pleistocene of the Quaternary period with expansion time of 0.038 Mya. The expansion time predate the last glacial maximum (LGM), which occurred approximately 0.020 Mya. Notably, the expansion time occurred during a period characterized by dramatic sea level fluctuations in the late Pleistocene. Approximate Bayesian computation (ABC) analysis supported scenario 1 as the most probable scenario, with posterior probabilities of 0.250 and 0.260 for the direct estimate approach and logistic regression approach, respectively ( Supplementary Figure S2 ). Model testing analysis demonstrated that summary statistics derived from observed data fell within the range of simulated datasets from the posterior predictive distribution based on scenario 1 ( Supplementary Figure S3 ), indicating robust model performance. Scenario 1 posits that ancient populations diverged into lineages B and C at time t3, and lineage A and D was originated from lineage C at time t2 and t1, respectively. Median divergence time were estimated at 73,700 (10,500-324,000), 155,000 (68,600-435,000), and 255,000 (208,000-493,000) generations for t1, t2, and t3, respectively. Assuming a generation time of 3 years for N. yoldii (Nakano and Nagoshi, 1980), these divergence time correspond to approximately 0.221 (0.032-0.927) Mya, 0.465 (0.206-1.305) Mya, and 0.765 (0.206-1.305) Mya for t1, t2, and t3, respectively.

The haplotype diversity and nucleotide diversity of N. yoldii from the NWP region were comparatively lower than those of N. scabricosta (h = 0.91, π = 0.005) and N. funiculata (h = 0.99, π = 0.008) from the temperate and subtropical East Pacific regions, respectively (Frey and Vermeij, 2008). N. yoldii originated approximately 10.3 Mya during the Miocene epoch (Frey and Vermeij, 2008), This relatively shorter evolutionary history compared to other Nerita species may have constrained the accumulation of genetic variation in N. yoldii. Another possible explanation is that following geographic isolation from the Pacific Ocean during the Pleistocene epoch, small founder populations may have emerged within each lineage of N. yoldii. Subsequent genetic drift within these isolated populations could have led to the reduction of intrapopulation genetic diversity (Hu et al., 2015). Rare variants in small populations are predicted to be eliminated by genetic drift. Recent population isolation and fragmentation during the Pleistocene glacial age may have increased the role of genetic drift in COI variation in N. yoldii (Luo et al., 2012). The dispersal pattern of Nerita’s ancestors from tropical to temperate seas on a large geographical scale is uncommon among gastropods and typically results in low speciation and high extinction rates (Williams, 2007). Additionally, the frequent cycles of glacial and interglacial periods during the Quaternary period likely led to unstable population sizes and low genetic diversity in N. yoldii populations. Despite these factors, only 1.16% of the variation was found among the populations of N. yoldii in this study ( Table 2 ), suggesting limited genetic differentiation among the populations of this species.

The distribution of COI haplotypes of N. yoldii exhibited a clear geographical pattern ( Figure 1B ), with significant differences among the northern, central and southern groups. Studies of larval diffusion simulations of N. yoldii along the coast of Jiangsu Province have revealed that larvae dispersed as far north as the Jiaodong Peninsula and as far south as the Yangtze River Estuary (Wang et al., 2020). However, there was little gene exchange between the northern and southern groups likely because of significant distances and temperature differences.

Mismatch analysis and neutral tests suggested that the lineage C of N. yoldii in China underwent expansion event. The expansion of the lineage C is estimated to have occurred approximately 0.038 Mya. The population expansion observed in our study occurred prior to the LGM. Traditionally, the LGM has been considered a pivotal factor driving the evolution and genetic structure of marine organisms (Hewitt, 2004). However, recent studies have highlighted the profound impact of Pleistocene glacial-interglacial cycles, including dramatic sea level fluctuations occurring 30-450 thousand years ago, which predate the LGM (Waelbroeck et al., 2002). These cycles have been shown to have substantial influences on population expansions in several species within the Northwestern Pacific (NWP) region (Ni et al., 2014). Despite the absence of a definitive explanation, substantial evidence of pre-LGM expansion has been documented in various marine realms, including the NWP (Marko et al., 2010), collectively emphasizing the importance of pre-LGM events in shaping the demographic trajectories of marine populations. During the Pleistocene glacial period, sea levels fluctuated due to the melting of glaciers, resulting in significant changes in coastal habitats (Wang, 1999). As a result, some marine species appear to have undergone dramatic demographic size changes during the Pleistocene when the sea level fluctuated (Ni et al., 2014). Intertidal species, such as N. yoldii likely experienced demographic expansion and rapidly occupied new habitats as the sea levels rose.

Our data suggested that the lineage divergence of N. yoldii in the NWP likely resulted from recurrent sea area isolation and connection during the mid to late Pleistocene. According to approximate Bayesian Computation (ABC), lineage A and D originated approximately 0.221 and 0.451 Mya, highlighting the significant influence of sea level changes. Lineage A and D are believed to have originated from lineage C during periods of low sea levels, which could have caused isolation between lineages (Waelbroeck et al., 2002). Therefore, the recurrent separations of these marginal seas during glacial-interglacial oscillations are considered significant drivers that influence the abundance and genetic structure of marine taxa (Hu et al., 2011; Cheang et al., 2012; Ni et al., 2014). The glacial period led to a decrease in sea level, causing the marginal sea of China to shrink to approximately 66% of its present size (Waelbroeck et al., 2002). Even when sea levels were low, there were still refugia in the NWP (Lee et al., 2012), in addition to the Okinawa Trough. The ancestors of lineages A and D were likely separated from lineage C during these periods of sea-level decline. These refugia provided suitable environments but restricted their dispersal, leading to relatively stable historical population sizes for lineages A and D compared to the more widespread lineage C. As a result, neither the mismatch distribution nor the extended Bayesian skyline plots indicated demographic expansion in the histories of lineages A and D. Various geological formations, such as lowland plains along the northwest coast of the Pacific Ocean and mountain ranges, including the Changshan Mountains of Vietnam and the Middle Cordillera Mountains of the Philippines, serve as barriers to coastal dispersal. Moreover, the South China Sea became largely landlocked to the south, west, and north, with its east coast blocked by the Sunda Isthmus (Wang, 1999). The only connection to the Pacific Ocean is through the Bashi Channel between Taiwan and the Philippines. These geographical features likely contributed to the isolation of populations in marginal seas from those in the open ocean, leading to genetic differentiation and historical expansion events following rising sea levels.

Previous studies have suggested that artificial structures can facilitate the invasion and dispersal of species in lithofacies intertidal zones (Dong et al., 2012; Adams et al., 2014). The haplotype diversity observed in the three populations from the northern group (NT, DYS, and SS) was low. This low diversity can be attributed to several factors including the geographical location of these populations and recent environmental changes. The Yangtze River Estuary, where these populations were sampled, has traditionally been regarded as the northern boundary of the distribution of N. yoldii (Wang et al., 2018). The haplotype diversity of these marginal populations was relatively low compared with that of populations in other regions. Climate change and the construction of artificial hard-bottom structures north of the Yangtze River Estuary have contributed to the expansion of the northern boundary of the distribution of N. yoldii (Wang et al., 2020). These changes may have affected the genetic diversity of populations in this region. Lineage B is predominantly distributed in the South China Sea. Analyses including mismatch distribution and extended Bayesian skyline plots indicate that its historical population size remained relatively stable. The South China Sea, being larger in area compared to the East China Sea, and maintaining connectivity with the Pacific Ocean during Pleistocene sea level decreases due to its deep waters (Wang, 1999), experienced lesser impact from Pleistocene sea level fluctuations. Consequently, ecosystems in the South China Sea were less affected by these fluctuations. The unique haplotype (Hap 23) observed in the NT population from the north of the Yangtze River Estuary suggests that this population has undergone in situ differentiation and is genetically distinct from other populations.

The impacts of climate change, including factors such as rainfall, drought, floods, and typhoons, as well as changes in seabed substrate hydrodynamics and water depth, can significantly influence species dispersal at various spatial scales (Pyron and Covich, 2003). Environmental changes can affect population reproduction, spawning, and larval dispersal dynamics (Lord et al., 2012). Specifically, changes in sea surface salinity, temperature, and other hydrological factors caused by diluted water in the Yangtze Estuary may act as potential barriers to larval spread and gene exchange among species along the estuary (Cheang et al., 2008; Ni et al., 2017). However, the findings of this study on the geographical distribution of Nerita indicate that the low-salinity area of the Yangtze River Estuary may not necessarily hinder the larval diffusion of Nerita. Interestingly, the lack of a hard substrate on the north coast of the Yangtze River and the temperature difference between the Yellow Sea and the East China Sea may be important factors limiting the further spread of N. yoldii to the north (Liu, 2013; Wang et al., 2018). As offshore water temperatures continue to rise in China, the population of N. yoldii will continue to extend northward, potentially forming a new distribution pattern. However, it is important to note that this study included only one mitochondrial gene sampled from nine populations. Therefore, the genetic diversity and faunal composition of Nerita populations distributed along the extensive coastline of the country may not be fully captured. Future studies should incorporate a larger number of samples and additional genetic markers, such as microsatellites and SNPs, to provide a more comprehensive understanding of the phylogenetic relationships among the N. yoldii lineages.