De novo transcriptome assembly and population genetic analyses of an important coastal shrub, Apocynum venetum L

Background Apocynum venetum L. is an important medicinal plant that is mainly distributed in the coastal areas and northwest of China. In addition to its high medical and economic value, its adaptation to saline-alkali and coastal saline lands makes A. venetum an ideal candidate for use in vegetation restoration. To date, the study of A. venetum has been limited in the northwest region of China, little attention has been paid to the genetic diversity and population structure of A. venetum populations in the coastal region. Here, we performed transcriptome sequencing of total RNA from A. venetum leaves and developed efficient expressed sequence tag-simple sequence repeat (EST-SSR) markers for analyzing the genetic diversity and population structure of A. venetum in the coastal region. Results A total of 86,890 unigenes were generated after de novo assembly, and 68,751 of which were successfully annotated by searching against seven protein databases. Furthermore, 14,072 EST-SSR loci were detected and 10,243 primer pairs were successfully designed from these loci. One hundred primer pairs were randomly selected and synthesized, twelve primer pairs were identified as highly polymorphic and further used for population genetic analysis. Population genetic analyses showed that A. venetum exhibited low level of genetic diversity (mean alleles per locus, NA = 3.3; mean expected heterozygosity, HE = 0.342) and moderate level of genetic differentiation among the populations (genetic differentiation index, FST = 0.032–0.220) in the coastal region. Although the contemporary (mean mc = 0.056) and historical (mean mh = 0.106) migration rates among the six A. venetum populations were moderate, a decreasing trend over the last few generations was detected. Bayesian structure analysis clustered six populations into two major groups, and genetic bottlenecks were found to have occurred in two populations (QG, BH). Conclusions Using novel EST-SSR markers, we evaluated the genetic variation of A. venetum in the coastal region and determined conservation priorities based on these findings. The large dataset of unigenes and SSRs identified in our study, combining samples from a broader range, will support further research on the conservation and evolution of this important coastal plant and its related species.


Background
Coastal habitats, which are located in the transition zone between terrestrial ecosystems and marine ecosystems, are characterized by unique, complex ecosystems with high ecological value [1]. Their unique environmental features, including poor soil structure with low water infiltration and poor drainage, a high salt content and PH level, give rise to unique plant diversity and many specialist species [2]. However, due to the rapid urbanization process, coastal habitats have been particularly affected by both land transformation and mass tourism, leading to the severe disturbance and loss of natural habitats [3]. Many coastal plants are adapted to specific costal environmental conditions and are thus highly vulnerable to habitat destruction [4][5][6]. Therefore, there is an urgent need to evaluate and protect the irreplaceable, vulnerable biodiversity found in the coastal zones, especially in the face of continuing anthropogenic pressure.
The coastal area of Jiangsu Province accounts for one quarter of the total coastal area in China, which mainly falls within Yancheng City [7]. The coastline of Yancheng City's is 582 km long, and the beach area covers 6.83 million hm 2 , with an annual growth rate of 10,000 hm 2 [7]. The Yancheng coastal region not only provides precious land resources, but is also an important wetland nature reserve. According to the plant surveys, there are 688 kinds of vascular plants belonging to 391 genera and 114 families in the Yancheng tidal flat wetland [8]. Multiple large-scale beach reclamation projects have been conducted in Yancheng city since 1949, and the natural habitats of this region have been significantly modified [7]. Many studies have shown that anthropogenic habitat alterations can affect both global biodiversity and genetic diversity, which jeopardizes the longterm survival of species and increases their risk of extinction [9][10][11][12][13]. There are many valuable plant resources, such as Limonium sinense (Girard) Kuntze, Apocynum venetum, Tamarix chinensis Lour., Tournefortia sibirica, and Salicornia europaea L. etc. in the Yancheng tidal flat, which are not only tolerant to salt stress but also present high medical and economic values [8]. Whether the genetic diversity and genetic structure of coastal plant populations have been affected by coastal habitat changes has seldom been evaluated.
Apocynum venetum L. is a perennial shrub that is specifically distributed in the coastal region of Jiangsu Province. It has also been referred to as "Luobuma" since it was first discovered on the Luopu plains of Xinjiang Province in the 1950s [14]. A. venetum not only provides precious fiber and nectar resources but is also used as an important medicinal plant for treating hypertension and hyperlipidemia. Its high stress resistance to high salt contents and poor soils also contribute to its great ecological value. To date, studies of A. venetum have mainly focused on its medicinal effects and physiological characteristics such as photosynthesis and water absorption [15][16][17]. However, few studies have examined the genetic diversity and genetic structure of natural A. venetum populations. To better protect and utilize this important plant species in the Yancheng coastal region, there is an urgent need to evaluate the genetic diversity and population structure of natural A. venetum populations.
Simple sequence repeat (SSR) markers have been widely used in population genetic analysis and molecular marker-assisted breeding because of their high polymorphism, repeatability and codominant inheritance [18,19]. A number of transcriptome analyses have been conducted on A. venetum [20][21][22], while SSR markers have not yet been developed and used to evaluate the genetic diversity of natural A. venetum populations yet. Using RAPD (Random Amplification Polymorphic DNA) and AFLP (Amplified Fragment Length Polymorphism) markers, researchers have detected moderate to high levels of genetic diversity in A. venetum populations from Xinjiang and Inner Mongolia regions [23,24]. Except in the arid region, there has been no study that has evaluated the genetic diversity of A. venetum populations in the coastal regions. Therefore, in this study, we conducted comprehensive transcriptome sequencing of A. venetum from the coastal region of Jiangsu Province using the Illumina sequencing platform. After data assembly and annotation, we developed a set of novel EST-SSR markers from the unigenes. By using these informative markers, we successfully evaluated the genetic diversity, population structure and demographic changes in six populations across the natural distribution of A. venetum in the Yancheng coastal region. We expect that the genetic information identified in this study will facilitate the management and conservation of natural A. venetum populations in the future.

Assembly of A. venetum transcriptome data from Illumina sequencing
Transcriptome sequencing generated 46,408,308 reads, totalling approximately 6.96 Gb for A. venetum in this study. After stringent quality filtering, 45,760,331 (98.6%) high-quality reads were obtained, exhibiting 98.31% Q20 bases and a GC value of 46.93% (Table 1). A total of 86,890 unigenes were successfully assembled using the Trinity software, with a mean length of 1767 bp and an N50 of 2580 bp. Among all the assembled unigenes, 3119 of which (approximately 3.58%) were less than 300 bp, and 14,657 unigenes (16.85%) were longer than 3000 bp, whereas most of the unigenes (69,204) (79.56%) ranged from 300 to 3000 bp (Additional file 1: Fig. S1A). The number of reads mapped to each unigene analysis revealed that 13,065 unigenes (about 47.90%) and 5506 unigenes (approximately 20.19%) consist of more than 100 and 1000 reads each, respectively, only a few unigenes (3.3%) were derived from less than 10 reads (Additional file 1: Fig. S1B).

Functional classification by GO and KOG
To further evaluate the functions of these unigenes, we used GO assignments to annotate and analyze each unigene. A total of 51,733 unigenes were assigned to 54 functional subgroups. Among the three ontology categories, the largest was biological process (47.6%), followed by cellular component (29.5%) and molecular function (22.9%) (Fig. 1). In the biological process group, the most frequent category was cellular process (22.04%), followed by metabolic process (20.99%). Cell (19.40%) and cell part (19.40%) were the most highly represented groups in the cellular component category. For the molecular function category, binding (45.49%) and catalytic activity (38.16%) accounted for the greatest proportions. Then we submitted all the assembled unigenes to the KOG database for further functional prediction and classification. Among the 25 KOG categories, posttranslational modification, protein turnover and chaperones consist the largest group (13.36%), general function prediction only (12.71%) and translation, ribosomal structure and biogenesis (10.60%) also showed high percentages (Additional file 2: Fig. S2).
Within the 14,072 SSRs, 10,243 primer pairs were successfully designed. We randomly selected 100 pairs from these primer pairs for amplification, and 35 were successfully amplified at expected sizes. Using twelve individuals from four A. venetum populations, these 35 primer pairs were applied to screen for polymorphism and twelve showed allelic polymorphism (Table 3). Using these 12 polymorphic EST-SSR markers, a total of 39 alleles were detected across the 103 samples from the six coastal A. venetum populations, with 2 to 5 alleles per locus. Using the Micro-Checker program, we found no evidence of null alleles and scoring errors in the dataset. Linkage disequilibrium analysis was performed between each pair of loci in each population. The results showed that 6 of 396 comparisons (Av39 and Av88 in population LD; Av02 and Av88 in population XY; Av08 and Av39, Av08 and Av63, Av39 and Av63, Av39 and Av75 in population SY) were significant after Bonferroni correction (P = 0.00013). Given that there is no overlap between these six pairs of loci, we treated these 12 EST-SSR markers as independent loci in the following analyses. For all loci, the expected (H E ) and observed heterozygosities (H O ) ranged from 0.030 to 0.651 and from 0.030 to 0.653, respectively. The polymorphism information content (PIC) was between 0.029 (Av55) and 0.575 (Av39), with an average of 0.297 (Table 3). Among 72 population-by-locus tests, departure from HWE was observed at locus Av02 in the LD population, loci Av08 and Av39 in the QG population after sequential Bonferroni correction (P < 0.004). Due to the significant P values of these three loci that were only present in a single population, their departures from HWE most likely reflect population-specific rather than locusspecific problems.  (Table 4). The F ST values of population pairs ranged from 0.032 to 0.220 (Table 5), with an overall value of 0.101, suggesting low to moderate levels of genetic differentiation across all the populations. Thirteen of the 15 pairwise comparisons were significant after sequential Bonferroni correction (Table 5).

Analyses of population genetic diversity and structure
In the structure analysis, we used LnP(D) and ΔK statistics to determine the most likely value of population genetic cluster K. Because the LnP(D) increased progressively from K = 1 to 6, it was difficult to determine the true number of genetic clusters (K). However, the ΔK statistic of Evanno et al. (2005) detected the highest peak at K = 2 ( Fig. 3a, b). Figure 3c exhibits the assignment of individuals to each cluster ('red' and 'blue' represent cluster I and cluster II, respectively) when K = 2. The individuals that contain both colors represent the mixed origin of two gene pools. Cluster I and II were present at a high frequency (78 and 88%) (Q > 0.80) in the DF and SY populations. For the QG and XY populations, nearly half of the individuals were assigned to each cluster respectively. For the LD and BH populations, 10 to 50% of all local samples were assigned to each cluster, and the remaining 30% consisted of mixed individuals. With increased values of K (K = 3), we observed that the SY population was assigned to a separated cluster III ('yellow') ( Fig. 3c), suggesting the potential genetic differentiation of this population from the rest of the populations in the Yancheng region.

Historical and contemporary gene flow
The MIGRATE-N results revealed that the mutation-scaled effective population size (Θ) for the six A. venetum populations ranged from 0.0342 to 0.0438. Historical gene flow (m h ) that inferred from the mutation-scaled migration rate (M) was highest from population LD to QG (m h = 0.169), and lowest from population LD to DF (m h = 0.041), with an average value of 0.106 across all the populations (Additional file 3: Table S1). Asymmetric gene flow was observed in two pairs of populations, with the predominant direction of gene flow occurring from population DF to LD and XY to BH. Using BAYE-SASS software, we found that the six populations were largely composed of individuals (the average 72%) that originated from within the same site, while approximately 6% of the individuals were exchanged with each other site (Additional file 3: Table S1). Although a moderate level of gene flow was detected among the six A. venetum populations, the Wilcoxon signed rank test indicated that the contemporary estimates (mean m c = 0.056) were significantly lower than the historical estimates of migration rates (mean m h = 0.106, P < 0.001),

Changes of effective population size
Under both the stepwise mutation model (SMM) and two-phase mutation model (TPM), our Wilcoxon test detected no significant heterozygote excess for most of the populations (Table 6). For the BH population, we detected an evidence of a historical bottleneck (P < 0.05; Table 6). Likewise, the mode-shift test revealed that most populations showed an L-shaped distribution of alleles, suggesting the absence of a recent bottleneck. The observation of a shifted mode in the two populations (QG, BH) suggested the occurrence of bottleneck events over the last few generations (Table 6).

Discussion
Characterization of A. venetum transcriptome and its potential use in germplasm resources evaluation With the decreasing costs of sequencing, the development of molecular markers that based on nextgeneration sequencing technology has become the most efficient method for molecular studies of non-model plants [25,26]. Using the Illumina sequencing platform, we acquired a well-assembled transcriptome sequencing data and developed a set of novel EST-SSR markers for A. venetum. As an important medicinal species and source of fiber, A. venetum has been the subject of extensive genetics and pharmacology studies. So far, ISSR (Inter-Simple Sequence Repeat), AFLP and RAPD molecular markers have been developed for A. venetum [23,24,27]. In recent years, a number of transcriptome analyses have been conducted on A. venetum [20][21][22]. Using leaf material, we acquired a similar amount of   A. venetum has been used to lower blood pressure and lipemia as a traditional medicine for a long time in China, and the roasted leaves of A. venetum have been commercialized as a sedative and anti-ageing supplement. A. venetum is rich in flavonoid compounds, which provide its broad pharmacological activities [31,32]. Flavonoids are synthesized through a long, complex pathway [33]. Identifying the regulatory mechanisms underlying flavonoid biosynthesis is essential for understanding the chemical composition or pharmacological activities of Apocynum. Previous studies have shown that A. venetum has a higher flavonoid content than its related species Apocynum. hendersonii, and some flavonoid components such as hyperoside have been identified as suitable chemical markers for the discrimination of the two species [34]. In a recent metabolome and transcriptome analyses of Apocynum, Gao et al. [21] found that the flavonoid biosynthetic pathway is responsible for a considerable proportion of the diversity between A. venetum and A. hendersonii, and identified anthocyanin as the key component that determines the phenotypic diversity of the stem and leaf color of A. venetum and A.
hendersonii. In our KEGG analysis, 28,505 unigenes were clustered into 130 pathways, including the flavone and flavonol biosynthesis (ko00944), flavonoid biosynthesis (ko00941), isoflavonoid biosynthesis (ko00943), and anthocyanin biosynthesis (ko00942) pathways. Among the 88 unigenes related to these pathways, 84 unigenes have an average FPKM (Fragments per kilobase of exon per million fragments mapped) value greater than one and 24 of these were above 50. In addition, 23 out of these 88 unigenes contain SSR loci (Additional file 5: Table S2). For example, unigenes that annotated as flavonol synthase (NR ID: BAD34463.1) and flavonoid 3-O-galactosyltransferase (NR ID: BAF49284.1), both contained SSRs and exhibited an average FPKM value of 50.38 and 22.16, respectively. Although ITS and cpDNA sequences have been used to identify A. venetum and its related species, these sequences are conservative within the genus and only provide limited genetic information [35]. The flavonoid biosynthesis related transcripts that contain SSRs identified in our study will serve as good candidates for the development of novel molecular markers in the future, which might aid in both the species identification in the Apocynum genus and the selection of A. venetum germplasm resources with high bast fiber and flavonoid content.

Low level of genetic diversity in A. venetum coastal populations
The disturbance of natural habitats poses a major threat to global biodiversity [36]. The ongoing process of habitat alteration generally has strong negative impacts on the species composition and genetic diversity of species [37][38][39]. The genetic evaluation of natural populations is an essential step in the conservation and utilization of plant resources, especially for vulnerable and threatened species. Using RAPD and AFLP markers, researchers have detected moderate to high levels of genetic diversity in natural A. venetum populations from the Xinjiang and Inner Mongolia regions [23,24]. They have also found that the genetic diversity of A. venetum is significantly correlated with environmental factors such as precipitation, elevation and latitude. Based on EST-SSR markers, we detected a low level of genetic diversity in natural A. venetum populations from the coastal habitats (H e = 0.342) ( Table 4). Although EST-SSRs are expected to be more conserved and exhibit a lower rate of polymorphism than genomic SSRs [40][41][42][43], there are also studies showing that there is no significant differences between these two types of markers [44][45][46]. Considering the small population size of the A. venetum populations in the coastal region, future studies that combine the analysis of genomic SSR markers or genome-wide SNPs could better represent the real level of genetic diversity in A. venetum from this region.
Coastal plant species are expected to be highly vulnerable to habitat alteration, due to their high specialization to coastal environments [6]. Through nuclear SSR analysis, a low level of genetic diversity has been observed in coastal herb, shrub and tree species [47][48][49]. Peng et al. [50] investigated the genetic diversity of nine wild A. venetum populations from eight provinces of China using RAPD markers, and their study showed that populations from coastal regions such as Jiangsu and Jilin exhibit lower genetic diversity than those from the arid northwest regions such as Xinjiang, Inner Mongolia and Ningxia. In a comparative study of Fraxinus angustifolia Vahl populations from the Continental region and the Mediterranean region, Temunović et al. [51] detected a significantly lower genetic diversity and higher population divergence in the latter region, revealing the influence of environmental heterogeneity in shaping the genetic variation between divergent habitats. Therefore, coastal A. venetum populations might represent a divergent ecotype, and environmental factors may cause differences in genetic diversity between arid and subhumid regions. By including environmental data in our future analyses, we could further test whether the differences are correlated with diverse natural conditions (subhumid vs. arid areas). Another plausible explanation for this low level of genetic diversity is that A. venetum populations in the coastal region have experienced severe reductions in population size, leading to the loss of genetic diversity and increased susceptibility. Using Bottleneck software, we found that two populations (QG and BH) showed distortions (modeshifts) in their allele frequency distributions, suggesting recent bottleneck events, and Wilcoxon's test also revealed significant heterozygote excess in the BH population (P < 0.05; Table 6). According to our field observations, for populations such as BH and SY, situated along the highway road, the natural habitat is severely being disturbed. In addition, the population size of most populations is quite small (N < 30). Based on scaled effective population sizes Θ (4N e μ), the estimated N e of the A. venetum populations in the Yancheng region is approximately 10-11 individuals. According to a genetic model predicting the proportion of initial heterozygosity retained per generation [1-(1/2N e )] [52], A. venetum populations would be expected to lose 4.5% of their heterozygosity per generation. Therefore, with declining genetic diversity and continuing demographic changes over time, the A. venetum populations in the Yancheng coastal region are of potential concern in terms of decreases in individual fitness and population viability, and an increased risk of extinction in the future.

Moderate level of gene connection with a decreasing trend in A. venetum coastal populations
Gene flow is crucial for population resilience and persistence following habitat disturbances and environmental changes. This process may provide a source of recruitment and maintain genetic diversity by introducing adaptive alleles from other populations [53][54][55]. In this study, we found little evidence of strong genetic structuring in the A. venetum populations across the Yancheng coastal region. Genetic structure analysis grouped all of the sampled populations into two clusters, and in the majority of populations, a number of admixed plants existed (Fig. 3c). Negative multilocus F IS values (mean value of − 0.08) were found for A. venetum, indicating an absence of inbreeding in the coastal populations of this species (Table 4). In accordance with this, contemporary migration rates among A. venetum coastal populations were found to be moderate (mean m c = 0.056) (Additional file 3: Table S1). Together, these results indicate that population connectivity in A. venetum coastal populations has not yet been greatly disturbed yet. In general, life-history traits of plant species such as their pollen and seed dispersal mode, often affect the genetic connection between populations [55,56]. The study of pollination biology has shown that Apis mellifera and Ichneumon sp. are the main pollinators of A. venetum [57]. Due to the limited foraging distances of these pollinators, gene dispersal through pollen alone seems insufficient to maintain the species' population connectivity across the coastal region. The seeds of A. venetum are less than 1 mm in size and are covered with pappus-like hair, which enables an anemochorous dispersion [58,59]. The high seed dispersal capabilities of A. venetum may play a key role in maintaining the moderate levels of ongoing gene flow and population connectivity across the coastal region. Although studies have shown that plants relying on wind for pollination or seed dispersal may be subject to less negative effects of habitat alteration on their genetic diversity or population connection [55,56], we still detected a decrease in gene flow over the last few generations compared to the historical level of migration rates (mean m h = 0.106). Under the persistent pressures of anthropogenic and climate disturbances, whether natural A. venetum populations, especially those with low genetic diversity and a small population size, have a sufficient capacity to survive following disturbance or adapt to future environmental change still requires long-term monitoring.

Conclusion
In this study, we acquired a high-quality transcriptome from A. venetum leaves by using the Illumina sequencing platform and successfully developed twelve polymorphic EST-SSR markers in A. venetum. Using these informative makers, we detected a low level of genetic diversity and bottleneck events in A. venetum populations across their natural distribution in the Yancheng coastal region of Jiangsu Province. Although population connectivity between A. venetum populations has been maintained by the high seed dispersal ability of the species, considering the ongoing anthropogenic activities, long-term monitoring and conservation strategies should be implemented to better protect these small populations. The conservation of available genetic diversity is essential to enable the continued utilization of this economically important plant. Based on the results of this study, we suggested that populations that have undergone bottlenecks (e.g., population BH and QG) and those with a high level of genetic diversity (population XY), should be given conservation priority. In addition, there are not enough comparative data available on the extent of genetic diversity in A. venetum across its global biogeographical ranges, and further studies combining samples from arid regions and more markers should be conducted to gain a more precise understanding of the genetic diversity, population structure, and evolutionary history of this important plant species.

Sample collection and DNA and RNA extraction
A total of 103 leaf samples of A. venetum were collected from 10 to 24 individuals in six wild populations, representing most of the distribution area of this species in Jiangsu Province (Fig. 4). Detailed information for all the populations is listed in Table 4. A. venetum is mainly distributed along the riversides or roadsides, and permission was not necessary to collect these samples. The collection of plant material complied with national guidelines. All sampled leaves were immediately dried in silica gel. The materials were identified by Jiangsu Province and Chinese Academy of Sciences, China. We used QIAGEN Plant DNA kit (Gaithersburg, MD) to extract  Table 4 the total genomic DNA according to the manufacturer's protocol. Besides, fresh leaves of three individuals from the LD population were immediately frozen in liquid nitrogen and stored at − 70°C for total RNA extraction. Total RNA from these individuals was extracted using the QIAGEN RNeasy Plant Mini Kit following the manufacturer's procedures (Chatsworth, CA). Then the quantified total RNAs were sent to Novogene Bioinformatics Institute (Beijing, China) for further processing.

cDNA library construction and transcriptome sequencing
First, the RNA concentration and integrity of the samples were measured using the Qubit RNA Assay Kit in a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA) and the RNA Nano 6000 Assay Kit with the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). Then, for each sample, cDNA library were construction and sequenced by using equal amounts of qualified RNA according to the standard procedures of Novogene Bioinformatics Institute (Beijing, China). The NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) were used to generate sequencing libraries.
After purification and quality assessment, the three library preparations were sequenced on the Illumina HiSeq platform (Illumina, USA).

Transcriptome assembly and functional annotation
Raw RNA-seq reads were processed using in-house Perl scripts to remove reads containing adapters and ploy-N sequences (greater than 5%) and reads with more than 20% low-quality bases (quality scores < 10). Then, we used Trinity software [60] with the default parameters to assemble the high-quality clean data. All unigenes were queried against the NCBI Nr (non-redundant protein) and Nt (non-redundant nucleotide) databases; the Swis-sProt protein database (http://www.expasy.ch/sprot); Pfam (Protein family database); the KOG (Clusters of Orthologous Groups of proteins) database and KO (KEGG Orthology) database. Gene ontology (GO) annotation [61] of the unigenes was performed using BLAS-T2go [62]. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (https://www.genome.jp/ kegg/) were determined with an E-value cut-off of 1e-5. The WEGO [63] were used to plot the distributions of level-2 GO terms with functional classification.

Identification of EST-SSRs
The MicroSAtellite (MISA) program [64] was used to detect transcripts containing EST-SSRs. A minimum repeat number of six for dinucleotide motifs and five for tri-, tetra-, penta-, and hexanucleotide motifs were set as detection criteria. The sequences containing SSRs were submitted to Primer premier 5.0 software (Premier Biosoft International, Palo Alto, CA, USA) to design primers. The parameter settings were as follows: product size ranging from 100 to 300 bp; primer length ranging from 18 to 25 bp; GC content between 40 and 60% and the annealing temperature between 55 and 65°C [65].
PCR amplification was carried out in a 15 μL total reaction volume containing 20-40 ng genomic DNA, 7.5 μl of 2 × Taq PCR MasterMix (Tiangen, Beijing, China) and 0.3 μM of each primer. The PCR procedure included an initial denaturation for 3 min at 94°C, followed by 35 cycles of 30 s at 94°C, 30 s at optimal annealing temperature for each locus, and 15 s at 72°C, followed by a final extension of 5 min at 72°C. The PCR products were checked using silver-stained nondenaturing polyacrylamide gels. Then the optimized SSR primers were further labelled with 6-FAM or HEX fluorescein dye (Sangon Biotech, Shanghai, China). After PCR amplification, allele identification and genotyping were performed with GeneMarker version 2.2.0 (SoftGenetics, State College, Pennsylvania, USA).

Genetic diversity analyses
The number of alleles (N A ), polymorphism information content (PIC), expected heterozygosity (H E ) and observed heterozygosity (H O ) of each EST-SSR locus were estimated with CERVUS 2.0 [66]. Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD) tests were performed using GENEPOP version 4.2 with a Bonferroni correction [67]. Null allele frequencies, stuttering, and large allele dropout were detected using the MICRO-CHECKER version 2.2.3 [68] program. For each population of A. venetum, the number of alleles (N A ), allelic richness (AR), expected heterozygosity (H E ), observed heterozygosity (H O ), inbreeding coefficient (F IS ) and differentiation among populations (F ST ) were calculated using FSTAT version 2.9.3.2 [69].
Population genetic structure, gene flow and demographic analyses Population structure was investigated by using the SRUC-TURE 2.3.4 program [70] implementing a model based Bayesian approach. The value of genetic clusters (K) was set from 1 to 6, assuming an admixture model and independent allele frequencies. Ten independent runs were conducted for each K with a burn-in of 10,000 and 100, 000 Markov Chain Monte Carlo replicates. The most possible K value was chosen by calculating ΔK [71] in STRUCTURE HARVESTER [72].
Contemporary inter-population migration between A. venetum populations was estimated using BAYESASS version 1.3 [73]. The delta values for allele frequencies, migration rates, and inbreeding coefficients were adjusted accordingly to ensure that the acceptance rates fell between 40 and 60% [73]. We performed the software for 10 7 iterations with a burn-in of 10 6 generations. Ten replicate runs were conducted with a different initial seed. We assessed the model convergence by comparing the posterior probability densities of parameter estimates across these ten runs. The results presented were from the best-fit run. We also estimated historical gene flow using a coalescent-based Bayesian method implemented in the program MIGRATE-N version 3.6 [74]. The process used the Brownian motion approximation as the mutation model following 10 7 iterations with a burn-in of 10 5 . The static heating scheme was used with four chains at different temperatures (1.0, 1.5, 3.0, and 100, 000.0) [75]. Two parameters, scaled effective population sizes Θ (4N e μ, where N e is effective population size, μ is the average mutation rate of microsatellites as 10 − 3 per generation) and scaled immigration rates M (m h /μ, where m h is historical migration rate) between pairs of populations over around 4N e generations were estimated simultaneously. The number of immigrants per generation Nm was estimated by the equation Nm = ΘM/4. A Wilcoxon signed-rank test was conducted to compare historical and contemporary gene flow estimates among the A. venetum populations.
To test for isolation-by-distance (IBD), the correlation of pairwise geographical distance (log geographical distance in km) and genetic distance (F ST /1-F ST ) values was evaluated. Statistical significance was tested with 1000 permutations of the Mantel test via the R package Vegan [76]. Wilcoxon's signed rank test and the mode-shift test in BOTTLENECK version 1.2.02 [77] were used to determine whether the A. venetum populations in the Yancheng coastal region have undergone significant reductions in the effective population size (N e ). The first methodology compares the heterozygosity expected (H E ) at Hardy-Weinberg equilibrium with the heterozygosity expected at mutation-drift equilibrium (H eq ) [78], which is suitable for detecting bottlenecks occurring in the last 2-4N e generations. The second methodology is based on the allele frequency distribution, which is more appropriate for detecting population declines that have occurred more recently (approximately the last few dozen generations) [79,80]. For population that has not experienced bottleneck event, a large proportion of alleles at a low frequency and a smaller proportion of alleles at intermediate frequencies distribution (L-shape distribution) presents. While in bottlenecked populations, a shifted mode of allele frequency distribution would be detected. We performed 10,000 simulations under both the stepwise mutation model (SMM) and the two-phase model (TPM) with 95% single-step mutations and 5% multistep mutations for each A. venetum population. P-values were assessed for statistical significance at the 0.05 level.