1 INTRODUCTION

Sea star (Echinodermata: Asteroidea) is an ancient group of marine invertebrates. Asterias amurensis Lütken, 1871 is a member of the family Asteriidae, originally inhabiting waters in the North Pacific, especially Russia, Japan, Korea, and China (Scherer, 2018). In China, A. amurensis is mainly distributed in the Yellow Sea and Bohai Sea. In addition, A. amurensis, as an invasive species, can also be found in southern Australia (Byrne et al., 1997). The gonads of A. amurensis are edible and rich in high-quality proteins, vitamins and lipids (Li et al., 2014). However, A. amurensis has been known as a severe "marine pest" due to high prey consumption. Their diverse diets may cause acute economic losses by preying on bivalves, such as mussels, scallops, and clams (Ross et al., 2002).

The massive and sporadic population outbreaks of A. amurensis have caused extensive economic and ecological losses to the local aquaculture industry and marine ecosystem over the past few decades. For example, outbreaks of A. amurensis devastated the shellfish industry in Japan (Ward and Andrew, 1995). In Tasmania, A. amurensis is the top predator in benthic communities, causing local extinctions of several species (Ross et al., 2004) and up to 50% of commercial scallop losses (Hutson et al., 2005). Similar cases also occurred in China. In Jiaozhou Bay, two massive outbreaks of A. amurensis were recorded in 2007 and 2021, respectively, which caused severe damage to the farming industry of Philippine clam (Ruditapes philippinarum). Several studies regarding the reproduction, development, and ecology of A. amurensis have been performed on the mechanism of population dynamics (Lee et al., 2004; Kashenko, 2005; Li et al., 2018). However, many questions remain unanswered. In particular, no study has been undertaken in A. amurensis to investigate the genetic diversity and relationships among populations. Understanding the genetic diversity and population structure of this species can also provide more-reliable estimates of population dynamics for resource management.

In recent years, mitochondrial DNA (mtDNA) and simple sequence repeat (SSR) microsatellite markers have been widely used to estimate the molecular variability, population genetic diversity, and genetic differentiation of many marine invertebrates (Cho et al., 2007; Meriam et al., 2015; Nie et al., 2015; Zheng et al., 2019; Zhang et al., 2020). The cytochrome oxidase subunit I (COI) gene is one of the most conserved protein-encoding genes in the mitochondrial genome. It is a widely used genetic marker for population genetic studies because of its fast evolution, maternal inheritance, and easy amplification and sequencing. In contrast, microsatellites are nuclear, biparentally inherited, highly polymorphic, and easy to isolate; therefore, they are particularly informative for population genetics studies. Overall, utilizing both fast-evolving and slow-evolving molecular markers together can comprehensively illustrate ancient and recent population histories. Thus, these two sets of molecular markers are complementary and indispensable in analyzing the population genetic structure of marine invertebrates.

To study the genetic diversity and population structure of A. amurensis, the samples of A. amurensis were collected from four geographic locations along the northern coast of China. A portion of the mitochondrial COI gene and ten SSR markers were employed to estimate the genetic diversity and population structure. The results will try to clarify two questions about A. amurensis, (ⅰ) whether there are differences in genetic diversity levels among different geographical populations; (ⅱ) whether there is genetic differentiation among starfish populations in and out of the outbreak area.

2 MATERIAL AND METHOD 2.1 Sampling

In 2021, a total of 120 individuals were collected by SCUBA diving from 4 locations (Qingdao (QD), Weihai (WH), Yantai (YT), and Dalian (DL); Fig. 1). Each location contains 30 individuals. Gonad tissue was kept in ethanol on site. Genomic DNA from gonad tissue was extracted by using the Plant Genomic DNA Kit (TIANGEN, Beijing, China). The integrity of genomic DNA was checked on 1% agarose gel electrophoresis, and the concentration was measured by using NanoDrop 1000 Spectrophotometer (NanoDrop, Wilmington, DE, USA).

2.2 COI amplification

The primers (AaCOI-F: 5'-GCTGGAACTGGCTGAACGAT 3'; AaCOI-R: 5'-TTCCCTGTAGGGTGGCCATT3') used for gene amplification were designed with reference to the complete mtDNA sequences of A. amurensis (NC_006665.1) and can amplify a 625-bp COI sequence. PCR amplification was carried out in a 50-μL reaction volume with predenaturation at 94 ℃ for 5 min, 35 cycles of denaturation at 94 ℃ for 45 s, annealing at 52 ℃ for 45 s, extension at 72 ℃ for 30 s, and further extension at 72 ℃ for 10 min. The PCR contained 25-µL 2×Taq plus MasterMix (Biosharp, China), 4-μL DNA template, 2-μL upstream and downstream primers (10 μmol/L), and 17-μL double-distilled water. The PCR products were checked by 1% agarose gel electrophoresis and then sent to Sangon Biotech Co. (Shanghai, China) for purification and sequencing.

2.3 SSR amplification

Ten SSR markers developed by Richardson et al. (2012) were used in this study. These 10 PCR sets contained the forward primer labeled with fluorescence dye (FAM) at the 5' end. PCR was carried out in a volume of 25 μL containing 1-µL DNA template, 0.5-μL dNTPs (10 μmol/L), 2.5-μL 10×PCR buffer, 0.5-μL upstream and downstream primers (10 μmol/L), 0.2-μL Taq DNA polymerase (Sangon, China) and 19.8-μL double-distilled water. Amplifications were performed by using touchdown PCR according to the following cycling program: denaturation at 94 ℃ for 30 s, annealing with temperatures decreasing 0.5 ℃ per cycle from 60 ℃ to 55 ℃ during the first 10 cycles and temperatures of 55 ℃ for the last 30 cycles for 30 s, and extension at 72 ℃ for 30 s. PCR products were separated and visualized using an ABI 3730XL automated sequencer by Sangon Biotech Co. (Shanghai, China), and the locus was genotyped with GeneMapper Software (Applied Biosystems, USA).

2.4 Statistical analysis

In the case of the COI sequence, the sequences were trimmed, edited and aligned in Bioedit 7.1.9.0 (Hall, 1999). Genetic diversity parameters were estimated in DnaSP 6.12.03 (Rozas et al., 2017). A haplotype network was further constructed using TCS method in PopART 1.7 software (Leigh and Bryant, 2015). Mismatch distributions analysis with the sudden expansion model (Rogers and Harpending, 1992) were performed using Arlequin 3.5.2.2 (Excoffier and Lischer, 2010). The calculation of pairwise F_st and analysis of molecular analysis of variance (AMOVA) were implemented in Arlequin 3.5.2.2 (Excoffier and Lischer, 2010). The significance of pairwise F_st values was determined through 1 000 permutations.

Based on SSR markers, the number of alleles (N_A), observed heterozygosity (H_o), expected heterozygosity (H_e), and inbreeding coefficient (F_is) were calculated using the R package hierfstat (Goudet, 2005). The polymorphic information content (PIC) and Hardy-Weinberg equilibrium (HWE) were calculated using the R package polysat (Clark and Jasieniuk, 2011) and pegas (Paradis, 2010), respectively. Considering the small sample size for each population, the contemporary effective population size (N_e) was only estimated for the overall populations by using NeEstimator V2.01 (Do et al., 2014). To evaluate the degree of genetic differentiation, the AMOVA was also implemented in Arlequin 3.5.2.2 (Excoffier and Lischer, 2010) and the pairwise F_st values were measured in the hierfstat package (Goudet, 2005) using the pairwise. neifst and boot.ppfst functions. The principal component analysis (PCA) was computed using the R package polysat based on Bruvo distance (Bruvo et al., 2004) and displayed by using the R package Scatterplot3d (Ligges and Mächler, 2003). In addition, the ancestry levels and admixture proportions were evaluated by using Bayesian clustering in STRUCTURE v2.3.4. The number of ancestral populations (K) ranging between 1 and 8 was tested, and 5 replicate analyses for each K were performed with a burn-in length of 20 000 and MCMC repetitions of 100 000. The optimal K was determined using the DeltaK method through Structure Harvester v.0.6.94.

3 RESULT 3.1 Genetic diversity based on COI sequence

From the initial 625 bp COI fragment amplified, a fragment of 543 bp of the A. amurensis COI gene was used after alignment and trimming. The average nucleotide frequencies of A, T, G, and C were 37.4%, 22.7%, 21.0%, and 18.9%, respectively. The genetic diversity parameters are presented in Table 1. A total of 42 polymorphic sites were detected. The nucleotide diversity (π) ranged from 0.004 (WH) to 0.007 (QD), with an average of 0.006. A total of 36 haplotypes were identified, and one main haplotype was found in four populations. The lowest value of haplotype diversity (H_d) was detected in YT (0.745), and the highest value was detected in QD (0.924). The haplotype network showed a high level of complexity, and no obvious spatial cluster pattern was found (Fig. 2). Mismatch analysis of all the samples showed unimodal distributions (Fig. 3).

3.2 Genetic structure based on COI sequence

The analysis of molecular variation (AMOVA) showed that 97.47% of the total molecular variance was distributed within populations, whereas 2.53% of the total molecular variance occurred among populations (Table 2). The average F_st value was 0.025 3 (P < 0.05), indicating a significant genetic variation among these four populations. The pairwise F_st ranged from -0.001 to 0.049 (Table 3). Of the 6 pairwise comparisons among all locations, 2 showed significant values of F_st.

3.3 Genetic diversity based on SSR markers

The diversity indices for ten SSRs were shown in Table 4. All the loci were polymorphic in each population. The number of alleles each locus ranged from a minimum of 6 (Aamr02 in QD, DL, and WH, and Aamr04 in WH) to a maximum of 22 (Aamr31 in YT) with an average count of 11.8 alleles. The H_o per locus and H_e ranged from 0.433 to 0.933 and 0.557 to 0.924, respectively. The F_is had a wide variation among combinations of population and SSR, with values from -0.058 to 0.515. The PIC ranged from 0.505 (Aamr02 in DL) to 0.902 (Aamr18 in QD). Eight significant deviations from HWE (P < 0.001) were observed in the 40 combinations of populations and SSR markers. Besides, the estimated N_e was 349.9 for all populations.

3.4 Genetic structure based on SSR markers

The AMOVA based on SSR presented that 98.94% of the total molecular variance was distributed within populations, and 1.06% of the total molecular variance occurred among populations (Table 2). The pairwise F_st values ranged from -0.001 to 0.017 (Table 3), and the highest level of differentiation was observed between QD and WH. Principal component analysis (PCA) derived from the SSR dataset revealed that most individuals clustered together, which was not related to the geographical pattern (Fig. 4). The admixture analysis presented similar results. Using Structure Harvester, the maximum DeltaK value was identified when K=3. However, all individuals had a similar admixture (Fig. 5).

4 DISCUSSION

Genetic diversity within and between populations could provide potential genetic resources for future adaptation, which is critical for the fitness of a population (Liu and Zhou, 2017). The current study is the first report on the population genetic diversity of A. amurensis in China seas. The results showed that all populations displayed high levels of haplotype diversity (H_d=0.745–0.924) and low levels of nucleotide diversity (π=0.004–0.007) based on the COI sequence. This pattern is consistent with previous studies about starfish, including Echinaster sepositus (Garcia-Cisneros et al., 2016) and Acanthaster spp. (Chen et al., 2021). Previously, Grant and Bowen (1998) classified marine fishes into four categories according to different combinations of H_d and π. According to their hypothesis, A. amurensis belongs to the second category in the current study. In this condition, A. amurensis in China seas had experienced rapid expansion and population growth after a period of low effective population size in their history. This can be confirmed by the mismatch distribution analysis. If populations had undergone demographic expansion, the distribution of pairwise differences between sequences would be unimodal (Rogers and Harpending, 1992).

Across the four geographical populations, the QD population displayed the highest values of H_d and π. Similar results can also be achieved from the diversity indices based on SSR markers. For example, the QD population displayed the highest value for PIC. These patterns might benefit from the rapid population expansion in Jiaozhou Bay, which enhanced the retention of new mutations in the QD population. In addition, all four populations showed high positive F_is values. Generally, positive F_is values are mainly caused by mating among relatives. The N_e value estimated for A. amurensis was lower than that of many aquatic species, including Fenneropenaeus chinensis (Wang et al., 2016) and Oncorhynchus tshawytscha (Larson et al., 2014). This suggested that all individuals might originate from a few ancestors, which may lead to inbreeding. Another factor may be asexual reproduction by the regeneration of arms attached to a portion of the central disc (Ward and Andrew, 1995). Furthermore, the life history trait of A. amurensis with free-spawned planktonic sperm can also contribute to the high positive F_is, considering that the F_is values were usually higher in the species with planktonic sperm than in those species that have direct sperm transfer to females or copulate (Addison and Hart, 2005).

The AMOVA and pairwise F_st values were firstly performed using COI sequences and SSR markers to assess the genetic differences among populations. Both results displayed significant genetic differentiation among populations, especially between QD and WH. Wright (1978) suggested that the F_st values ranging from 0 and 0.05 correspond to a little genetic differentiation. Hence, the current study may indicate that little genetic differentiation has developed between QD and WH. The current genetic difference may be shaped by environmental factors including ocean temperature or nutrient concentrations (Fontaine et al., 2007; Chen et al., 2021), but may also result from random processes in the transmission of alleles from one generation to the next. Meanwhile, it' s worth noting that the current F_st value estimated by using SSR may be deflated due to the high degree of polymorphism and higher mutation rates of microsatellites (Balloux and Lugon-Moulin, 2002). Recently, genome-wide single nucleotide polymorphisms (SNPs) had been developed to evaluate the genetic structure of Pacific crown-of-thorns starfish (Iguchi et al., 2021; Horoiwa et al., 2022). Hence, population genomics approach based on genome-wide SNPs and copy number variations (CNVs) may provide powerful candidate to evaluate more systematically population genetic parameters of A. amurensis in the future.

In the current study, the PCA and admixture analysis showed no apparent spatial cluster pattern. Still, several individuals from YT and DL have little genetic association with other individuals. A similar phenomenon was also recorded in previous studies. For example, some seahorses from New Zealand had no similar ancestries to local populations. Interestingly, these individuals displayed the same ancestry as the seahorses from the Australian Sea (Ashe and Wilson, 2020). Similarly, two Pacific cod from coastal areas displayed the same ancestry as the fish from the Salish Sea (Drinan et al., 2018). Therefore, in the current study, several differentiated individuals from YT and WH may imply that there is a genetically differentiated population in adjacent sea. In particular, it is worth exploring whether the population genetic structure of A. amurensis from the Bering Sea to the Japan Sea is different from that of the Chinese seas, considering that the two sea areas are separated by the Korean Peninsula, which might form a dispersal barrier. Thus, investigations on a larger spatial scale should be carried out to evaluate the population genetic structure of A. amurensis in the North Pacific region.

5 CONCLUSION

Having evaluated the genetic diversity and population structure of A. amurensis from four populations along the northern coast of China by analyzing COI sequence and SSR markers, a total of 36 haplotypes were identified, and a main haplotype was found in four populations. The QD population displayed the highest genetic diversity. A small but statistically significant population differentiation were observed among the four populations, especially between QD and WH, by the AMOVA and pairwise F_st. In addition, the PCA and admixture analysis identified several individuals in YT and DL with little genetic association with other individuals. Overall, the current results provide useful information on the genetic diversity and population structure of A. amurensis and also will contribute to the resource management of A. amurensis in China.

6 DATA AVAILABILITY STATEMENT

The COI sequence supporting the findings of this study are openly available in GenBank (reference numbers: ON171504–ON171623), and the SSR dataset are available from the corresponding author upon request.