Development of 12 Microsatellite Markers in Dorcus titanus castanicolor (Motschulsky, 1861) (Lucanidae, Coleoptera) from Korea Using Next-Generation Sequencing

In the present study, we used next-generation sequencing to develop 12 novel microsatellite markers for genetic structural analysis of Dorcus titanus castanicolor (Lucanidae; Coleoptera), a popular pet insect in China, Korea, and Japan. We identified 52,357 microsatellite loci in 339,287,381 bp of genomic sequence and selected 19 of the loci based on their PCR amplification efficiency and polymorphism. The 19 selected markers were then tested for the presence of null alleles and linkage disequilibrium. We did not detect any evidence of null alleles; however, four pairs of loci (DT03 and DT11, DT05 and DT26, DT08 and DT26, DT26 and DT35) exhibited linkage disequilibrium. Thus, we assessed the genetic diversity of a D. titanus castanicolor population from the Daejeon region of Korea (n = 22) using 13 markers. Among them, one marker (DT17) deviated from Hardy-Weinberg equilibrium. Therefore, 12 markers may be useful for further analyzing the genetic diversity of D. titanus castanicolor.


Introduction
Dorcus titanus castanicolor [1] is a subspecies of stag beetle that is distributed in northeast China, Korea, and Tsushima Island, Japan [2][3][4] and is distinguished from the other subspecies by a distinctly bilobed clypeus and flattened last abdominal sternite in males and moderately punctured elytra and straight protibiae in females [4] (Figure 1). In Korea, the subspecies is distributed throughout the country and is commonly found during local fauna and monitoring surveys [3,4].
Marketing campaigns from the insect industry have contributed to a rapid increase in the number of pet insects [5][6][7]. D. titanus castanicolor is a popular pet insect [7], and the rearing method of the subspecies is well documented [8]. In China and Japan, many stag beetle breeders have reared large-sized individuals of D. titanus [9]. Japan has 12 subspecies of D. titanus, which were interbred to select for large individuals. Consequently, four foreign subspecies (D. titanus platymelus, D. titanus pilifer, D. titanus palawanicus, and D. titanus titanus) were introduced to Japan without prior ecological assessment. Thus, ecological problems, such as decreases in the populations of native subspecies, have occurred, owing to competition with the introduced subspecies [9]. populations of native subspecies, have occurred, owing to competition with the introduced subspecies [9]. Recently, several Southeast Asian subspecies, such as D. titanus palawanicus, were introduced into Korea. Kim and Kim [4] conducted a comparative analysis of the male genitalia through a taxonomic review of the Korean Lucanidae and reported that the Korean subspecies (D. titanus castanicolor) might actually be a distinct species. A different case was reported for Dorcus hopei [10]. Thus, seven microsatellite markers were developed in order to help conserve the Korean population and used to compare the genetic structure of the Korean Dorcus hopei with other regional populations [11]. Recently, the Korean Quarantine Agency also used the markers for a quarantine inspection of foreign Dorcus hopei. Therefore, the goal of this study was to develop microsatellite markers for D. titanus castanicolor that could be used for population genetic analysis and origin testing.

Next-Generation Sequencing and Microsatellite Loci Identification
Through Illumina sequencing of D. titanus castanicolor genomic DNA, the number of contigs with length >1000 bp was 123,291. The maximum length was 103,231 bp, and the mean length was 2729 bp. We obtained 339,287,381 bp of genomic DNA sequence with a GC content of 38.98% (Table 1). Searching the genomic DNA sequence yielded 52,357 microsatellite loci with 14,644 di-, tri-, and tetra-repeat sequence motifs, 186 penta-repeat loci, and 37 hexa-repeat loci. The relative Recently, several Southeast Asian subspecies, such as D. titanus palawanicus, were introduced into Korea. Kim and Kim [4] conducted a comparative analysis of the male genitalia through a taxonomic review of the Korean Lucanidae and reported that the Korean subspecies (D. titanus castanicolor) might actually be a distinct species. A different case was reported for Dorcus hopei [10]. Thus, seven microsatellite markers were developed in order to help conserve the Korean population and used to compare the genetic structure of the Korean Dorcus hopei with other regional populations [11]. Recently, the Korean Quarantine Agency also used the markers for a quarantine inspection of foreign Dorcus hopei. Therefore, the goal of this study was to develop microsatellite markers for D. titanus castanicolor that could be used for population genetic analysis and origin testing.

Next-Generation Sequencing and Microsatellite Loci Identification
Through Illumina sequencing of D. titanus castanicolor genomic DNA, the number of contigs with length >1000 bp was 123,291. The maximum length was 103,231 bp, and the mean length was 2729 bp. We obtained 339,287,381 bp of genomic DNA sequence with a GC content of 38.98% (Table 1). Searching the genomic DNA sequence yielded 52,357 microsatellite loci with 14,644 di-, tri-, and tetra-repeat sequence motifs, 186 penta-repeat loci, and 37 hexa-repeat loci. The relative frequencies of the specific repeated sequences were as follows: ATT, 24%; AAT, 23%; AT, 8%; and remaining motifs, below 5% ( Figure 2). frequencies of the specific repeated sequences were as follows: ATT, 24%; AAT, 23%; AT, 8%; and remaining motifs, below 5% ( Figure 2). The 14,644 searched microsatellite loci were filtered with two conditions. First, the loci composed of AT repeat were excluded because of difficulty in PCR. Second, the loci composed of 24-96 bases were selected. Using these two filtering conditions, we selected 316 microsatellite loci, and using the primer design software PRIMER3 version 0.4.0 [12,13], we designed primer sets for 301 loci. Eighty-five markers were selected based on the GC contents of the product sequences and melting temperatures of the designed primer sets and were subsequently tested for specificity and polymorphism. Finally, 19 markers were selected and the forward primers of each selected loci were labeled with a fluorescent dye (FAM) at the 5′ end [14] (Table 2).  The 14,644 searched microsatellite loci were filtered with two conditions. First, the loci composed of AT repeat were excluded because of difficulty in PCR. Second, the loci composed of 24-96 bases were selected. Using these two filtering conditions, we selected 316 microsatellite loci, and using the primer design software PRIMER3 version 0.4.0 [12,13], we designed primer sets for 301 loci. Eighty-five markers were selected based on the GC contents of the product sequences and melting temperatures of the designed primer sets and were subsequently tested for specificity and polymorphism. Finally, 19 markers were selected and the forward primers of each selected loci were labeled with a fluorescent dye (FAM) at the 5' end [14] (Table 2).

Microsatellite Marker Assessment
The genetic diversity of the Korean D. titanus castanicolor (n = 22) was assessed, using the 19 selected markers. We determined the genotypes of 22 individuals using a 3730XL DNA analyzer (Applied Biosystems ® , Foster City, CA, USA) (Table S1). Each marker was tested for PCR errors and the presence of null alleles (Table 3). Then, the linkage disequilibrium analysis was conducted for each locus. From these results, linkage disequilibrium was detected in four pairs of loci (DT03 and DT11, DT05 and DT26, DT08 and DT26, and DT08 and DT35), which were excluded from further analysis. Thus, we considered that the 13 remaining microsatellite loci might be suitable for analyzing the population genetic structure of Korean D. titanus castanicolor.

Genetic Diversity of the Korean Dorcus titanus castanicolor Population
Genetic diversity metrics for the 13 selected markers were calculated using PowerMarker version 3.25 [15] ( Table 4). The genotype number, allele number, expected heterozygosity, and observed heterozygosity ranged from 3-11, 3-9, 0.3140-0.7789, and 0.3636-1.0000, respectively. Significant departure of genotype frequencies from Hardy-Weinberg equilibrium (HWE) was determined after sequential Bonferroni corrections for the significance level (critical p = 0.003). We found that only one locus (DT17) significantly deviated from HWE (Table 4). Increasing the benefit as a pet insect, the wild population of the Korean D. titanus castanicolor is decreasing because of over-collecting by the Korean breeders. Some Korean breeders surreptitiously made crosses between the Korean and exotic subspecies (such as D. titanus palawanicus) for body size and mandible morphology, like in Japan [9]. For these reasons, we anticipated that there could be multiple non-HWE markers. Although our data did not suggest this, further surveys are necessary using the 12 markers and additional populations of Korean D. titanus castanicolor.

Sample Collection and Genomic DNA Extraction
We collected 23 Dorcus titanus castanicolor specimens from Mt. Jangtae, which is located in Jangan-dong, Seo-gu, Daejeon, Korea (Table S2). We extracted the thoracic muscles of the specimens for genomic DNA isolation and deposited the dried samples in the National Institute of Biological Resources (Incheon, Korea). Genomic DNA for microsatellite loci amplification and genotyping was isolated using a DNA purification kit (Prime Prep Genomic DNA Isolation Kit, GeNet Bio, Daejeon, Korea) according to the manufacturer's instructions. Genomic DNA for next-generation sequencing was extracted from one individual using the NucleoSpin ® Tissue Kit (Macherey-Nagel GmbH and Co. KG, Düren, Germany), and its quality was checked, using both spectrophotometry and electrophoresis on a 1% agarose gel.

Next-Generation Sequencing, Microsatellite Loci Identification, Marker Selection, and Genotyping
For next-generation sequencing, 10 µg of genomic DNA from a single specimen was used for 2 × 300 paired-end sequencing with a HiSeq Sequencer (Illumina, San Diego, CA, USA). Assembly was performed using the de Bruijn graph algorithm in CLC Genomics Workbench version 7 (CLC Bio, Aarhus, Denmark). Assembly mapping was performed after removing Illumina adapters and low-quality sequences using the CLC trimmer function (default limit = 0.05). The assembly procedure used the parameters Length Fraction (LF) and Sequence Similarity (SIM) between DNA reads, as described by the CLC Genomics Workbench software, with maximum stringency (0.50 LF and 0.80 SIM) and a minimum contig length of 70 bp.
Microsatellite loci were identified in the assembled partial draft genome using tandem repeat search software, Phobos version 3.3.12 [16,17]. The searching criteria were 10% mismatch di-, tri-, and tetra-nucleotide repeats with a sequence length between 24 and 96 bp, and AT repeats were excluded from the screening results because of difficulty in PCR. Primer sets for amplifying of the microsatellite loci were designed using PRIMER3 version 0.4.0 [12,13] with the following criteria: GC content >30%, primer length 18-22 bp, optimal temperature of 56-60 • C, and default settings for the remaining parameters. PCR for the qualification test on the designed primer sets was conducted in 20 µL reactions using AccuPower PCR PreMix (Bioneer, Daejeon, Korea), 30 ng template DNA, and 10 pmol each primer. Extra MgCl 2 was not added. PCR amplification was carried out in a Veriti 96-Well Thermal Cycler (Applied Biosystems ® , Foster City, CA, USA), using the following program: an initial denaturation step of 5 min at 94 • C; followed by 35 cycles of 10 s at 94 • C, 10 s at 56 or 60 • C, and 20 s at 72 • C; and a final extension step of 5 min at 72 • C. The specificity of the primer sets and polymorphism of the microsatellite loci were examined through electrophoresis of the PCR amplicon using QIAxcel (Qiagen, Leipzig, Germany). Each forward primer of the selected markers was labeled with 6-carboxyfluorescein at the 5 end for genotyping [14]. To reduce polymerization error, high fidelity polymerase (SuPrime HF Premix (2×), GeNet Bio, Daejeon, Korea) was used in PCR for genotyping, which was carried out using the same conditions as those for the qualification test of the designed primer sets. The sequences of microsatellite loci for the developed markers were submitted to NCBI GenBank (Table 3).

Data Analysis
The presence of null alleles and PCR errors in the genotyping results was checked using Micro-Checker version 2.2.3 [18], and pairwise linkage disequilibrium was analyzed using Arlequin version 3.1 [19]. The genetic diversity of the Daejeon population was estimated using the major allele frequency (MAF), number of Genotypes (Gn), sample size (SS), number of alleles (An), expected heterozygosity (He), observed heterozygosity (Ho), polymorphism information content (PIC), and Hardy-Weinberg equilibrium (HWE), all of which were calculated using PowerMarker version 3.25 [15]. Deviation from HWE was determined after sequential Bonferroni correction for the significance level.