Microsatellite marker development and population genetic analysis revealed high connectivity between populations of a periwinkle Littoraria sinensis (Philippi, 1847)*

Mengyu LI , Yuqiang LI , Tengfei XING , Yulong LI , Jinxian LIU ,**

1 CAS Key Laboratory of Marine Ecology and Environmental Sciences, Institute of Oceanology, Chinese Academy of Sciences,Qingdao 266071, China

2 Laboratory for Marine Ecology and Environmental Science, Pilot National Laboratory for Marine Science and Technology(Qingdao), Qingdao 266237, China

3 Center for Ocean Mega-Science, Chinese Academy of Sciences, Qingdao 266071, China

4 University of Chinese Academy of Sciences, Beijing 100049, China

Abstract Periwinkle Littoraria sinensis is a widely distributed gastropod in rocky intertidal zone of the northwestern Pacific. To examine genetic diversity and genetic connectivity among coastal populations of L. sinensis in China, 1 636 pairs of primers were successfully designed using whole-genome shotgun sequencing, de novo assembly, and a bioinformatics pipeline QDD. Twelve highly variable polymorphic markers were selected to genotype 351 individuals from 15 populations. Data of nine microsatellite loci were retained for population genetic analysis, and weak genetic differentiation among populations were detected,suggesting high gene flow among populations. The long planktonic larval duration of L. sinensis might have played an important role in the high gene flow among populations. A tendency of genetic differentiation between north and south populations of L. sinensis was detected, which might be resulted from isolation due to lowered sea level in the last glacial maximum. Furthermore, the newly founded populations along the coast of Jiangsu Province were closely related to populations to the south of the Changjiang (Yangtze) River estuary, suggesting that the main source of the newly founded populations is from natural rocky populations south of the estuary.

Keyword: population genetics; Littoraria sinensis; gastropod; simple-sequence repeats

1 INTRODUCTION

The rocky intertidal zone, one of the most intensively studied ecosystems globally, is subject to environmental challenges posed by both aquatic and aerial climatic regimes (Raffaelli and Paine, 1995). It is an important coastal system, which is among the most artificially and heavily disturbed marine ecosystems (Menge, 2000). Periwinkles are important players in shaping the structure of intertidal communities (Apolinário et al., 1999). They can prey on other mobile predators and compete for space with algae and mussels on the hard-bottom intertidal zone(Cubit, 1984). The littorinid genusLittorariais closely associated with rocky shores and mangroves(Reid et al., 2010).Littorariasinensis(Gastropoda:Littorinidae), a periwinkle species, is one of the most conspicuous and abundant gastropods in the northwestern Pacific and distributed on the coast of Korea, Japan, and China.L.sinensisis a newly described species fromLittorariaarticulatabased on a morphology analysis (Reid, 2001). As an intertidal dominant species,L.sinensiscould play a significant ecological role in the intertidal communities (Reid,1989).

Littorariasinensisinhabits hard bottom habitat including rocky shores, trunks at the seaward edge of mangrove forests and artificial structures in the intertidal zone (Dong et al., 2016). The estimated planktonic larval duration ofL.sinensiswas 3-10 weeks (Reid, 2001). In the distribution range ofL.sinensis, the increasing of artificial structures along the coasts of Jiangsu Province, China, provided new suitable habitats forL.sinensisand other intertidal marine invertebrates; and changed the intertidal communities (Dong et al., 2016). Before the construction of the artificial structures, the large mudf lat of this area was probably an effective barrier to larval dispersal. For most benthic marine species,dispersal and exchange ofindividuals among populations occur primarily during the pelagic larval stage (Cowen and Sponaugle, 2009). Defining dispersal patterns of marine organisms is challenging since it could be affected by many factors including survival and development rates of larvae (Doropoulos et al., 2016), hydrodynamic regimes, and seascapes(Cowen et al., 2007; Cowen and Sponaugle, 2009).These challenges could be overcome by using molecular markers to evaluate the population dynamics of marine benthic species (Selkoe et al.,2016). Phylogenetic studies includingL.sinensishave been conducted using molecular sequences,such as nuclear (28S rRNA) and mitochondrial (12S rRNA, COI) (Reid et al., 2010). In addition, a study using mitochondrial COI sequences demonstrated significant genetic differentiation between north and south populations ofL.sinensis(Dong et al., 2016).The newly formed populations on the artificial structures south of Wanggangzha (33°10′N, 120°50′E)were significantly different from north natural populations and shared high genetic similarity with populations south of the Changjiang (Yangtze) River estuary (Dong et al., 2016). In addition, the artificial structures have promoted the northward distribution range shift ofL.sinensisacross the Changjiang River estuary. However, the power of the molecular markers adopted in Dong et al. (2016) for tracing the source ofindividuals in the newly formed populations is limited due to the low intraspecific variability and maternal inheritance of mitochondrial COI.

Higher-resolution molecular markers, such as microsatellites / simple sequence repeat (SSRs), may help to better understand the population patterns and draw convincing conclusions of phylogeographical structure (Selkoe and Toonen, 2006). Microsatellites are widely distributed throughout the genome, and are co-dominant and highly polymorphic (Estoup et al.,1993). It is one of the most widely used markers in molecular genetic applications, including population genetic analysis, germplasm conservation, forensic genetics, parentage analysis and the study of genealogy (Sharma et al., 2007; Schoebel et al., 2013;Flanagan and Jones, 2019). Microsatellite development was promoted by the high throughput and low-cost Next-Generation Sequencing (NGS)technologies, especially in non-model organisms with limited genomic information (Schoebel et al., 2013;Lopez-Marquez et al., 2018; Rosa et al., 2019).Microsatellites have been developed for several Littorinidae species (Tie et al., 2000; Kennington et al., 2008; Carvalho et al., 2015). However, most of them adopted the traditional methods which are laborious, time-consuming, and expensive. For the genusLittoraria, no microsatellite markers have been developed and the genomic resources of this genus are also lacking.

We investigated the population genetic structure ofL.sinensisusing newly developed microsatellite markers developed by the whole-genome shotgun sequencing approach, and conducted population genetic analyses of 15 geographic populations ofL.sinensisalong the coast of China using 12 polymorphic microsatellites, to clarify the population genetic structure across the distribution ofL.sinensis,to check the source of the newly founded populations along the coast of Jiangsu Province, to improve our understanding of the larval dispersal process of marine intertidal organisms, and to provide references to the ecological connectivity of this region.

2 MATERIAL AND METHOD

2.1 Sampling and DNA extraction

A total of 351 individuals ofL.sinensiswere collected from 15 geographic populations from natural rocky intertidal habitats and artificial structures along the coast of China (Table 1). Foot-muscle tissue was dissected and preserved in 95% ethanol. Genomic DNA was isolated using standard phenol-chloroform extraction protocol and DNA quality was checked using 1% agarose electrophoresis and Nanodrop 2000c spectrophotometer.

2.2 Whole-genome sequencing and de novo assembly

Foot-muscle tissue of one individual from Zhangjiatai (Rizhao, Shandong Province) was selected for whole-genome shotgun sequencing. A high-quality DNA sample was cut into fragments of about 500 bp using the Covaris S220 (Covaris,Woburn, MA, USA). A genomic library with insert size ～350 bp was constructed following the Illumina DNA sample preparation protocol using NEB Next Ultra DNA Library Prep Kit (NEB, USA). Finally,sequencing was performed to generate 150-bp pairedend reads on the HiSeq X Ten platform (Illumina) at Novogene Bioinformatics Technology (Tianjin,China).

Table 1 Sampling information of L. sinensis

The raw reads generated by Illumina sequencing were preprocessed to remove adapter fragments,ambiguous reads (i.e., reads with more than 10%unidentified nucleotides) and low-quality sequences(reads with ＞50% bases having phred quality ＜20).Then, de novo assembly was carried out using Celera Assembler 8.3 with default settings (Denisov et al.,2008). The scaffolds were used for further microsatellite discovery.

2.3 Microsatellite identification and primer design

Microsatellite discovery and primer design were performed with QDD 3.1.2 (Meglécz et al., 2010) in a local galaxy platform. Sequences containing microsatellites which were defined as pure tandem repeats of di- to hexa-nucleotide motif with at least five uninterrupted repeats, and with 200 bases of f lanking regions on both sides were detected and used for primer design in Primer3 (Rozen and Skaletsky,2000). Primer parameters were set as follows: length from 18 bp to 27 bp with 20 bp as the optimum length;melting temperature between 57 °C and 63 °C with 60 °C as the optimum, and PCR product size between 90 bp and 320 bp. In order to improve the success rate, only primers and loci that passed the following criteria were retained: one primer pair for each locus;only in design category A (target regions without multiple microsatellites, nanosatellite, and homopolymers); and alignment score below 10.

2.4 Validation of microsatellites

Polymorphism analyses were performed on a total of 60 pairs of primers, which were randomly selected out of the 1 636 retained pairs of primers. We used four individuals from QD (Taiping Bay, Qingdao,Shandong Province) and two individuals from ZJT(Zhangjiatai, Rizhao, Shandong Province) to test whether the primers could be PCR amplified with regular primers or not. The successfully amplified primers were then selected for M13-tailed primer method PCR. The cost-effective M13-tailed primer genotyping protocol (Schuelke, 2000) was then adopted to check the polymorphism of microsatellites.The primer mix contains a forward primer with an M13-tail (5′-GGAAACAGCTATGACCATG-3′) on the 5′ end, a universal fluorescent-labelled M13-tail and a reverse primer. A 10-μL PCR amplification volume was adopted, which contains 10-ng genomic DNA, 5-μL 2× PCRmix (Dongsheng Biotech Co.,China), 0.025-μmol/L M13-forward primer,0.25-μmol/L reverse primer, 0.25-μmol/L M13 primer labeled with different fluorescent dyes (FAM, HEX,TAMRA, ROX), and 3.25-μL ultrapure water.Amplifications were performed with the following cycling conditions: initial activation step at 94 °C for 4 min, 30 cycles consisting of denaturation for 30 s at 94 °C, 40 s of gradient annealing at 50-60 °C, and 40 s of extension at 72 °C, followed with 8 cycles consisting of denaturation for 30 s at 94 °C, 40 s of annealing at 53 °C, and 40 s of extension at 72 °C, and a final extension at 72 °C for 10 min. The annealing temperature of 53 ℃ for the eight cycles was based on that of the fluorescent dye-labeled M13(-21)primer (Schuelke, 2000). The PCR products were electrophoresed on a 1.5% agarose gel and only microsatellite loci that produced specific products were sent to Tsingke Biological Technology (Qingdao,China) for capillary electrophoresis using 3730xl DNA Analyzer with the GeneScan™ 500 LIZ™ dye Size Standard.

2.5 Microsatellite genotyping, genetic diversity,and population genetic structure analyses

After the test of validation, 12 microsatellite loci were retained for population genetic analyses.GeneMarker 2.2 (SoftGenetics, State College, USA)was used for allele calling and manual corrections were performed to minimize genotyping errors.Genetic diversity indices including the number of alleles per locus (Na), observed heterozygosity (HO),expected heterozygosity (HE), polymorphism information content (PIC), allelic richness (AR), and inbreeding coefficient (FIS) were calculated using the Excel Microsatellite Toolkit (Park, 2002) and FSTAT(version 2.9.4). Deviations from Hardy-Weinberg equilibrium (HWE) and genotypic linkage disequilibrium (LD) were determined using Genepop 4.5.1 (Rousset, 2008). The significance levels of tests were corrected with standard Bonferroni correction.

Population genetic differentiation was assessed between pairs of populations by calculating pairwiseFSTandDestvalues (two types of population differentiation statistic) with ARLEQUIN 3.5.2 and GenAlEx (Excoffier and Lischer, 2010; Peakall and Smouse, 2012). GlobalFSTwas also calculated with Genepop 4.5.1. Then, a multidimensional scaling analysis (MDS) was performed with SPSS 23 (SPSS Inc., Chicago, IL, USA) to visualize patterns of population structuring. Spatial analysis of molecular variance (SAMOVA) was performed in SAMOVA 2.0 (Dupanloup et al., 2002) to define groups of populations that are geographically homogeneous and maximally differentiated from each other, which was run with two to five groups and 1 000 permutations.Bayesian clustering of microsatellite genotypes was performed in STRUCTURE 2.3.4 (Pritchard et al.,2000). Ten replicate runs were performed with the values ofKranging from 2 to 10 (instead of pursuing genetic structure between all of the 15 populations,Kis only set to 10 to save time), 1 000 000 Markov Chain Monte Carlo iterations for each run, and a burn-in period of 100 000 iterations with prior knowledge of sampling locations. StructureSelector(Li and Liu, 2018) was used to determine the number ofKthat best fit the data by calculating Ln Pr(X|K)(Pritchard et al., 2000) and ΔK(Evanno et al., 2005).The software POPTREE 2 (Takezaki et al., 2010) was also used to construct a UPGMA phylogenetic tree based on theDSWdistance (Shriver et al., 1995) which takes into account the mutational pattern of microsatellite loci. The number of bootstrap replicates was set at 1 000. Isolation-by-distance (IBD) analyses were conducted for the seven natural rocky populations using IBDWS (Jensen et al., 2005).Google Earth version 4.3 was used to estimate the shoreline geographic distances between sampled populations. Correlation between geographic distance and pairwiseFST/(1-FST) were tested with 10 000 randomizations of the data.

We investigated whether the newly founded populations along Jiangsu Province were affected by genetic drift resulting from bottleneck and/or founder effect. The occurrence of genetic bottleneck was tested with the program Bottleneck 1.2.02 (Piry et al.,1999). Both the Wilcoxon sign-rank test for heterozygosity excess and the sign test that compares the number of loci that present a heterozygosity excess with the number of such loci expected by chance alone was conducted under three mutational models: the infinite alleles model (IAM), the stepwise mutation model (SMM), and the two-phase model(TPM). In the TPM, the default settings of 30%variation from the IAM model and 70% from the SMM model were adopted. Furthermore, effective population size for each population was estimated with NeEstimator v2 by using the linkage disequilibrium method (Do et al., 2014).

3 RESULT

3.1 De novo assembly and microsatellite isolation

In total, 50 579 079 paired-end reads (15.17-Gb data) were obtained. After quality filtering and de novo assembly, 254 132 scaffolds and 259 738 contigswere successfully assembled. The mean bases in scaffolds were 2 045 bp with an N50 of 2 137 bp, and the GC content was 42.29%. Scaffolds were then used for microsatellites detection. A total of 132 044 sequences possessing microsatellite motifs were identified and 103 971 pairs of primers were successfully designed. After filtering, 1 636 primer pairs were retained (Supplementary File S1). The most common repeat motifs of microsatellites were dinucleotide (67.91%), followed by trinucleotide(22.68%), tetranucleotide (8.31%), pentanucleotide(0.86%), and hexanucleotide (0.24%).

Table 2 Characterization of 16 polymorphic microsatellite loci in 6 L. sinensis individuals

Among the 60 primer pairs used for primer testing,33 were successfully amplified by PCR(Supplementary File S2), of which 16 were successfully amplified microsatellites, showing specific and polymorphic amplification products using labelled M13-tailed primer approach (Table 2).According to the preliminary validation results, we randomly selected 12 loci for subsequent polymorphism survey in the 351 individuals from 15 populations.

3.2 Genetic diversity

The 12 microsatellite loci displayed a high level of polymorphisms. Significant deviation from HWE was detected in three loci in all populations (ZH30, ZH34,ZH53), and no LD was detected. These three loci were removed and we retained 9 loci for downstream population genetic analyses (Supplementary File S3).Among populations, the average number of alleles per locus per population (Na) ranged from 9.40 to 12.70. The average expected (HE) and observed (HO)heterozygosity across loci ranged from 0.75 to 0.81 and from 0.44 to 0.64, respectively. The average PIC ranged from 0.71 to 0.77.FISin all 15 populations were positive and ranged from 0.21 to 0.43 (Table 3).

Table 3 Summary statistics for 15 populations of L. sinensis screened with 9 polymorphic microsatellites

3.3 Population genetic differentiation and structure

For comparison purposes, the 15 populations were grouped as north group (QD, ZJT), central group(LYG, ZDZ, SYG, SYHK, ZGD, HG, DFG, and JD populations located along Jiangsu Province, north of the Changjiang River estuary) and south group (YS,YQ, QZ, ZZ, and ST populations located south of the Changjiang River estuary) (see Table 1). The globalFSTwas 0.01 and pairwiseFSTranged from -0.005 to 0.035 andDestvalues ranged from -0.050 to 0.099.Populations that are more geographically distant tended to show higher and significantFSTandDestvalues, especially those between natural intertidal populations from north and south of the Changjiang River estuary (Table 4). TheFSTbetween the north group and central group, north group and south group,central group and south group are 0.014 3, 0.009 8,and 0.001 2, respectively.FSTandDestbetween thenorth group and populations from south of the Changjiang River estuary (south group) were generally low but significant. The pairwiseFSTandDestwithin north group was not significant (P＞0.05).TheFSTvalues between five populations from the central group (SYHK, DFG, HG, ZGD, and JD) and the south group were less than 0.02 and not statistically significant, which were much lower than those between population from north group and south group. Meanwhile, theDestvalues between the central group and the south group also showed a significant decrease. The MDS plots showed a trend of population differentiation between the north group and all the other populations (Fig.1a).

Table 4 Degree of genetic differentiation F ST values (below diagonal) and D est values (above diagonal) among 15 populations of L. sinensis

Fig.1 Analysis of population structure across 15 geographic populations

The result of STRUCTURE analysis (Fig.1c) was similar to the results of MDS. The most likely value ofK=3 was supported by a peak in values of ΔK(Fig.1b). But the values of Ln Pr(X|K) showed two plateaus and the peak values atK=10 (Supplementary File S4), which indicated the high genetic polymorphism between different populations.Furthermore, genetic divergence between the newly founded populations in the northern (LYG, ZDZ,SYG, and SYHK) and southern (ZGD, HG, DFG, and JD) parts along Subei Shoal (Central group) was detected. As for the UPGMA tree based onDSWdistance, the bootstrap values were extremely low,except the clustering of QD and ZJT (bootstrap: 67)(Fig.2), which could only indicate that there was a trend of differentiation between the north group (QD,ZJT) and other populations. The genetic and geographic distance matrices were significantly associated in the IBD analysis (R=0.79,P=0.000 6)(Fig.3).

Fig.2 UPGMA phylogenetic tree based on the D SW distance among 15 populations

3.4 Population demography

The results of the two tests of bottleneck effect were consistent (Table 5). In the Sign test, all of the populations were found to be at mutation-drift equilibrium under IAM and TPM. As for the analysis under SMM, the bottleneck effect was found in ninepopulations. In the Wilcoxon test, significant bottleneck in SYG and ZGD were found under IAM,while a significant bottleneck effect found was found in nine populations under SMM. The estimates of effective population size were infinite for nine populations, five of which were natural populations(Table 3). The six remaining populations showed finiteNe, and four of them were along the coast of Jiangsu Province. However, the upper bound of the confidence interval in all of the populations was inf inity. Small sample size or sampling error can lead to infinite estimates ofNe(Do et al., 2014).

Table 5 Results of the BOTTLENECK tests in 15 populations of L. sinensis

4 DISCUSSION

4.1 Marker development

The NGS approach has been widely used for developing microsatellite markers in ecological and evolutionary important non-model species including plants, mammals, and mollusks (Göl et al., 2017; Liu et al., 2017; Lopez-Marquez et al., 2018). In the present study, we developed abundant candidate microsatellites for the periwinkleL.sinensisin a costeffective manner by using an NGS-based method.Compared with the traditional microsatellite development approaches, the NGS-based method could enhance the marker development efficiency in marine gastropods and promote evolutionary and ecological studies in mollusks. The microsatellite markers developed here, to our knowledge, were the first reported microsatellite resources in this ecologically and evolutionarily important family,which could facilitate the applications of microsatellites in other species of the family by crossspecies amplification (Pérez et al., 2008; Schoebel et al., 2013).

Fig.3 Plot of pairwise estimates of geographic distance and F ST/(1- F ST) between populations of L. sinensis

ForL.sinensis, the present study provided a plentiful resource of microsatellites. A total of 1 636 primer pairs were designed and retained after quality control. We successfully amplified 33 of 60 (55.00%)microsatellite loci and 16 (26.67%) worked for the M13-tailed primer method, which could reduce the cost of fluorescent primer labelling (Hayden et al.,2008). However, the amplification success rate was not very high compared with studies that used similar approaches (Xue et al., 2017; Zhang et al., 2019). For species with large population sizes including mollusks, genomic features including unstable f lanking sequences and occurrences of cryptic repetitive DNA could be responsible for PCR interference or inconsistencies (McInerney et al.,2011). Although challenges in assembling the draft genome ofL.sinensisincluded high levels of polymorphism and abundant repetitive sequences as in many marine invertebrates (Zhang et al., 2012,2019), our results still reinforced that NGS is an important and efficient tool for developing genomewide microsatellites.

The tested loci have a relatively high value for PIC(0.65). Even though significant deviation from the HWE was observed in three out of the 12 selected loci caused by heterozygote deficiency, and the expected heterozygosity was higher than the observed heterozygosity.FISwas positive and very high, which indicated a heterozygosity deficiency that may result from the presence of null alleles and/or inbreeding.Obvious heterozygosity deficiency might result from the presence of null alleles, inbreeding or the Wahlund effect (Andrade et al., 2005; Plutchak et al., 2006; Hui et al., 2017). Deviations from HWE and the presence of a large number of null alleles were also observed in other marine gastropods includingLittorinasaxatilis(Panova et al., 2008),Littorinascutulata(Holborn et al., 2015),Gibbuladivaricate(Lopez-Marquez et al.,2016) and many other marine mollusks including the pacific oysterCrassostreagigas(Hedgecock et al.,2004) and the Wedge ClamDonaxtrunculus(Rico et al., 2017), which seems to be common in marine invertebrates (Panova et al., 2008). As demonstrated in other studies, deviations from HWE can be an indicator of various biological processes such as selffertilization, allele scoring bias, aneuploidy, parental imprinting, spatial and temporal Wahlund effects,age-effects, and direct natural selection on the marker loci (Plutchak et al., 2006). Null alleles at microsatellite loci could have limitations on the estimate of population differentiation (Rico et al., 2017). Null alleles are found in a wide range of taxa (Dakin and Avise, 2004), but it is particularly common in populations with high effective population sizes(Chapuis and Estoup, 2007), including mollusks.Simulation studies and empirical data sets showed that null alleles with frequencies between 5% and 8%should have only minor effects on classical estimates of population differentiation, although higher null allele frequencies may inflateFSTand genetic distances (Chapuis and Estoup, 2007). Microsatellites are widely and routinely used as molecular markers in diverse genetic studies, the effectiveness of microsatellites, particularly forL.sinensis, needs further validation and the results should be interpreted cautiously.

4.2 Population genetic differentiation and structure

Both the results ofFSTandDestand MDS plots supported genetic differentiation between populations in the north group and the other populations, which was consistent with results of the previous study based on mitochondrial COI sequences (Dong et al.,2016), suggesting significant differentiation between the north group and populations south of Lianyungang,Jiangsu Province. Similar patterns have been observed in marine fishes, mollusks, and crustaceans (Ni et al.,2014), where the phylogeographic structure splits could have been induced by historical glaciation in the Northwest Pacific. During the Last Glacial Maximum,L.sinensispopulations might have been isolated in two different marginal seas when most of the East China Sea was exposed (Wang, 1999), and hence, two distinct lineages diverged in that time.However, considering the result of STRUCTURE analysis and the low bootstrap values in the phylogenetic tree, the results may be inconclusive and therefore the conclusion should be taken with caution.

For COI data, the dominant haplotypes were shared by both the north and south populations (Dong et al.,2016). In the present study, the population genetic structure revealed byFST,Dest, STRUCTURE analysis and phylogenetic clustering of populations was weak,which implied high gene flow among populations ofL.sinensis. The estimated planktonic larval duration ofL.sinensiswas 3-10 weeks and the pelagic larvae could be transported by ocean currents to considerable distance during the long pelagic period (Becker et al.,2007). Additionally, many marine invertebrates are able to raft on buoyant macroalgae, which could enhance dispersal over long distances (Nikula et al.,2011). The long planktonic dispersal phases may promote the high population connectivity and maintain weak genetic structure among populations ofL.sinensis.

4.3 New populations on the artificial structures

Both the long planktonic larval stage and new habitat could play important role in structuring populations ofL.sinensis. Along the coast of Jiangsu Province, recently constructed artificial structures have provided suitable habitat for intertidal species andL.sinensiswas the pioneer species (Dong et al.,2016). SAMOVA analysis showed the maximized genetic differences among groups (FCT=0.034 87;P=0.000) was found when populations were divided into three groups (group 1: QD, ZJT; group 2: LYG,ZDZ, SYG, SYHK, ZGD, HG, DFG, JD, YS, YQ,QZ, ZZ; group 3: ST). The results of genetic distance and the clustering analyses indicated that these newly founded populations (central group) were closely related to the natural populations south of the Changjiang River estuary. Consistent with the conclusion in Dong et al. (2016), these newly founded populations demonstrated the northward spread of populations south of the Changjiang River estuary,which could be accomplished by larval dispersal. Our result suggests that the Changjiang River diluted water did not serve as a barrier for the northward shift ofL.sinensis, which was different from previous research suggested (Wang et al., 2015). The coastal current in the Yellow Sea and the East China Sea could have played an important role in the dispersal of planktonic larvae. The energy of flow regime from south to the north could transport the planktonic larval ofL.sinensisand enhance the population connectivity across the Changjiang River estuary. Moreover, the small effective population sizes, the bottleneck effect,and the low genetic diversity in some populations on the artificial structures were possibly the results of the founder effect resulting from the colonization of the new habitats. In addition, the significant bottleneck in QD and ZJT under the SMM model might represent characteristics of peripheral populations (Vucetich and Waite, 2003; Göstring et al., 2010).

Little is known of reproductive patterns ofL.sinensis, but according to the research ofLittorariapallescensandLittorariamelanostoma(Ng and Williams, 2012), we speculated that this species spawn frequently or rapidly returns to reproductive condition after spawning. Therefore, the overlap between generations could be common. The population sizes in our study were small and the linkage disequilibrium method that we used assumes discrete nonoverlapping generations (Do et al., 2014),therefore, the precision of theNeestimates may be reduced by overlapping generations, sample size(Antao et al., 2011), and age structure (Robinson and Moyer, 2013). For wildlife populations, it is much difficult to accurately estimate the effective population size owing to the lack of relevant information (Wright et al., 2021). Although the accuracy of population size estimation is uncertain,Neis recommended primarily as a metric to detect changes over time (Pierson et al.,2018). The strong genetic drift effect in the newly colonized populations with short history could underestimateNe(Wang and Whitlock, 2003).

Newly-formed populations could provide new sources of larvae, create new dispersal pathways and enhance gene flow (Bulleri and Airoldi, 2005; Glasby et al., 2007; Adams et al., 2014). Results of the serpulidPomatocerostriqueterand the limpetPatellacaeruleahave shown that newly formed hard-bottom populations on artificial structures can affect the genetic composition of natural populations (Fauvelot et al., 2009, 2012). Our results show that some populations in the central group were not significantly differentiated from natural populations in both north and south of the Changjiang River estuary. It is thus necessary to investigate the exact contribution of north and south natural populations to the formations of these newly formed populations using highresolution genetic markers, such as genome-wide genetic markers (Sandoval-Castillo et al., 2018;Cornwell, 2020). However, marker selection should be made with more caution, as putative selected SNP marker would not be efficient for detecting the phylogeographic structure and the founder effects.

5 CONCLUSION

In conclusion, the present study provided abundant microsatellite markers and genotyped 15 populations ofL.sinensisusing nine newly characterized microsatellites. This NGS-based microsatellite development approach is a powerful technique for developing microsatellite markers in non-model species. Low genetic differentiation was detected between north and south populations, which suggested high genetic connectivity. Furthermore, populations on artificial structures located along Jiangsu Province showed a close genetic affinity with natural populations south of the Changjiang River estuary,suggesting their southern origin. To our knowledge,we provided here the first microsatellites database and population genetic analysis ofL.sinensisfrom almost the whole distribution range. This study can potentially enhance the understanding of population connectivity in the northwestern Pacific and provide insight into the ecological and evolutionary study ofLittoraria.

6 DATA AVAILABILITY STATEMENT

The raw whole genome shotgun sequencing data were deposited in NCBI SRA data set: SRR13107012.And the de novo assembly sequences are accessible via FigShare through https://figshare.com/s/7dcc021c5b25e3ebffd7. Other data generated or analysed during this study are included in this published article and its supplementary information files.