Development of simple sequence repeat markers and genetic diversity of Sipunculus nudus in the Beibu Gulf of China

Materials

A total of 219 samples of S. nudus were collected from eight geographic locations, namely Tieshangang District, Beihai City (TSG, 27 samples); Yingpan Town, Beihai City (YP, 28 samples); Yinhai District, Beihai City (BH, 27 samples); Hepu County, Beihai City (HP, 27 samples); Rhinoceros Horn Town, Qinzhou City (QZ, 29 samples); Guangpo Town, Fangchenggang City (FCG, 27 samples); Caotan Town, Zhanjiang City (ZJ, 27 samples) and Guangcun Town, Danzhou City, Hainan Province (HN, 27 samples) (Table 1; Fig. 1). All samples were transported back to the laboratory. The body wall of each sample was washed with PBS and used to extract genomic DNA via TIANamp Marine Animals DNA Kit (DP324, Tiangen Biotech; Beijing; China).

Development of genome-wide SSRs

The MIcroSAtellite (MISA) tool was utilized to identify SSR loci in the genome deposited in the Genome Warehouse in National Genomics Data Center (BioSample: SAMC4151241). The criteria for SSRs in this study were as follows: a minimum of 20, 12, six, five, four and three repeat units for one to six nucleotide motifs, respectively (Han et al., 2024). The primer3 was employed to design the corresponding primer pairs by default parameters.

A hundred primer pairs (15–18 primer pairs were selected for each motif type, with five to seven pairs evenly distributed on each chromosome) were selected and synthesized by Sangon Biotech (Shanghai) Co., Ltd.

PCR amplification and polymorphic SSRs screening

The synthesized SSR primer pairs were screened using six random DNA samples as templates. PCR amplification was performed in a volume of 10 µL, which contained 2 × SanTaq PCR Mix (B532061; Sangon Biotech (Shanghai) Co., Ltd.) five µL, 0.2 µL each of forward and reverse SSR primers (10 µmol L⁻¹), one µL of DNA template and 3.6 µL of ddH₂O (Han et al., 2024). The 2 × SanTaq PCR Mix contains Mg²⁺, dNTPs, TaqDNA polymerase, PCR buffer, loading, and PCR enhancer.

The amplification program was as follows: pre-denaturation at 95 °C for 4 min; 30 cycles of 95 °C for 45 s, the annealing Tm (depending on primers) for 45 s and extension at 72 °C for 45 s; followed by a final extension of 7 min at 72 °C.

The amplification products were separated by electrophoresis on 6% non-denaturing polyacrylamide gels and visualized by silver staining according the previous report (Han et al., 2024). The bands were recorded by hand according the molecular weight for future analysis.

Genetic diversity

GeneAlEx6.501 (Peakall & Smouse, 2012) was employed to calculate the parameters of genetic diversity including the number of alleles (Na), the effective number of alleles (Ne), expected heterozygosity (He), observed heterozygosity (Ho), Shannon’s index (I), and fixed index (F). Additionally, Cervus V3.0 software was used to calculate the polymorphic information content (PIC) (Kalinowski, Taper & Marshall, 2007).

Genetic distances and gene flows

GeneAIEx 6.501 was also used to assess the parameters of population genetics including gene flow (Nm), genetic differentiation (Fst), genetic distance (Nei’s) and genetic similarity (Nei’I) among the populations. Molecular AMOVA was performed using Arlequin V3.5.2.2 software (Excoffier & Lischer, 2010) to detect the distribution of genetic variation within and among populations of S. nudus.

DARwin6 was employed to generate a phylogenetic tree file in the format of PAU via neighbor-joining method with 1,000 bootstrap replications. The clustering dendrogram was drawn using the iTOL web platform.

Genetic structure

Structure 2.3.4 software (Sheehan et al., 2001) was utilized to infer genetic structure using the admixture model. A burn-in period of 10,000 and Markov chain Monte Carlo (MCMC) replicates of 100,000 were set. Initially, the number of subpopulations (K) ranged from 2 to 10. The online software Structure Harvest was employed to determine the theoretical population number based on the ΔK algorithm (Evanno, Regnaut & Goudet, 2005).

Distribution and characteristics of genome-wide SSRs

A total of 277,264 SSRs (273,552 SSRs on chromosomes and 3,712 SSRs on scanffolds) was identified in the genome of S. nudus (1,466.684 Mb) (Fig. 2). Of them, there were 51,025 compound SSR loci. In the genome, dinucleotide motifs (213,846, 77.127%) predominated (Table 2), while pentanucleotide motifs (1,584, 0.571%) were the least. The average length of the pentanucleotides was the longest (73.045 bp), whereas that of the mononucleotides was the shortest (22.722 bp). The density was 189.041 SSRs/Mb in the genome.

There were 175,735 SSRs consisting of motifs more than 16 times, accounting for 63.382%. The SSRs with the motifs of nine times were relatively fewer (Table 3).

There was a total of 135 different motifs in the genome. Of them, the AT/AT motif dominated (170,807, 61.604%).

Screening of polymorphic SSRs

A total of 198,902 primer pairs were successfully designed. Among them, 100 primer pairs were selected for synthesis and screened. The results showed that five primer pairs did not amplify bands, while 95 primer pairs could amplify bands including two primer pairs with single band (Fig. 3). This suggested that the successful rate of amplification was 95% and polymorphic rate was 93%.

Genetic diversity

Thirty highly polymorphic primer pairs were employed to evaluate the genetic diversity of 219 samples of S. nudus. The results (Table 4) indicated that 370 alleles were detected and the Na ranged from five to 23 with an average of 12.333; the Ne varied from 2.432 to 13.743 with an average of 6.950; the Ho ranged from 0.758 to 1.000 with its average of 0.905; the He ranged from 0.589 to 0.927 with an average of 0.833; the I ranged from 1.050 to 2.783 with an average of 2.042; the PIC ranged from 0.506 to 0.923, with an average of 0.811. All these findings suggested that the genetic diversity of the population was relatively high in this study.

As for the genetic diversity of eight populations, the results (Table 5) indicated that the ZJ population had the maximum genetic diversity with the parameters of 6.233, 4.271 and 1.479 in Na, Ne and I. In contrast, the YP population had the minimum values of 1.933, 1.596, 0.558, respectively. The Ho ranged from 0.382 (YP) to 0.94 (QZ) with a mean of 0.790. The He ranged from 0.312 (YP) to 0.718 (ZJ) with the mean of 0.590. The values of Ho were more than that of He in all population. All the fixation indices (F) of the populations were negative with the mean of −0.373. According to the He values, the eight populations in descending order were YP, HN, QZ, HP, FCG, BH, TSG, and ZJ.

Analysis of molecular variance analysis

The results of Analysis of molecular variance (AMOVA) analysis showed that 3.830% of the genetic variation existed among populations, while 96.170% were observed within the population. The genetic variation mainly existed within the population rather than that between the populations. The fixation index was 0.0383.

Genetic differentiation

The study revealed genetic distances between populations ranging from 0.262 to 2.998 (average 1.310), with the largest divergence between HP and HN and the smallest between HP and ZJ. Genetic similarity values spanned 0.050–0.769 (average 0.332), showing HP’s lowest similarity with HN and highest with ZJ. For S. nudus, pairwise Fst values ranged from 0.049 to 0.581 (average 0.239), with the lowest Fst between HP and ZJ and the highest between YP and HN. Nm values varied from 0.180 to 4.842 (average 1.434), peaking between HP and ZJ and reaching the minimum between YP and HN.

The clustering results (Fig. 4) indicated that the samples of eight populations of S. nudus were divided into five groups. Group 1 contained the samples from QZ, TSG and FCG; Group 2 consisted of the individuals from HP, ZJ and BH; Group 3 included the samples from YP; Group 4 harbored the individuals from HN. Group 5 only had a single sample from TSG. In Group 1, two samples from QZ appeared in samples from FCG; similarly, two individuals from HP clustered with samples from ZJ.

Genetic structure

The results demonstrated that the K with the maximum of the ΔK at the turning point was seven. Evidently, the whole population was divided into seven subpopulations in this study (Fig. 5). Subpopulation 1 included 12 individuals (TSG1∼TSG12, green); Subpopulation 2 consisted of 15 individuals (TSG13∼TSG27, orange); the subpopulation 3 (red, 56 individuals) was composed of 27 individuals (FCG1-FCG27) and 29 individuals (QZ1-QZ29); Subpopulation 4 had 22 individuals (YP1∼YP15, YP18, YP20∼YP24, and YP28, blue); Subpopulation 5 harbored 33 samples (YP16, YP17, YP19, YP25, YP26, YP27, and BH 1∼27, yellow); Subpopulation 6 (violet, 54) consisted of 27 samples from HP (HP1-27) and 27 samples from ZJ (ZJ1-27); and the samples of Subpopulation 7 come from HN (azure, HN1-27).

Genome-wide SSRs

Now, various types of molecular markers have been extensively utilize to the identification of marine animal genetic resources (Kong, Li & Qiu, 2007). To date, only few types of molecular markers including RAPD (rapid amplification polymorphic DNA) markers, DNA fragments on mitochondrial genome and SSRs (Wang et al., 2012) were used to evaluate the genetic diversity in S. nudus. However, RAPD markers had some limitations including low instability, and inability of identifying heterozygotes, while conserved mitochondrial genes offered insufficient polymorphism. High polymorphism and codominance are key features of ideal molecular markers. SSR markers, meeting these criteria, are widely regarded as the most practical. However, only few SSRs have been developed for S. nudus (Guo et al., 2012; Wang et al., 2012), which significantly limit molecular marker-aided breeding. To address this, large-scale SSR markers were developed in the present study.

The characterization of genome-wide SSRs revealed the presence of various repeat motifs in the genome. Among them, dinucleotide dominated (70.18%) which was consistent with the SSRs in most aquatic animal genomes such as Protosalanx chinensis (Tang & Zhou, 2023), Schizothorax biddulphi (Ren, 2020), and Puffer fish (Edwards et al., 1998). Further studies revealed that the AT/AT was predominant in the genome of S. nudus, while AC/TG was more frequent in the genome of the Schizothorax biddulphi. This might be due to the differences of species and SSR criteria.

The distribution of SSRs was non-uniform, with regional variations in density and sporadic gaps detected across all chromosomes, likely attributable to taxonomic differences in SSR composition (Srivastava et al., 2019).

Usability of SSRs markers

This study developed numerous codominant SSR markers and validated 100 in S. nudus. Of these, 95 primer pairs amplified products, with 93 exhibiting polymorphism, demonstrating their utility for assessing genetic diversity and enabling marker-assisted breeding in S. nudus.

Genetic diversity

Genetic diversity, the fundamental basis of populations, forms the foundation for both environmental adaptation and genetic improvement (Cui et al., 2012).

The 30 polymorphic primer pairs showed Na and Ne values of 12.333 and 6.950, respectively, in the entire population—approximately four-fold higher than those reported (Wang et al., 2012) and three-fold greater than values documented (Guo et al., 2012). The average of PIC was 0.811, which was higher than that in other reports. All the results above suggested that the genetic diversity here was abundant (Botstein et al., 1980). In addition, the Shannon’s index of the 30 polymorphic primer pairs in entire population were 2.042, also surpassed 0.269 and 0.394 reported previously (Wang, Du & Li, 2006; Song et al., 2011), which may be closely related with the sample size and the types of molecular markers in another part. The results demonstrated higher genetic diversity in the population of S. nudus in Beibu Gulf.

In the present study, the Ho s were significantly higher than He s, which was consistent with the report (Wang et al., 2012). Both Ho and He were relatively similar and lower than those in this study (Guo et al., 2012). The observed heterozygosity exceeded the expected value, suggesting selection against homozygotes. The higher proportion of heterozygous individuals, carrying two distinct alleles with different functional roles, enhanced the population’s adaptive capacity to environmental changes (Jiang & Ma, 2014).

Genetic differentiation and gene flows

Genetic differentiation arises from restricted gene flow due to reproductive and geographic isolation, with key indicators like Fst and Nm quantifying divergence levels (Ren, 2020). Genetic differentiation levels are categorized based on Fst thresholds: little to none (0 <Fst <0.05), moderate (0.05–0.15) and high (Fst >0.15) between populations (Wright & Curnow, 1978). Therefore, moderate genetic differentiation existed among 10 pairwise populations, whose Fst s ranged from 0.05 to 0.15. The higher genetic divergences were among 17 pairwise populations. The genetic divergence was minimal (0.049) between ZJ and HP, while maximum (0.581) between YP and HN. The spatial segregation between the northern YP and southeastern HN populations in the Beibu Gulf restricted gene flow, leading to greater genetic differentiation. This has been confirmed (Kawauchi & Giribet, 2014). However, subpopulations were rarely fully isolated (occasional gene flow remains possible).

Genetic structure

The genetic composition of a population indicates its diversity, which in turn determines a species’ ability to adapt to environmental changes (Niu et al., 2019). Analyzing population structure through allele frequency patterns helps trace evolutionary origins, identify distinct lineages, and reveal genetic exchange between populations (Liu et al., 2017).

This study grouped the entire population into seven distinct genetic groups, demonstrating clear genetic structure. Moreover, the identified genetic clusters closely matched their predefined classifications.

The genetic structure is influenced by many factors including mating system, gene flow, selection and genetic drift et al. (Wang et al., 2017; Hu, Chen & Francis, 2021). Mating systems play a key role in shaping genetic structure. For S. nudus, sperm and eggs are released into seawater, where they are mixed evenly with current waves. This result in random fusion—a hallmark of random mating. hence, we speculated that the mating system did not drive the genetic structure in this species. This study found that both Subpopulation 5 and Subpopulation 6 consisting of various provenances showed signals of gene flow. Gene flow occurs via gamete dispersal (e.g., ocean currents) or human-mediated transfer. Restricted gene flow over time drives population divergence, while intense market demand for S. nudus (food/medicine) drives overharvesting, causing population declines and potential genetic structure alterations. Our findings suggested human activities were the primary driver of the observed genetic differentiation in this species. The observed genetic structure in S. nudus likely stemmed from human-driven overexploitation, which accelerated population decline and directly reduced genetic diversity via shifts in gene frequencies (Cheng, Kao & Dong, 2020). Higher observed than expected heterozygosity indicated an excess of heterozygotes. As heterozygotes carry two distinct alleles with potentially divergent functions, this suggests that populations with elevated heterozygosity may exhibit greater adaptive capacity in complex environments. Therefore, this entire population have greater adaptive capacity (Jiang & Ma, 2014).