1 INTRODUCTION

Microsatellites, also known as simple sequence repeats (SSRs), are randomly repeated DNA sequences composed of 1-6 bp long units (Tautz and Renz, 1984). They are ubiquitous in prokaryote and eukaryote genomes, present in both protein encoding and noncoding regions (Tóth et al., 2000 ; Subramanian et al., 2003 ; Hoffman and Nichols, 2011). SSRs can generate and maintain extensive length polymorphism because they have a high enough mutation rate of between 10 - 3 and 10 - 4 mutations per gamete per generation (Weber and Wong, 1993 ; Bhargava and Fuentes, 2010), which makes them important for genome evolution (Tautz et al., 1986 ; Kashi et al., 1997). Additionally, because of their locus-specificity, co-dominant inheritance, reproducibility, and abundance SSR markers have become powerful molecular tools for a wide range of applications in genetic analysis and breeding, such as population genetics, parentage analysis, genetic mapping, association mapping, comparative mapping, and quantitative trait loci (QTL) analysis (Bruford and Wayne, 1993 ; Jarne and Lagodav, 1996 ; Cruz et al., 2005 ; Varshney et al., 2005 ; Bohra et al., 2011). Although the advent of SNP markers has attracted more attention and funding for their development and applications in various animal and plant genera, many studies have indicated that SSR markers are still being extensively developed as an irreplaceable molecular tool for various purposes (Barchi et al., 2011 ; Ansari et al., 2013 ; Barzegar et al., 2013 ; Doulati-Baneh et al., 2013 ; Guo et al., 2013 ; Liu et al., 2013 ; Smee et al., 2013 ; Yuan et al., 2013 ; Zitouna et al., 2013 ; Vukosavljev et al., 2013 ; Lee et al., 2014 ; Parobek et al., 2014).

However, conventional experimental approaches to develop SSRs for non-model species are laborious, time consuming, and expensive (Zane et al., 2002). Fortunately, recent developments in DNA sequencing technology and the falling costs of next-generation sequencing (NGS) have provided novel methods for SSR identification with high efficiency and lower costs (Castoe et al., 2010 ; Jennings et al., 2011 ; Triwitayakorn et al., 2011 ; Wang et al., 2011, 2013 ; Nybom et al., 2014 ; Tadano et al., 2014 ; Králová-Hromadová et al., 2015). Genome survey sequencing (GSS) based on the NGS platform, proven particularly useful in identifying genome-wide SSRs in non-model species with the advantages of low costs and high throughput (Miller et al., 2007 ; Barchi et al., 2011 ; Rowe et al., 2011 ; Zhou et al., 2013 ; Xu et al., 2014). Previous studies have demonstrated that SSR markers distributed in noncoding regions have practical implications (Hancock, 1995 ; Hosseini et al., 2008). Although genic-SSRs may have an insufficient degree of polymorphism for genetic analysis because they are often more conserved, SSR markers developed in functional genes can also be useful for association analysis and marker-aided selection (Cheng et al., 2007 ; Shikano et al., 2010 ; Guichoux et al., 2011). The GSS approach could solve the above issue, because it is not only productive and efficient for identifying SSRs in non-coding regions, but it also efficiently predicts putative gene functions and identifies potential exon-intron boundaries because it possesses high similarity and long sequence lengths (Strong and Nelson, 2000).

The kuruma prawn, Marsupenaeus japonicus, is one of the most cultivated and consumed shrimps with a native range encompassing an Indo-West Pacific distribution from Japan, through Southeast Asia to East Africa and the Red Sea (Holthuis, 1980 ; Hewitt and Duncan, 2001). However, very few molecular genetic/genomic resources are publically available for M. japonicus. The characterization and distribution of SSRs remains ambiguous for this species, and the use of SSR markers in genomic studies and marker-assisted selection is limited. Consequently, the aim of the present study was to characterize and develop genome-wide SSR markers in M. Japonicus by GSS. The newly identified SSRs would be useful for extending our current knowledge of M. japonicus genome organization and for genome mapping, genome-wide association studies, markeraided selection, comparative genomic analysis, and population genetics.

2 MATERIAL AND METHOD 2.1 Genome survey sequencing

An M. japonicus individual, the maternal parent of an F 1 mapping population, was chosen for genome survey sequencing. Genomic DNA was extracted from its muscles following the methods described by Wang et al.(2006). Two paired-end DNA libraries were constructed with insert sizes of 180 and 500 bp, and then sequenced using the Illumina HiSeq2000 platform following the manufacturer's protocol, which has been described by Etter et al.(2011).

Because Illumina HiSeq2000 produces some lowquality reads, a quality trim was performed on the first raw reads to make the further analysis more accurate, as follows:(1) removal of adapter contamination;(2) removal of bases that were not A, G, C, and T from the beginning of 5′ end;(3) trimming the ends of the lowquality reads (sequencing quality value ＜20);(4) removal of reads with more than 10% “N” bases;(5) abandon reads ＜100 bp following adapter removal and quality trimming. The clean reads were assembled into contigs in SOAPdenove (v2.04) with a K-mer of 17 by applying the de Bruijn graph structure (Jiao et al., 2014). The paired-end information was then used to join the unique contigs into scaffolds. After that, paired-end extracted reads that had one read uniquely aligned on a contig and another read located in a gap region were used to fill the intra-scaffold gaps (Jiao et al., 2014). The details for the SOAPdenove assembly are described in Li et al.(2010). Repetitive element occurrence was detected with CENSOR, which screens query sequences against a reference collection of repeats (http://www.girinst.org/censor)(Kohany et al., 2006 ; Bohra et al., 2014), adopting default parameters and using Viridiplantae as the target database.

2.2 Microsatellite mining and primer design

To identify SSR markers, the assembled scaffolds were screened for SSRs in SciRoKo (Kofler et al., 2007). To assess the distribution and characterization of microsatellite abundance, the mono-, di-, tri-, tetra-, penta-, and hexanucleotide motifs were screened from the scaffolds. In the search for an SSR standard, we defined SSRs as mononucleotides ≥10 bases, dinucleotides ≥12 bases, trinucleotides ≥15 bases, tetranucleotides ≥20 bases, pentanucleotides ≥25 bases, and hexanucleotides ≥30 bases (Cardle et al., 2000); the maximum number of bases interrupting two SSRs in a compound microsatellite was 100 bases. Additionally, because DNA is double stranded and an SSR start site could be arbitrarily chosen, reverse-complement repeat motifs and translated or shift motifs were also considered and grouped together in the analysis (e.g., AG/CT representing AG and GA, TC and CT)(Jurka and Pethiyagod, 1995). The software Websat (http://wsmartins.net/websat/) was used to design SSR primers (Martins et al., 2009). Primer pairs were designed to meet the following restrictions: the target amplification size range was 100-500 bp, the primer annealing temperature was restricted to 55-60℃ and the primer length was 19- 27 bp. We termed such loci with PCR primers as potential microsatellite markers (Potentially Amplifiable Loci or PAL).

2.3 Functional annotation of genome-wide SSR markers

For homology annotation of all of the SSRcontaining sequences, non-redundant sequences were subjected to the public database GO (Gene Ontology, http://www.geneontology.org/). The sequence annotation was based on sequence similarity according to gene ontology (Barchi et al., 2011) with the annotation parameters at an E-value cut-offof 10^-10.

2.4 Laboratory verification of the SSR markers in a natural population

We randomly selected 30 primer pairs from the microsatellite loci with relatively longer repeats for laboratory verification, containing six di-, tri-, tetra-, penta-, and hexanucleotide loci. Forty-two M. japonicus individuals collected from a natural population on the coast of Fujian Province (China) were used for verification of the SSR primers. The total genomic DNA was also extracted from muscle tissue following the methods described by Wang et al.(2006). The PCR reactions were conducted in a 10 μL PCR mix, containing 1×PCR buffer, 1.5 mmol/L gCl₂, 0.2 mmol/L dNTPs, 0.2 U Taq DNA polymerase (Promega, USA), approximately 20 ng of genomic DNA, and 0.25 μmol/L of each primer pair. PCR amplification was performed in a MJ-PCR-200 thermal cycler (Bio-Rad, USA) with the following procedure: initial denaturation at 94℃ for 5 min, 35 cycles of denaturation at 94℃ for 30 s, annealing temperature for 30 s (Table 1), and extension at 72℃ for 40 s, followed by a final step at 72℃ for 10 min. The PCR products were evaluated by 20 g/L agarose gel electrophoresis and then loaded into an Applied Biosystems 3130 sequencer and scored using GeneMapper version 3.7 software (Applied Biosystems) with GeneScan-500 LIZ Size Standard as an internal size standard. For each polymorphic SSR locus amplified, polymorphic information content (PIC) was calculated with the following formula:

where, p _i s the proportion of the i th allele (Biswas et al., 2014).

3 RESULT 3.1 Genome survey sequencing

Genome survey sequencing was performed for a single M. japonicus individual by sequencing two paired-end DNA libraries with insert sizes of 180 and 500 bp based on the Illumina HiSeq2000 platform. All of the raw data are accessible in the Short Read Archive (SRA) of the National Center for Biotechnology Information (NCBI) under the accession number SRX1162702. According to the Lander-Waterman algorithm, the M. japonicus genome size was calculated by the following equation:

where, k _num is the number of K-mer, k _depth is the expected depth of K-mer, b _num is the number of bases, and b _depth is the expected depth of bases. According to the equation, the size M. japonicus genome was estimated as ~2.4 Gb. In total, 79.75 Gb of sequencing data were produced, equivalent to ~33.23 × genome coverage (Table 2). The assembled data with a K-mer of 29 was selected for the final results according to the comprehensive evaluation of technique indexes (such as the length of all the scaffolds, number of scaffolds, and scaffold N50). The assembly statistics of the genome survey sequencing data is summarized in Table 3. A total of 3 244 404 scaffolds were assembled from the genome survey sequencing data, among which there were 127 538 ones larger than 1 000 bp with the largest length of 15 860 bp.

3.2 SSR distribution in the M. japonicus genome

The general information on the SSR-containing sequences is summarized in Table 4. The total number of scaffold sequences with SSRs was 289 962, among which 126 210 contained more than one SSR. The lengths of SSR-containing sequences ranged from 201 bp to 15 960 bp, with an average of 468.1 bp. A total of 396 061 SSRs were detected, including 69 116 present in compound formation. Among the 326 945 perfect SSRs, dinucleotide repeats were dominant (44.08%), followed by mononucleotide repeats (29.67%), trinucleotide repeats (18.96%), and tetranucleotides (5.66%). Interestingly, we found that the hexanucleotide repeats (1.07%) were two-fold higher than pentanucleotide repeats (0.56%)(Fig. 1).

A more detailed investigation of individual repeat types was performed and is presented in Fig. 2. The results showed that the AT-rich repeats were more abundant than GC-rich repeats. The most frequent mono-, di-, and trinucleotide repeat motifs were A (88.05%), AG (54.12%), and AAT (29.14%), respectively (Fig. 2a); however, the motifs G (11.95%), CG (0.21%), and CCG (0.49%) were very rare. ATrich motifs were also predominant in tetra-, penta-, and hexanucleotide repeats (Fig. 2b), e.g., AAAC and AAAG were the most common tetranucleotides and AAGAG and AAAAG were the dominant hexanucleotides. However, the hexanucleotide repeats had a relatively lower degree of preference for A/T than the motifs of other sizes (Fig. 2b).

3.3 SSR distribution in the M. japonicus genome

The number of repeat motifs within the same and among different size units were distributed differently. The details of the repeat number distributions of different size units (mono-, di-, tri-, tetra-, penta-, and hexanucleotides) are listed in Fig. 3 and Table 5. The results revealed that shorter motifs and AT-rich types had relatively larger repeat numbers. The number A/T in mononucleotide repeats ranged from 10 to 59, but that of C/G ranged from 10 to 26, and the most common repeat number (10-12) accounted for 62.81% of the total mononucleotide repeats (Fig. 3a). The number AG/CT and AC/GT repeats ranged from 6 to 29, but that of CG/CG ranged from 6 to 7, and the most common repeat number (6-9) accounted for 64.69% of the total dinucleotide repeats (Fig. 3b). The number of AAT/ATT repeats ranged from 5 to 16, but CCG/CGG was mainly distributed in repeats of five (Table 4). The most common repeat number (5-7) accounted for 80.11%, 82.36%, 79.92%, and 78.23% of the total tri-, tetra-, penta-, and hexanucleotide repeats, respectively (Fig. 3c, d).

3.4 Functional annotation of SSRs

Among all of the SSR-containing sequences, 58 256 exhibited a significant GO annotation with three categories. The biological process category covered most of the functional genes (51.12%), followed by the molecular function (24.70%) and cellular component (24.18%) categories. In the biological process category, the dominant subcategories were metabolic process (19.44%), single-organism process (17.40%), and cellular process (17.22%)(Fig. 4a). Among the molecular function terms, a main proportion of clusters assigned to binding (42.15%) and catalytic activity (37.16%)(Fig. 4b). In the cellular component category, the major subcategories were membrane (22.27%), cell (16.02%), and cell part (15.83%)(Fig. 4c).

3.5 Development and laboratory verification of genome-wide SSR markers

In the present study, primers were designed for the di- to hexanucleotide repeats to develop genomewide SSR markers in M. japonicus. With the exceptions of compound repeats and mononucleotide loci, primers were successfully designed for 63.60%, 72.05%, 67.03%, 64.46%, and 46.74% of the di-, tri-, tetra-, penta-, and hex-nucleotide loci, respectively, proving themselves to be promising candidates for either PCR amplification or PAL (Fig. 1).

A subset of 30 SSR primer pairs were synthesized and tested in 42 individuals from a natural population to assess the PAL availability rate. Among the primer pairs tested, 90.0% could be amplified with prominent PCR products of the expected size and 80.0% were polymorphic (Table 1). The hexanucleotide (47 alleles for four loci with an average of 11.75) and pentanucleotide (52 alleles for five loci with an average of 10.40) exhibited the highest polymorphism, followed by trinucleotide (56 alleles for six loci with an average of 9.33) and dinucleotide loci (43 alleles for five loci with an average of 8.60). However, tetranucleotide loci exhibited the lowest polymorphism (33 alleles for four loci with an average of 8.25). A total of 231 alleles were recorded from the 24 SSR loci, and the number of alleles ranged from 5 to 17(with an average of 9.63). The PIC values for the polymorphic SSR loci ranged from 0.692 to 0.912 with an average of 0.796.

4 DISCUSSION

Despite their abundance and wide geographical distribution, little effort has been devoted to decoding the M. Japonicus genome, which has limited the use of SSR markers in genomic studies and markerassisted selection. In the present study, we analyzed the genome-wide SSR distribution and frequency of 1 to 6 bp sequences in the M. japonicus genome via GSS. Our study represents the first step towards decoding the M. japonicus genome, and also enables the development of SSR markers for genetic analysis. An extensive set of 326 945 perfect SSRs were detected without previously knowing anything about the genome sequence. As expected, the majority of the SSR markers had no Gene Ontology assignment because the SSRs were from M. japonicus, but a total of 58 256 of the SSR-containing scaffold sequences still had significant GO hits. The newly identified SSRs would be useful for extending our current knowledge of M. japonicus genome organization, and for genetic mapping, genome-wide association studies, marker-aided selection, comparative genomic analyses, and population genetics.

In the kuruma prawn, dinucleotide repeats were dominant, which was consistent with previous studies (Tóth et al., 2000 ; Somridhivej et al., 2008 ; Wang et al., 2011 ; Xu et al., 2011 ; Iranawati et al., 2012 ; Ji et al., 2012 ; Thanh et al., 2014). However, some studies have reported that the most abundant repeat motifs were trinucleotides (Stàgel et al., 2008). The higher frequency of trinucleotides in other studies might have been because they were identified from EST sequences, trinucleotide repeats prevail in the proteincoding exons of all taxa (Tóth et al., 2000). This difference might also have arisen from using different species, search parameters, and algorithms (Cavagnaro et al., 2010 ; Biswas et al., 2014). Interestingly, the M. japonicus genome had particularly high counts of hexanucleotide compared to pentanucleotide repeats, which was consistent with other studies (Tóth et al., 2000 ; Xu et al., 2011 ; Castoe et al., 2012 ; Vukosavljev et al., 2012). The higher counts of hexanucleotides might be because they were selected and changed during the process of evolution, because variation in hexanucleotides leads to increased morphological and functional proteins (Wahba et al., 1963 ; Touriol et al., 2003 ; Turanov et al., 2009).

As shown in Fig. 2, there is a remarkable variation in the frequency of individual repeat patterns in M.japonicus. The base composition of SSR patterns in M. japonicus is strongly biased toward A/T in all unit sizes, which may be the result of a high A+T content (61.79%) in the M. japonicus genome. Additionally, Tóth et al.(2000) also reported that poly (A/T) tracts were more abundant than poly (C/G) sequences in some species (e.g., primates, rodents, other mammals, non-mammalian vertebrates, arthropods, Caenorhabditis elegans, plants, yeast, and other fungi). The most abundant dinucleotide motif in M. japonicus was AG (61.74%), which was consistent with previous studies (Nagy et al., 2007 ; Zeng et al., 2010). However, this contradicted observations in other reports, where AT was the predominant motif (Shirasawa et al., 2010 ; Barchi et al., 2011 ; Wang et al., 2011). AAT was the dominant trinucleotide motif in M. japonicus, which was consistent with the research of Shirasawa et al.(2010) and Tóth et al.(2000) ; however, it was in contrast to other studies, where AAC was the predominant motif (Barchi et al., 2011 ; Wang et al., 2011). The same situation was detected for tetra-, penta-, and hex-nucleotide repeats. Large variation in motif abundance among species has also been reported in other studies (Cruz et al., 2005 ; Tanguy et al., 2008). One of the reasons for the differences in dominant motifs might be that they were from different resources (Genome vs EST), because characteristic differences between species can be observed in intergenic regions and introns (Tóth et al.(2000)).

The common trends in this study were, that AT-rich motifs were dominant and the CG-rich motifs, such as CG (2.12%), CCG (1.02%), and CCCG (0.070%), were very rare in the M. japonicus genome, which has also been reported in other species (Tóth et al.(2000) ; Kumpatla and Mukhopadhyay, 2005 ; Barchi et al., 2011 ; Wang et al., 2011). Additionally, the AT-rich motifs had relatively longer repeat lengths than CGrich motifs (Fig. 3). One of the main causes for the rarity and short length of the CG-rich motifs might have been that the three hydrogen bonds between C/G enable more stability and reduce the chance of DNA polymerase slippage (Tóth et al.(2000)). Another reason might be that CG-rich motifs are common in exons with less variance, whereas other regions have fewer of them (Tóth et al.(2000) ; Subramanian et al., 2003). Moreover, CG-rich repeats might form higher structures such as hairpin loops, which might stabilize the slipped structure (Schlötterer and Tautz, 1992 ; Sinden, 1999). Importantly, the variation in SSR lengths may be a good indication of polymorphism.

Because mononucleotide repeats are not suitable for marker development (Cavagnaro et al., 2010), only di- to hexanucleotide repeats were considered for primer design; primers were successfully designed for 151 541 SSR loci. We obtained a high PCR amplification efficiency and a high level of polymorphic loci in this study (90% of primer pairs yielded specific products, 80% of which were polymorphic), which likely resulted from the long repeats (Table 1). Previous studies have indicated that the SSR loci with long repeats are more prone to mutations both in plants and animals (Cavagnaro et al., 2010 ; Biswas et al., 2014). Our results also allowed us to understand the possible relationships between the degree of polymorphism and particular features of microsatellites, from which we found that hex- and pentanucleotides are more polymorphic than tetranucleotide repeats. These results might have been biased by the selection of the SSR loci analyzed, because we selected loci with long repeats for verification. In a future study, we will validate more SSR markers for linkage mapping and association analysis. The many polymorphic loci now available could be used for evolutionary and population genetic research, as well as high-resolution chromosome linkage mapping studies of M. japonicus.

5 CONCLUSION

Overall, the present study has contributed a detailed characterization and development of genome-wide SSR markers in M. japonicus by GSS. The kuruma prawn genome was relatively rich in microsatellites, and dinucleotide repeats and AT-rich SSRs were prevalent. Primers were successfully designed for almost half of the SSRs and a surprisingly high percentage (80%) of the polymorphic markers were verified in the laboratory. A total of 58 256 SSRcontaining sequences had significant GO annotation, these are good functional molecular marker candidates for association studies and comparative genomic analysis. The newly identified genome-wide SSR markers enhance M. japonicus genomic resources and will facilitate many genetic and genomic studies.

6 ACKNOWLEDGEMENT

We thank Majorbio Pharm Technology Co., Ltd.(Shanghai, China) for assisting with the RAD-Seq technique and Novogene Bioinformatics Technology Co., Ltd.(Beijing, China) for data analysis.