Genetic Diversity and Population Genetics of Cuttlefish Sepiella japonica based on Newly Developed SSR Markers

Cuttlefish Sepiella japonica , an economically valuable fishery resource that is considered as one of the four most famous sea products in China, has dramatically declined since the 1990s to near extinction, presumably due to over exploitation and environmental disruption and particularly because of the deterioration of their spawning grounds. In this study, we developed primers for nine microsatellite loci resulting from an enrichment process based on the Wenzhou population. Eight out of nine microsatellite loci were used to analyze population genetic diversity of cuttlefish in China Sea. The results showed high levels of genetic diversity in all populations with highest mean number of alleles and mean effective alleles estimated in population Ningde and lowest in Qingdao. Heterozygosity at species level was high and F st ranged between 0.105 and 0.271, overall, without signs of inbreeding except for Qingdao population, with F st = 0.271. The studied populations did not show a clear genetic structure with an average differentiation estimated between populations of 0.014. Knowledge obtained from this study has important implications for the management of existing genetic resources for aquaculture, artificial releasing project and conservation in cuttlefish .


INTRODUCTION
S epiella japonica (Decapoda, Sepiidae), is an economically important cephalopod species that inhabits all over the Chinese coastal waters from Yellow Sea to South China Sea. The biological characteristics, life cycle (not longer than half a year), high natural mortality and year to year increasing abundance variability, which make this marine species vulnerable to anthropogenic pressures. The wild populations of this species dramatically declined since the 1990s and the species has become nearly extinct, presumably due to over exploitation, environmental disruption, and in particular, the deterioration of its spawning grounds (Wu et al., 2006). Aiming to protect and recover this valuable cuttlefish wild resource, our research team has made many efforts to enhance the stocks by releasing cuttlefish eggs and larvae to native waters since 2006. Genetic concerns related to stock enhancement include a potential disturbance of spatial population structure along with negative effects on fitness and genetic diversity of wild stocks (Kitada et al., 2009;Laikre et al., 2010), thus irresponsible stock enhancement can bring harmful genetic effects on wild stocks (Blankenship and Leber, 1995;Lorenzen et al., 2010). Knowledge of the genetic structure of supplemented populations is essential to guide stock enhancement actions and to design appropriate spatial management plans.
So far, there were only few recorded population genetic studies about cuttlefish, using allozyme electrophoresis and the COI gene (Zheng et al., 2001) or 16S rRNA gene (Zheng et al., 2005) analysis. These studies showed low genetic diversity in cuttlefish, and no significant genetic differentiation was detected among different geographic populations in China Sea areas. Similar results were also obtained by COI gene data analysis (Li et al., 2014).

O n l i n e F i r s t A r t i c l e
Considering that these results could be due to intrinsically characteristics of the used markers, it is essential to use higher resolution molecular markers to further evaluate the genetic diversity and genetic variation of cuttlefish in China.
Microsatellites or SSRs are defined as highly variable DNA sequences composed of tandem repeats of 1-6 nucleotides with codominant inheritance and have become the markers of choice for a variety of applications. Their abundant presence in genomes, high mutation rates and multi-allelic nature (Ellegren, 2000) make such markers among the best choices for analysis of continental and regional level questions of population genetics (Rosenberg et al., 2002;Sahoo and Kashyap, 2005). During the past several years, the development of SSRs technology from PAGE based to sequencing based, which improved the resolution of the SSRs technology. Although a few microsatellite markers have been developed for cuttlefish (Wu et al., 2010a;Guo et al., 2013), their polymorphism is limited.
The main objectives of the present study were to develop new polymorphic microsatellite markers for cuttlefish; provide information about the genetic diversity and genetic structure of cuttlefish at regional scale; and guide the selection of parents of cuttlefish for subsequent stock enhancement to support the management and conservation of cuttlefish.

Sampling and study area
A total of 179 cuttlefish individuals from 5 populations were collected from September 2011 to November 2012 across the Chinese coastal areas throughout the Yellow Sea (Qingdao, QD, 120°45´N, 36°18´E), East China Sea (Zhoushan, ZS, 122°48´N, 30°30´E; Wenzhou, WZ, 121°35´N, 28°25´E; Ningde, ND, 120°51´N, 27°45´E) and South China Sea (Zhanjiang, ZJ, 109°50´N, 21°25´E) ( Fig.  1). No specific permits were required for the described field studies. Samples for this study were contributed by several collectors with the help of a trap net. A small part (about 0.5 g) of the cuttlefish wrist muscle was removed by scissor and stored at -80°C until it was used for DNA extraction, and all cuttlefish were released afterwards. Total genomic DNA was extracted from mantle muscle tissue of periopods in all the 179 samples using the DNeasyTM Tissue Kit (Qiagen).

Development of microsatellite markers
The isolation and characterization of the microsatellite loci were conducted with slight modifications for enriched genomic libraries, as described by Guo et al. (2013).
Briefly, total genomic DNA, obtained from 39 individuals from population Wenzhou, was digested with the restriction enzyme Tru II (Fermentas) and the resulting fragments, ranging from 400 to 1000 bp, were ligated to the specific adaptors, namely, Tru II A 5'-TACTCAGGACTCAT-3' and Tru II B 5'-ACGATGAGTCCTGAG-3'. The DNA ligated to the adaptor was amplified using the Tru IIA adaptor as the primer. The adaptor-ligated DNA fragments were selected by hybridization to a biotin-(CA) 12 probe and captured with streptavidin-conjugated magnetic beads (Dynabeads, Dynal Biotech). After two rounds of PCR amplification, the adaptors were cleaved internally with MboI (Fermentas), following, the microsatellite-enriched library was cloned into a PMD18-T vector (TaKaRa) and transformed into Escherichia coli DH5a competent cells. After the screening for the SSRs, the positive clones were sequenced and chosen for the primer designation.
The loci with optimal polymorphism were selected for further investigation and, in these cases, forward primers were labeled at the 5' end with a fluorescent dye (HEX or 6-FAM) to generate fluorescently labeled microsatellite fragments.

O n l i n e F i r s t A r t i c l e Preliminary microsatellite selection and full-scale genotyping
Nine new developed microsatellite loci as listed in Table I were chosen for a preliminary analysis of the genetic variation in the species using population Wenzhou. Microsatellite loci PCR amplifications were conducted in 15μl reactions with 10 pmol of primers each (forward primer in each pair was 5' labeled with a fluorescent dye), 10-100 ng of genomic DNA, 1×PCR buffer (TaKaRa), 1.5 mM MgCl 2 , 0.2 mM of each dNTP (TaKaRa) and 0.5 U Taq DNA polymerase (TaKaRa) on a MJ Mini Thermal-Cycler (Biorad). The thermal cycling parameters were: initial denaturation at 94°C for 5 min, followed by 35 cycles of denaturation at 94°C for 1 min, primer-specific annealing temperature (Table I) for 45 s, extension at 72°C for 1 min and final extension at 72°C for 7 min. Following PCR amplification, the extension products were resolved on 1.5% agarose gel. The final concentration of agarose gel may in range of 2~3% after microwave heating. Genotyping was performed on an ABI PRISM 310 Genetic analyzer (Applied Biosystems), and the alleles were sized relative to an internal size standard (ROX GS 400HD; Applied Biosystems) using GENESCAN 3.7 (Applied Biosystems). Full-scale genotyping of all samples and populations followed the same conditions as described for the preliminary analysis.

Statistical analysis
The raw data were initially analyzed using the Microchecker 2.2.3 software (Oosterhout et al., 2004) to check the microsatellites for null alleles and scoring errors. The characteristics of the microsatellite loci, including the number of alleles, observed (H o ) and expected (H e ) heterozygosity, deviations from the Hardy-Weinberg equilibrium (HWE) and linkage equilibrium, were determined using GenAlEx 6.5 (Peakall et al., 2012). The polymorphic information content (PIC) was calculated by the formula in Botstein et al. (1980). A pairwise F st was calculated to assess the genetic differentiation among the populations in Arlequin 3.5 (Excoffier and Lischer, 2010); comparing the pairwise F st values consisted of randomizing multi-locus genotypes between each pair of samples (10,000 permutations) to check the significance. F st and the genetic structure of the populations were further quantified using a two hierarchical analysis of molecular variance (AMOVA) (Michalakis and Excoffier, 1996) with Arlequin version 3.5 (Excoffier and Lischer, 2010). The neighbour-joining algorithm (Saitou and Nei, 1987) tree was constructed by NEIGHBOR and CONSENSE within the populations 1.2.32 software package (Langella, 1999) based on the Cavalli-Sforza's and Edwards genetic distance (D C , 1000 bootstrap replications), which was computed by Microsatellite Analyser software (Dieringer and Schlötterer, 2003). We used individual assignment tests implemented in STRUCTURE 2.3.4 (Pritchard et al., 2000) to estimate the number of genetically discrete populations. Population groups were simulated for K=1 to 8, assuming possible mixed ancestry and correlated allele frequencies, and using 1,000,000 MCMC steps, discarding the first 100,000 steps. The best K values were analyzed using the Evanno et al. (2005) method in STRUCTURE Harvester (Earl and von Holdt, 2012). Evidence of recent bottlenecks in each population was examined using BOTTLENECK version 1.2.02 (Piry et al., 1999). The probability distribution was estimated using Wilcoxon's signed-rank test with 10,000 simulations. The infinite allele model (IAM), stepwise mutation model (SMM), and two-phase model (TPM) of microsatellite mutations were implemented. For the TPM, we used a model of 95% single-step mutations, and for the remaining 5%, multistep mutations with a variance of 10. The two-phase model (TPM) of microsatellite mutations with 95% single-step mutations and 5% multiple-step mutations, and a variance among multiple steps of mutations were implemented (Wilcoxon's test two tails for heterozygosity excess and deficiency with 10,000 replications).

Development of microsatellite markers
Some difficulties were experienced when developing the new microsatellite loci, such as short sequences, the absence of flanking regions in one or both ends, PCR instability, stutter and difficulties in genotyping contributed to a reduced number of employed markers. However, nine selected microsatellite loci tested were highly polymorphic and amplified reliably, in a preliminary analysis using only the Wenzhou population (Table I). No evidence of linkage disequilibrium was observed, and an initial test for Hardy-Weinberg equilibrium (HWE) showed significant deviation at locus S.m17, S.m26 and S.m35 (Table I). The Microchecker analysis detected scoring errors and null alleles at locus S.m26, which was excluded in the full genotyping analyses.

Genetic diversity analyses
The final eight microsatellite loci examined were polymorphic in all five populations. We observed a total of 234 alleles across eight loci in five different geographic regions. The total number of alleles per locus ranged  (Table II)  For the bottleneck tests, under the infinite allele model (IAM), the two-phase model (TPM) and stepwise mutation model (SMM), no significant heterozygote excess nor deficiency was found in any of the populations (Table III). Table II  Genetic structure analyses STRUCTURE Harvester indicated K= 2 as the best value (Fig. 2) and several populations showed a high level of admixture, such as ND, ZS and QD (Fig. 3).

Table III. detection of bottlenecks based on Wilcoxon's test for three mutation models of five S. japonica populations in Chinese coastal waters (see
The neighbor-joining tree generated from D C values (Fig. 4) shows a polytomy between populations WZ, QD and a clade formed by ND, ZS and ZJ populations. The ND and ZS populations clustered together, in a group sister to the ZJ population.
F st values ranged between 0.003 and 0.024 (Table   IV) and overall, pairwise population differentiation was 0.0142. A hierarchical AMOVA revealed that most of the genetic variance was found within populations (98.48%, P < 0.001) ( Table V).

O n l i n e F i r s t A r t i c l e
Y. Ye et al. Fig. 4. Neighbour-joining clustering of samples based on D C , percentage at node was a bootstrap value obtained from 1000 replicates and population codes are as in Fig. 1.

Genetic diversity of cuttlefish populations
Genetic diversity is influenced by many factors including historical factors, anthropogenic activity, habitation and a low rate of mitochondrial evolution (Avise, 1994). Historical factors may play an important role in determining the patterns of genetic variability (Yamaguchi and Tsukaya, 2010). Thus, a high level of genetic diversity is essential for long-term survival of populations (Nie et al., 2015). Microsatellites have been a useful molecular marker for revealing both genetic variability and genetic structuring for populations in different species within the cephalopodan class, such as the cuttlefish Sepia officinalis (Perez-Losada et al., 2002) and the European squid Loligo vulgaris (Shaw et al., 1999;Garoia et al., 2004). In this study, we developed nine new highly polymorphic microsatellite loci for cuttlefish from CA dinucleotide enriched genomic libraries. Despite many positive clones derived from the enriched genomic libraries some difficulties were experienced, similarly to what was reported by other researchers working with marine invertebrate species, mainly because microsatellite repeats in invertebrates are typically less abundant and shorter than in vertebrates (Chambers and MacAvoy, 2000). The microsatellites generally exhibited high allelic diversity and significant deviations from the Hardy-Weinberg equilibrium recorded for three loci could be connected to the high heterozygosity values estimated. The allele method of detection seem to influence the number of alleles since, for O. vulgaris (Cabranes et al., 2008;Moreira et al., 2011), the number of alleles detected with primers labeled with fluorescent dye was much higher than those detected using a silver staining method, maybe due to microsatellite allele sizing methodology.
Higher levels of genetic diversity were displayed in this study, for all five populations of cuttlefish, when compared with previous studies for the same species, using samples collected before population enhancements began. These studies used molecular techniques sequence divergence of the COI gene (Zheng, 2001) and of the 16S rRNA gene (Zheng et al., 2005). However, in the study of Li et al. (2014), conducted after stocking enhancement, the COI gene data also showed a low genetic diversity for cuttlefish. This is not an unexpected result given the proposed higher rates of mutation and improved resolution of microsatellite loci relative to allozyme loci and agrees well with the levels of microsatellite variability reported for squid Loligo forbesi (Shaw et al., 1999) and Ommestrephidae (Adcock et al., 1999). Since a population genetic analysis of cuttlefish with microsatellite data before the stocking enhancement was not performed, a final conclusion as to whether nearly ten years of stocking enhancement has affected the wild population and at what extent cannot be fully determined. However, the microsatellite data in our research demonstrated that the genetic variance of cuttlefish is optimistic, and the releasing of eggs and larvae has aided the recovering cuttlefish population and the practice of selecting the parent cuttlefish is appropriate. For instance, the ND population showed the highest level of genetic diversity and is the preferred selection for parent cuttlefish to recover the valuable resource; in addition, in fact, the release of eggs and larvae of ZS was the offspring of the ND population in this study.

Genetic structure of cuttlefish populations
Pelagic organisms in the open ocean generally do not show high levels of population differentiation, most likely resulting from many opportunities of dispersal in various life cycle stages and the lack of physical barriers in the environment. However, since cuttlefish is an annual species, naturally occurring gene flow is expected to be reduced (Zhang, 2007), probably due to the very short pelagic period and the inability of long-distance migration O n l i n e

F i r s t A r t i c l e
Genetic Diversity of Sepiella japonica in China (Wu et al., 2010b). In our study, although no clear population genetic structure appears to exist in the target area, with only two genetic groups estimated, the high levels of admixture occurring in three of the screened populations, seems to contradict reduced gene flow. Notwithstanding, before the decline of this natural resource, Ni and Xu (1986) indicated that cuttlefish populations in the East China Sea divided into two groups: the Zhongjieshan-Yushan group (equivalent to Zhoushan-Wenzhou populations) and the Mindong-Beiji group (equivalent to Ningde population) based on morphological characteristics. Furthermore, after the sharp decline in the resource, intensive stock enhancement activities were practiced with the release of larvae and fertilized eggs. The detection of two original population groups is thus in accordance to our genetic structure analysis (K= 2) and previous genetic studies (Li et al., 2014), while the present level of admixture is most likely a result of human-mediated translocation of eggs and larvae during the population enhancement actions.

Implications of the genetic differentiation of cuttlefish for sustainable use of fishery resources
The level of genetic diversity found among the cuttlefish populations in our study is an important variable to be considered in the sustainable management of the cuttlefish fishery industry because the decline of a given genetic population due to overfishing may limit its recovery. Further steps have to be implemented in identifying new breeding areas and consequential genetic stock boundaries to establish management units that will be pivotal to the long-term sustainable use of this valuable marine resource. This manuscript suggests that the source of the cuttlefish released to the sea should come from populations having higher genetic diversity (for example, the ND population) to avoid low genetic diversity and a possible repeat of a sharp population decline.

CONCLUSIONS
The current Sepiella japonica resource exhibited a high degree of genetic diversity, which revealed that our stock enhancement has not resulted in evident negative effects. However, due to the special biological character of an annual species, which requires particular attention due to species vulnerability, continuous genetic monitoring of the populations should be conducted to maintain a longterm sustainable use of this valuable marine resource.

O n l i n e F i r s t A r t i c l e O n l i n e F i r s t A r t i c l e
Genetic Diversity of Sepiella japonica in China