Development and Characterization of Novel Microsatellite Markers for the Peach Fruit Moth Carposina sasakii (Lepidoptera: Carposinidae) Using Next-Generation Sequencing

The peach fruit moth Carposina sasakii is an economically important pest on dozens of fruits from Rosaceae and Rhamnaceae in Northeast Asia. We developed novel microsatellite markers for C. sasakii from randomly sequenced regions of the genome using next-generation sequencing. In total, 95,153 microsatellite markers were isolated from 4.70 GB genomic sequences. Thirty-five polymorphic markers were developed by assessing in 63 individuals from two geographical populations. The allele numbers ranged from 2 to 9 with an average value of 4.60 per locus, while the polymorphism information content ranged from 0.075 to 0.696 with an average value of 0.407. Furthermore, the observed and expected heterozygosity varied from 0.000 to 0.677 and 0.062 to 0.771, respectively. The microsatellites developed provide abundant molecular markers for investigating genetic structure, genetic diversity, and existence of host-plant associated biotypes of C. sasakii.


Introduction
The peach fruit moth Carposina sasakii Matsumura (Lepidoptera: Carposinidae), is an important orchard pest in Northeast Asia [1,2]. In China, this pest was distributed throughout the country, except for Tibet. Its larvae can inflict direct damage on dozens of fruits, including peach, apple, pear, jujube, wild jujube, apricot, hawthorn, and pomegranate [3][4][5][6][7] by boring into fruitage. Differences in number of generations, emergence time of overwintering and diapause generation were found among populations on different host species [8,9], likely leading to low gene flow among host-plant populations. Thus, several studies attempted to reveal the differentiation of those moths occurring on different host plants [10,11].
To date, three types of molecular marker have been used to examine the existence of host biotypes in C. sasakii. Using esterase isozyme,   [11] reported that there is nearly no differentiation in isozyme-spectra between C. sasakii collected from jujube and wild jujube; however, populations collected from above two hosts were obviously different from those collected from apple orchard. Using RAPD (random amplified polymorphic DNA) to compare the populations collected from six kinds of host plants, including apple, hawthorn, peach, apricot, jujube, and wild jujube, Xu and Hua revealed that there were remarkable genetic differentiation between populations from apricot and those from other hosts [12]. Although a recent study using one region of mitochondrial DNA (mtDNA) sequences of cytochrome coxidase subunit I (COI) found that there was no evidence for associations between the variation of populations and host plants, the genetic differentiation showed significant correlation with the geographical distance [13]. Varied genetic markers used in the studies obviously lead to different results. To address this issue, more polymorphic and stable molecular markers are required.
Microsatellite is a kind of special sequence comprised by tandem repeats of one to six nucleotides. It always has high polymorphism, widely dispersed in both coding and noncoding regions of all prokaryotic and eukaryotic genomes [14]. Due to their codominant inheritance, high polymorphism, easy detection by polymerase chain reaction (PCR), and broad distribution in the genome, microsatellites are widely used for population genetic studies [15][16][17].
The traditional approach of microsatellite development, such as an enriched library followed by gene cloning, is time-consuming and labor-intensive. New approaches based on next-generation sequencing can be a good alternative. With the advantage of this technology, it is possible to develop a huge number of microsatellites, which are capable of generating tens of millions of short DNA sequence reads at a relatively low cost [18][19][20][21].
In the present study, we aimed to isolate microsatellites for C. sasakii from randomly obtained genomic sequences. This is the first report of novel microsatellites for C. sasakii. The markers developed will be helpful in investigating genetic structure, genetic diversity, and existence of host-plant associated biotypes of C. sasakii.

Microsatellite Marker Development
We generated 4.70 GB paired-end (PE) sequences with read length of 300 base pairs (bp), including 15,725,132 reads from a 500 bp insert DNA library constructed by Illumina MiSeq system. Raw data sequences were submitted to the National Center for Biotechnology Information (NCBI) Short Read Archive under accession number SRP068817. After removing low quality reads using SolexaQA software [22], the remaining high-quality reads were assembled into 1,902,994 contigs by SOAPdenovo2 [23]. They were with mean size of 252 bp and N50 of 286 bp, which are much shorter compared to a similar study in Dorcus hopei (Coleoptera) (N50 = 1218) [24]. This might be due to the method of assembly and the coverage of sequencing reads. However, the number of primer pairs designed in our study is reasonable (totally 8074 primer pairs / 479 Mb), as in other studies [19,[24][25][26].
A total of 95,153 microsatellite loci were discovered using MSDB version 2.4.3 software (http://msdb.biosv.com/) from the assembled contigs, which will be provided upon request. The detected microsatellites included 54,559 (57.34%) dinucleotide, 34,957 (36.74%) trinucleotide, 5591 (5.88%) tetranucleotide, and 46 (0.05%) pentanucleotide repeats (Table 1). There are no hexanucleotide repeats found under our searching conditions of microsatellite loci (a minimum of 25, 5, 5, 5, 5 and 5 repeats were used to identify the mono-, di-, tri-, tetra-, penta-, and hexanucleotide motifs, respectively). Dinucleotide repeats are more than the higher order motif, which is in agreement with the previous report of Arthropoda in insects like Aphis glycines (Hemiptera) [27] and Coccinella septempunctata (Coleoptera) [28], and in species of Arachnida [29]. According to the distribution of microsatellite (Table 1), it seems that the quantity of loci decreases followed with the increase of corresponding motif repeats.
Sixty-four primer pairs designed according to the sequences that are flanking trinucleotide repeats were selected for initial validation in eight individuals. Of them, 35 loci have polymorphic amplifications, 16 loci were monomorphic and 13 primer pairs did not produce any visible amplicon. These polymorphic loci can serve as candidate markers for future research, such as genetic diversity and relatedness analysis of different populations.

Characteristics of Validated Microsatellite Loci
The polymorphic loci obtained were assessed with two C. sasakii natural populations, including 31 individuals from Beijing and 32 individuals from Hubei province, China (Table 2 and Table S1). The 35 microsatellite markers had allele numbers ranging from 2 to 9 with an average value of 4.60 per locus. The polymorphism information content (PIC) revealed a range from 0.075 to 0.696 with an average value of 0.407. The observed (H O ) and expected (H E ) heterozygosity ranged from 0.000 to 0.677 and 0.062 to 0.771, respectively. The inbreeding coefficient (F IS ) ranged from´0.240 to 1.00.
The significantly high F IS in locus CS21, CS38 and CS82 might be caused by the low H O , rather than sampling bias since most loci showed low F IS in the two populations. The loci CS31 and CS33 showed significant linkage disequilibrium only across Beijing population (corrected by Holm's correction, p < 0.05). It is speculated that the linkage disequilibrium observed at certain loci in some populations may be due to substructure of population or bottleneck [30]. Eight loci in Beijing population and 13 loci in Hubei population significantly deviated from Hardy-Weinberg equilibrium (HWE), while 5 loci (CS05, CS17, CS21, CS29 and CS82) showed significant value in the both tested populations. The loci deviated from HWE might be resulted by heterozygote deficiency, because H O is much lower than H E in these loci ( Table 2). Heterozygote deficiency can be caused by the Wahlund Effect [31] or the presence of null alleles, for which Lepidoptera species are notorious [32][33][34][35][36]. It was considered that the present of null alleles is very common in this order due to the flanking region with repetitive sequences and multiple copies of loci [37][38][39]. Random sequences of C. sasakii genome obtained by the Illumina MiSeq system may cover coding regions. Thus, they are probably linked to sites under selection, which cannot reflect facticity of population diversity and structure. A neutrality test was done with all of the 35 loci. Interestingly, all of the loci were under neutral expectations ( Figure 1). Therefore, deviating from HWE is not necessarily due to the characteristics of loci. It may imply the distinct population structure, biological property of the species, or just sampling error, e.g., examined individuals from the same egg brood can also lead to deviation from HWE [40].
The population structure of C. sasakii was inferred with the dataset of 35 microsatellite markers. The 63 individuals from two geographic populations were divided into two clusters. As can be clearly seen in Figure 2, there are genetic differences between two populations, indicating that the microsatellite markers validated could be used to discriminate geographic populations and other genetic study of C. sasakii. showed significant linkage disequilibrium only across Beijing population (corrected by Holm's correction, p < 0.05). It is speculated that the linkage disequilibrium observed at certain loci in some populations may be due to substructure of population or bottleneck [30]. Eight loci in Beijing population and 13 loci in Hubei population significantly deviated from Hardy-Weinberg equilibrium (HWE), while 5 loci (CS05, CS17, CS21, CS29 and CS82) showed significant value in the both tested populations. The loci deviated from HWE might be resulted by heterozygote deficiency, because HO is much lower than HE in these loci (Table 2). Heterozygote deficiency can be caused by the Wahlund Effect [31] or the presence of null alleles, for which Lepidoptera species are notorious [32][33][34][35][36]. It was considered that the present of null alleles is very common in this order due to the flanking region with repetitive sequences and multiple copies of loci [37][38][39]. Random sequences of C. sasakii genome obtained by the Illumina MiSeq system may cover coding regions. Thus, they are probably linked to sites under selection, which cannot reflect facticity of population diversity and structure. A neutrality test was done with all of the 35 loci. Interestingly, all of the loci were under neutral expectations ( Figure 1). Therefore, deviating from HWE is not necessarily due to the characteristics of loci. It may imply the distinct population structure, biological property of the species, or just sampling error, e.g., examined individuals from the same egg brood can also lead to deviation from HWE [40]. The population structure of C. sasakii was inferred with the dataset of 35 microsatellite markers. The 63 individuals from two geographic populations were divided into two clusters. As can be clearly seen in Figure 2, there are genetic differences between two populations, indicating that the microsatellite markers validated could be used to discriminate geographic populations and other genetic study of C. sasakii.   showed significant linkage disequilibrium only across Beijing population (corrected by Holm's correction, p < 0.05). It is speculated that the linkage disequilibrium observed at certain loci in some populations may be due to substructure of population or bottleneck [30]. Eight loci in Beijing population and 13 loci in Hubei population significantly deviated from Hardy-Weinberg equilibrium (HWE), while 5 loci (CS05, CS17, CS21, CS29 and CS82) showed significant value in the both tested populations. The loci deviated from HWE might be resulted by heterozygote deficiency, because HO is much lower than HE in these loci ( Table 2). Heterozygote deficiency can be caused by the Wahlund Effect [31] or the presence of null alleles, for which Lepidoptera species are notorious [32][33][34][35][36]. It was considered that the present of null alleles is very common in this order due to the flanking region with repetitive sequences and multiple copies of loci [37][38][39]. Random sequences of C. sasakii genome obtained by the Illumina MiSeq system may cover coding regions. Thus, they are probably linked to sites under selection, which cannot reflect facticity of population diversity and structure. A neutrality test was done with all of the 35 loci. Interestingly, all of the loci were under neutral expectations ( Figure 1). Therefore, deviating from HWE is not necessarily due to the characteristics of loci. It may imply the distinct population structure, biological property of the species, or just sampling error, e.g., examined individuals from the same egg brood can also lead to deviation from HWE [40]. The population structure of C. sasakii was inferred with the dataset of 35 microsatellite markers. The 63 individuals from two geographic populations were divided into two clusters. As can be clearly seen in Figure 2, there are genetic differences between two populations, indicating that the microsatellite markers validated could be used to discriminate geographic populations and other genetic study of C. sasakii.

Sample Collection and DNA Extraction
A total of 63 larvae were collected from two geographic regions in China, of which 31 samples were from Yanqing of the Beijing (N 40˝27 1 20.05", E 115˝58 1 8.14"), named BJYQ, and 32 specimens came from Yichang of Hubei province (N 30˝41 1 39.43", E 111˝16 1 50.77"), named HBYC. Additionally, eight individuals from eight sampled sites were used for the initial test. Samples were stored in ethanol absolute and frozen at´80˝C prior to use. Genomic DNA were extracted from half of an individual larva using DNeasy Blood & Tissue Kit (QIAGEN, Hilden, Germany), according to the manufacturer's instructions.

Sequencing, Microsatellites Searching and Primer Design
One larva of C. sasakii from Beijing was used to prepare the library with the Illumina TruSeq DNA PCR-Free HT Library Prep Kit (Illumina, San Diego, CA, USA), and then sequenced on a Illumina MiSeq Sequencer using the MiSeq Reagent Kit v3 (Illumina, San Diego, CA, USA). Generated genomic sequences were assembled by SOAPdenovo program [23].
The microsatellite isolation from the genomic sequences and primer design for loci was conducted in the software QDD [41]. The searching criteria were as follows: at least six motif repeats for target microsatellites, and PCR product lengths ranged between 90 and 350 bp. For primer design, the annealing temperature ranged from 52 to 68˝C, and the difference in annealing temperature in one pairwise primer was <5˝C. The remaining parameters were at default settings.

Primer Testing and Polymorphism Detection
Firstly, in order to improve efficiency and lower cost, we added a PC tail (Primer tail C) (5 1 CAGGACCAGGCTACCGTG 3 1 ) to the 5 1 end of the candidate forward primer [42]. Eight larvae of C. sasakii from eight different populations were used for the initial test. Amplification was carried out in a final volume of 10 µL, containing 0. The amplification program was as follows: 4 min at 94˝C; 35 cycles of 30 s at 94˝C, 30 s at 56˝C, and 45 s at 72˝C, with a final 10-min extension at 72˝C. PCR products were visualized on agarose gel (1.5%) electrophoresis. This step was taken to screen primers that can amplify PCR fragment.
Secondly, primers selected in previous steps were tested using a capillary sequencer. Amplification was performed in a final volume of 10 µL, containing 0.5 µL (12.5 ng) of template DNA, 5 µL of Master Mix (Promega, Madison, WI, USA), 0.08 µL of forward primer (modified by the PC tail) at a final concentration 0.08 µM, 0.16 µL of reverse primer at a final concentration 0.16 µM, 0.32 µL of PC tails modified by fluorescence (FAM (blue), HEX (green), and ROX (red)) including different color at a final concentration 0.32 µM, and 3.94 µL of ddH 2 O. The amplification program was the same as above. The ABI 3730xl DNA Analyzer (Applied Biosystems, Foster, CA, USA) was used to analyze the amplified PCR fragments with the GeneScan 500 LIZ size standard (Applied Biosystems).
Finally, marker primers screened out by the first two steps were validated in 63 samples from two regions. Amplification mixture, amplification program, and analysis of PCR fragments were the same as the second step.

Statistical Analysis
Genotyping data was identified, and errors were corrected by MICRO-CHECKER [43]. Diversity statistics including allele frequencies, Ho, He and PIC were estimated by the macros Microsatellite Tools [44]. Tests for linkage disequilibrium among loci within each population and deviation from HWEat each locus/population pair, and estimation of F IS for each population, were performed in GENEPOP v4.0 (Applied Biosystems). Additionally, the null allele test was conducted with FREENA [34]. The program LOSTAN [45] was used to detect putative loci potentially under selection with two options: neutral mean F ST ' and force mean F ST '. Corresponding sequences of polymorphic loci were screened using BLASTx and BLASTn in the NCBI database (http://www.ncbi.nlm.nih.gov/). Population differentiation was investigated using the Bayesian clustering approach implemented in the program STRUCTURE, version 2.3.3 [46]. Simulations were run for 200,000 Markov chain Monte Carlo with a burn-in of 100,000 iterations under admixture ancestry and correlated allele frequency models. We performed 15 independent runs for each K (from 1 to 6) to confirm consistency across runs. The most accurate number of groups (K) was visually examined when plotting K against delta-K and using the Evanno method in the online program STRUCTURE HARVESTER [47].

Conclusions
We characterized and developed microsatellite markers for C. sasakii from random regions of the genome generated by using next-generation sequencing. The loci assessed in our study could reveal the genetic structure in two geographical populations. This method provides fast way for high throughput development of microsatellite markers from non-model species without reference genome.