Comprehensive Analyses of the Complete Mitochondrial Genome of Figulus binodulus (Coleoptera: Lucanidae).

Figulus binodulus Waterhouse is a small stag beetle distributed in East Asia. We determined the first mitochondrial genome of F. binodulus of which is 16,261-bp long including 13 protein-coding genes, two ribosomal RNA genes, 22 transfer RNAs, and a single large noncoding region of 1,717 bp. Gene order of F. binodulus is identical to the ancestral insect mitochondrial gene order as in most other stag beetle species. All of 22 tRNAs could be shaped into typical clover-leaf structure except trnSer1. Comparative analyses of 21 Lucanidae mitochondrial genomes was conducted in aspect of their length and AT-GC ratio. Nucleotide diversities analyses provide that cox1 and cox2 in Lucanidae are less diverse than those of Scarabaeoidea. Fifty simple sequence repeats (SSRs) were identified on F. binodulus mitochondrial genome. Comparative analysis of SSRs among five mitochondrial genomes displayed similar trend along with SSR types. Figulus binodulus was sister to all other available family Lucanidae species in the phylogenetic tree.

Figulus binodulus Waterhouse is a small stag beetle species of the tribe Figulini found in east Asian countries, such as Korea, Japan, China, Taiwan, and Vietnam (Hangay and De Keyzer 2017). They have a rather unusual life cycle for a stag beetle: While most stag beetles are herbivorous throughout their life time, the adult F. binodulus turn carnivorous, preying on beetle larvae or other insects living in the decaying wood (Mori and chiba 2009). Moreover, F. binodulus is also unique in that they live in small groups, where the adults pulverize the wood for the young ones, making it easier for the larvae to consume wood (Mori and Chiba 2009). As a result of this subsocial life style, they rarely leave the dead wood, with the only exceptions seen in the breeding seasons (Mori and Chiba 2009).

One of the critical issues in scientific researches is disproportionate research efforts. Only some groups of insects, such as Papilionidae, Lucanidae, and Cicindelinae, are well researched in that aspect (Stork 2018). Even in such well-examined taxa, the uneven study interest continues. For instance, while there are 20 complete mitochondrial genomes of Lucanidae available in the NCBI (As of February 2020), 15 of them originate from subfamily Lucaninae covering only three tribes, Dorcini, Lucanini, and Aegini (Table 2).

Table 2.

List of 21 available Lucanidae mitochondrial genomes including Figulus binodulus used in this study

Sub-family	Species name	NCBI accession	Total length (bp)	Number of PGCs	Number of tRNAs	Number of rRNAs	GC ratio (%)	AT ratio (%)	GC Skew	AT Skew	Reference
Lucaninae	Figulus binodulus	MN180051	16,261	13	22	2	30.71	69.29	−0.3206	0.1125	This study
	Dorcus hansi	NC_043928	18,130	13	22	2	29.62	70.38	−0.3151	0.0552	Unpublished
	Dorcus hopei	MF612067	16,026	13	22	2	31.95	68.05	−0.3043	0.0915	Chen et al., 2018b
	Dorcus seguyi	NC_038212	17,953	13	22	2	28.78	71.22	−0.3066	0.0612	Chen et al., 2018b
	Dorcus seguyi	MF612069	18,503	13	22	2	28.76	71.24	−0.3031	0.0607	Chen et al., 2018b
	Dorcus parallelipipedus	KT876887	17,561	13	22	2	30.88	69.12	−0.2763	0.0662	Linard et al., 2016
	Lucanus cervus	MH595464	22,714	13	22	2	30.10	69.90	−0.2536	0.0802	Chen et al., 2019
	Lucanus mazama	NC_013578	15,258	13	22	2	32.87	67.13	−0.2723	0.0744	Sheffield et al., 2009
	Lucanus sp.	KT876903	20,631	13	22	2	30.16	69.84	−0.2663	0.0683	Linard et al., 2016
	Neolucanus maximus	NC_039652	16,601	13	22	2	32.09	67.91	−0.3321	0.1494	Linard et al., 2016
	Odontolabis cuvera fallaciosa	MF908524	19,614	13	22	2	29.50	70.50	−0.2572	0.1335	Wang et al., 2018
	Prosopocoilus astacoides blanchardi	KF364622	21,628	13	22	2	32.99	67.01	−0.2655	0.1069	Kim et al., 2015
	Prosopocoilus confucius	NC_036038	16,951	13	22	2	32.15	67.85	−0.3616	0.0767	Lin et al., 2017
	Prosopocoilus gracilis	NC_027580	16,736	13	22	2	33.91	66.09	−0.3318	0.1062	Wu et al., 2016
	Rhaetus westwoodii	MG159815	18,131	13	22	2	32.50	67.50	−0.3222	0.0755	Jing et al., 2018
	Serrognathus platymelus	NC_044096	16,790	13	22	2	31.23	68.77	−0.2837	0.0785	Unpublished
Syndesinae	Ceruchus minor	NC_042613	18,601	13	23	2	31.61	68.39	−0.2799	0.1292	Unpublished
	Sinodendron rugosum	NC_042614	18,126	13	22	2	26.97	73.03	−0.2979	0.0638	Unpublished
	Sinodendron yunnanense	NC_036157	16,921	13	22	2	24.94	75.06	−0.2706	0.0347	Lin et al., 2017
Aesalinae	Aesalus sp.	MH120282	17,743	13	22	2	26.48	73.52	−0.2129	0.0672	Unpublished
	Himaloaesalus gaoligongshanus	NC_042922	17,013	13	22	2	25.61	74.39	−0.2564	0.1078	Unpublished

To extend our understanding of family Lucanidae especially in the aspect of mitochondrial genomes, we completed the first complete mitochondrial genome of F. binodulus in the tribe Figulini belongs to family Lucanidae. We compared the sequence with all 20 available mitochondrial genomes of Lucanidae in various aspects, including transfer RNA structure, mitochondrial genome configuration, nucleotide diversity throughout complete mitochondrial genomes, simple sequence repeats (SSRs), and phylogenetic analyses. Comparative data generated in this study would be useful for further understanding of phylogenetic relationship (Cameron et al. 2009, Li et al. 2019, Liu et al. 2018), for developing molecular markers to distinguish species or even populations within species based on nucleotide diversity and SSRs (Simon et al. 1994, Mousson et al. 2005), and for identifying cryptic species (Burger et al. 2014).

Materials and Methods

Sample Preparation and DNA Extraction of F. binodulus

Total DNA of F. binodulus was extracted from an adult individual collected in Gageodo, Jeollanam province (34°03′04.0″ N, 125°07′48.6″ E), Republic of Korea, using DNeasy Blood &Tissue Kit (QIAGEN, Hilden, Germany). DNA sample and specimen (95% ethanol) are deposited in the InfoBoss Cyber Herbarium (IN; J. Lee, INH-00021).

Genome Sequencing and de novo Assembly of Mitochondrial Genome of F. binodulus

Raw sequences were obtained from Illumina HiSeqX with constructing 350-bp insertion pair-end library at Macrogen Inc., Korea. Under the environment of Genome Information System (GeIS; http://geis.infoboss.co.kr/; Park et al., in preparation), raw sequence were filtered by Trimmomatic v0.33 (Bolger et al. 2014) and subjected to de novo assembly process done by Velvet v1.2.10 (Zerbino and Birney 2008) with k-mers ranging from 61 to 75 in order to obtain the complete mitochondrial genome sequences. Filling gap sequences as well as circular test were conducted with SOAPGapCloser v1.12 (Zhao et al. 2011). After that, all assembled bases were manually investigated using BWA v0.7.17 (Li et al. 2009) and SAMtools v1.9 (Li 2013) to correct misassembled sequences bases.

Mitochondrial Genome Annotation

Geneious R11 11.1.5 (Biomatters Ltd, Auckland, New Zealand) was used to annotate the mitochondrial genome based on sequence alignment with other Lucanid mitochondrial genomes. To confirm location and structure of transfer RNAs (tRNAs), the annotated GenBank format file of F. binodulus mitochondrial genome was subjected to the MITOS web server with genetic code ‘05-invertebrate’ (Bernt et al. 2013), and ARWEN server with default option (Laslett and Canbäck 2008). The prediction results were reviewed manually and drawn into final tRNA structure. The annotated GenBank format file of F. binodulus mitochondrial genome was used to draw the circular map using CGView with default options (Grant and Stothard 2008)

Identification of SSRs on F. binodulus Mitochondrial Genome

SSRs were identified on the chloroplast genome sequence using the pipeline of the SSR database (SSRDB; http://ssr.pe.kr/; Park et al., in preparation). Based on conventional definition of SSR on organelle genomes: monoSSR (unit sequence length is 1 bp) to hexaSSR (unit sequence length is 6 bp), that over 10-bp long. Since various criteria of SSRs was used for organelle genomes (Gandhi et al. 2010, Chen et al. 2015, Cheng et al. 2016, Shukla et al. 2018, Jeon and Kim, 2019, Li et al. 2019), we adopted the criteria used in organelle genomes of Dysphania ambrosioides (Kim et al. 2019) and Arabidopsis thaliana (Park et al. 2020b) monoSSR (unit sequence length is 1 bp) to hexaSSR (6 bp) are used as normal SSRs and heptaSSR (7 bp) to decaSSR (10 bp) were defined as extendedSSRs. Among normal SSRs, pentaSSRs, and hexaSSRs of which unit number was 2 were classified as potentialSSRs.

Nucleotide Diversity Analysis of F. binodulus Mitochondrial Genome

Nucleotide diversity of the 21 Lucanid mitochondrial genomes was calculated based on the method proposed by Nei and Li (Nei and Li 1979) using the perl script, one of analysis tools implemented in the GenomeArchive; http://www.genoomearchive.info/ (Park and Xi 2018, Park et al., in preparation). Window size and step size of sliding-window method were set as 500 and 200 bp, respectively. Genomic positions of each windows were compared with gene annotations of the mitochondrial genome.

All available 88 Scarabaeoid mitochondrial sequences including the 21 Lucanids mitochondrial genomes that contained all 13 protein-coding gene (PCGs) were retrieved from NCBI and sequences of each PCGs were extracted. Multiple sequence alignments for each PCG sequence for two datasets (21 Lucanids mitochondrial genomes and 88 Scarabaeiod mitochondrial genomes) were conducted with MAFFT v7.450 (Katoh and Standley 2013). Nucleotide diversity for each alignment was calculated based on Nei and Li (1979) using the perl script without sliding-window option.

Construction of Phylogenetic Trees

Thirteen PCGs and 2 ribosomal RNAs (rRNAs) were extracted from 21-stag beetle mitochondrial genomes and an outgroup cockchafer species (Rhopaea magnicornis; NC_027602). The 15 genes were first aligned individually using MAFFT v7.450 (Katoh and Standley 2013), then were concatenated to construct the phylogenetic trees. Maximum likelihood (number of bootstrap repeats is 1,000) and neighbor-joining (number of bootstrap repeats is 10,000) phylogenetic trees were constructed using MEGA X (Kumar et al. 2018). During the ML analysis, a heuristic search was used with nearest-neighbor interchange (NNI) branch swapping, the Tamura-Nei model, and uniform rates among sites. All other options were set to their default values. Bootstrap analyses with 1,000 pseudoreplicates were conducted with the same options. Bayesian inference (number of generations is 1,100,000) tree was constructed by Mr. Bayes v3.2.6 (Huelsenbeck and Ronquist 2001) under the environment of Geneious R11 1.1.5 (Biomatters Ltd, Auckland, New Zealand). The GTR model with gamma rates was used as a molecular model. A Markov-chain Monte Carlo (MCMC) algorithm was employed for 1,000,000 generations, sampling trees every 200 generations, with four chains running simultaneously. Trees from the first 100,000 generations were discarded as burn-in.

Results and Discussions

Complete Mitochondrial Genome of F. binodulus

Figulus binodulus mitochondrial genome (GenBank MN180051) was 16,261-bp long and contained 13 PCGs, 2 rRNAs, and 22 tRNAs, which is the typical configuration of an insect mitochondrial genome (Figure 1; Table 1). Its GC ratio is 30.71% which was within the range of previously known species with the highest found in Prosopocoilus gracilis (33.91%) and the lowest in Sinodendron yunnanense (24.94%; Table 2). The order of the 37 genes found in F. binodulus mitochondrial genome was identical to the ancestral insect mitochondrial gene order (Cameron 2014). This order was conserved throughout all stag beetle lineages with only some minor modifications: 1) trnTyr duplication was found in Cherucus minor and 2) relocation of trnLeu2 was observed in Sinodendron species (Lin et al. 2017).

Fig. 1

Complete mitochondrial genome sequence of Figulus binodulus. Black graph inside circular diagram presents GC ratio of mitochondrial genomes. Colorful bars on outer circular form indicates CDSs (blue), tRNAs (light brown), rRNAs (pink), and control region (gray).

Table 1.

Summary of Figulus binodulus mitochondrial genome

Gene	Strand	Start (bp)	End (bp)	Size (bp)	Start codon	End codon	Anticodon	Intergenic length (bp)
trnIle	J	2	63	62	ACA	TA	GAU	N/A
trnGln	N	61	129	69	TAC	TAA	UUG	−3
trnMet	J	129	196	68	AAA	TA	CAU	−1
nad2	J	197	1,207	1011	ATG	TAA		0
trnTrp	J	1,218	1,282	65	AAG	TA	UCA	10
trnCys	N	1,275	1,335	61	AGC	T	GCA	−8
trnTyr	N	1,335	1,398	64	GGT	A	GUA	−1
cox1	J	1,400	2,930	1531	AAC	T		1
trnLeu	J	2,931	2,995	65	TCT	AA	UAA	0
cox2	J	2,996	3,680	685	ATA	T		0
trnLys	J	3,681	3,751	71	CAT	GA	CUU	0
trnAsp	J	3,751	3,812	62	AAA	TA	GUC	−1
atp8	J	3,813	3,968	156	ATT	TAA		0
atp6	J	3,962	4,630	669	ATG	TAA		−7
cox3	J	4,630	5,413	784	ATG	T		−1
trnGly	J	5,414	5,476	63	ACT	GTA	UCC	0
nad3	J	5,477	5,830	354	ATA	TAG		0
trnAla	J	5,829	5,892	64	AGG	A	UGC	−2
trnArg	J	5,892	5,955	64	AAA	A	UCG	−1
trnAsn	J	5,956	6,018	63	TTA	AAA	GUU	0
trnSer	J	6,019	6,085	67	GGG	T	UCU	0
trnGlu	J	6,086	6,147	62	ATT	TA	UUC	0
trnPhe	N	6,146	6,207	62	ACT	TA	GAA	−2
nad5	N	6,208	7,921	1714	ATT	T		0
trnHis	N	7,922	7,984	63	ACT	GTA	GUG	0
nad4	N	7,986	9,321	1336	ATG	T		1
nad4l	N	9,315	9,602	288	ATG	TAG		−7
trnThr	J	9,605	9,666	62	GTT	CT	UGU	2
trnPro	N	9,667	9,728	62	CAA	GA	UGG	0
nad6	J	9,733	10,224	492	ATG	TAA		4
cytb	J	10,224	11,366	1143	ATG	TAG		−1
trnSer	J	11,365	11,429	65	AGT	TT	UGA	−2
nad1	N	11,455	12,405	951	ATT	TAG		25
trnLeu	N	12,406	12,468	63	ATT	ATA	UAG	0
16S ribosomal RNA	N	12,469	13,729	1261	GTT	T		0
trnVal	N	13,730	13,799	70	CAA	A	UAC	0
12S ribosomal RNA	N	13,800	14,544	745	AAA	A		0

Total length of F. binodulus mitochondrial tRNAs ranged from 61 bp (trnCys) to 71 bp (trnLys; Table 1). The length of nucleotide-amino acid accepter arms (AA arm) was uniformly 7 bp in all 22 tRNAs, while the anticodon arms (AC arm) varied from 3 to 5 bp. Compared with the AA and AC arms, length of dihydrouridine arms (D arm), and TΨC arms (T arm) were much more variable, with the length varying from 1 to 5 bp (Yang et al. 2018). All tRNAs could be shaped into typical cloverleaf structure except trnSer1, which had a shortened D arm (Fig. 2). This, however, is a well-known phenomenon for metazoan mitochondrial genomes, repeatedly reported in many animal species (Wolstenholme 1992).

Fig. 2.

Twenty-two tRNA structure originated from Figulus binodulus mitochondrial genome. Structure of 22 tRNAs with base pair of tRNAs. Names of tRNAs and anticodon were displayed bottom right of each structure. Mismatches are indicated in red bases, and wobble-pairs (G-U) are indicated as blue.

Twenty-eight mismatched base pairs were identified from the predicted tRNA secondary structures. G-U wobble pairs accounted for most of the mismatches (19 out of 28, 67.86%), which is a common feature in tRNAs (Varani and McClain 2000). UA-U mismatch base pair was found in the AA arm of trnIle, which were supported by the both predictions of ARWEN and MITOS (Fig. 2). In addition, two G-A, two C-A, and four U-U mismatches were also identified (Fig. 2).

Leucine (Leu) was the most frequently used amino acid in F. binodulus mitochondrial PCGs (Fig. 3), congruent to previous insect mitochondrial researches (Negrisolo et al. 2011, Wei et al. 2014, Zhang et al. 2014, Xin et al. 2017, Chen et al. 2018a). Serine (Ser) was the second most frequent amino acid and was followed by Isoleucine (Ile) and Phenylalanine (Phe; Fig. 3A). Interestingly, the most frequent two amino acids (Leu and Ser) were the only amino acids that had two different mitochondrial tRNA genes. In addition, these four amino acids showed different ratio of plus and minus strands (Fig. 3A): for example, Ile had the lowest proportion of minus strand and Leu and Ser are the next (Fig. 3A).

Fig. 3.

Codon usage of Figulus binodulus mitochondrial genome. (A) X-axis presents amino acids and Y-axis presents number of amino acids along with direction of genes (+ strand is blue color and – strand is red color). (B) Frequency of codons along with amino acids. Codons are displayed on the bars or around bars with arrows. Red colored codons indicate exceptional cases caused by RNA editing event.

Each codon in one amino acid displayed different proportions, indicating codon bias of each amino acid (Fig. 3B). In total, 63 codons including the exceptional start codon, ATT, were found in F. binodulus mitochondrial PCGs. The ratio of each codon varied however, usually with A-T-based codons out numbering the G-C based codons (e.g., TTT vs TTG in Phe; Fig. 3B), which is a feature steadily found in insect mitochondrial genomes (Dai et al. 2015).

Twelve out of 13 PCGs started with a typical start codon, ATN: 7 with ATG (nad2, atp6, cox3, nad4, nad4l, nad6, and cytb), 3 with ATT (atp8, nad5, and nad1), and 2 with ATA (cox2 and nad3). Cox1, on the other hand, had AAC as a start cordon which is known to be an alternative start codon in Polyphaga beetles (Sheffield et al. 2008). Three stop codons types were identified: the typical TAA and TAG stop codons and an abnormal T stop codon. Four genes (nad2, atp8, atp6, and nad6) had TAA codons, and other four genes (nad3, nad4l, cytb, and nad1) had TAG codons. The rest of the five genes (cox1, cox2, cox3, nad5, and nad4) ended with an incomplete T codon which is thought to be completed into a TAA stop codon by a posttranscriptional polyadenylation (Boore 1999, Meng et al. 2016).

Comparisons of F. binodulus Mitochondrial Genomes with 21 Complete Lucanid Mitochondrial Genomes

Including F. binodulus mitochondrial genome sequenced in this study, 21 complete mitochondrial genomes were available in family Lucanidae (Table 2). F. binodulus was third from the shortest, only longer than two species, Lucanus Mazama (NC_013578, 15,258 bp) and Dorcus hopei (MF612069, 16,026 bp; Fig. 4A). No correlation was found between the length of the mitochondrial genomes and the phylogenetic position of each species (see different bar colors in Fig. 4A) as the length of each sequence were mostly determined by species specific intergenic insertions (e.g., Prosopocoilus astacoides blanchardi: KF364622, 21,628 bp; Fig. 5B) or control region variations (e.g., Odontolabis curvera fallaciosa: MF908524, 19,614 bp).

Fig. 4.

Mitochondrial genome length, AT% vs AT-Skew, and GC% vs GC-Skew graphs in Lucanidae mitochondrial genomes. (A) X-axis indicates 21 Lucanidae mitochondrial genomes, and Y-axis means length of their complete mitochondrial genome. (B) X-axis present GC ratio and Y-axis show GC-skew. (C) X-axis present AT ratio and Y-axis show AT-skew. Orange dotted box indicates Figulus binodulus assembled in this study, and each color of symbol indicates subfamily

Fig. 5.

Nucleotide diversity and structural comparison of 21 Lucanidae mitochondrial genomes including Figulus binodulus. (A) X-axis indicates start position of sliding window (500 bp) used for calculating nucleotide diversity. Y-axis shows nucleotide diversity value. Blue line means nucleotide diversity of each window and yellow dotted line is average nucleotide diversity value. (B) Color arrow indicates PCGs (Orange), tRNAs (Pink), and rRNAs (Red). Orange dotted boxes indicate insertions.

Basal subfamilies, Syndesinae and Aesalinae, showed high AT ratio (from 73 to 76%; Table 2 and Fig. 4B) compared with that of Lucaninae except Ceruchus minor (Syndesinae; 68.39%); however, AT skew of the two subfamilies was similar to those of the Lucaninae mitochondrial genomes (Fig. 4B). This peculiarity of C. minor led to the dispersed positions of Syndesinae in both AT skew/AT ratio and GC skew/GC ratio graphs, where C. minor was clustered with the members of Lucaninae and Syndesinae (Fig. 4B and C). Excluding the exceptional C. minor, the three subfamilies could be roughly distinguished (Fig. 4B and C).

Nucleotide Diversity Analysis of F. binodulus Mitochondrial Genomes With 21 Complete Lucanid Mitochondrial Genomes

Based on multiple sequence alignments of available 21 stag beetle mitochondrial genomes (Table 2), nucleotide diversity was calculated. The total nucleotide diversity was 0.1331 (Fig. 5A), which was higher than those of family Ptinidae (0.08368; Park et al., under revision) and family Aphididae (0.0432; Park et al., in preparation). Cox1 presented the lowest nucleotide diversity among all PCGs (Fig. 5A). This phenomenon was also consistent with that of mitochondrial genomes of superfamily Scarabaeoidea including family Lucaenidae (Fig. 6). The nucleotide diversities of rRNA genes, region known to be well conserved in insect mitochondria (De Mandal et al. 2014) were very low (from 0.090 to 0.145; Fig. 5), were still higher than those of family Ptinidae (0.027 to 0.040; Park et al. under revision) and Aphididae (0.014 to 0.043; Park et al. in preparation). Nad2, nad6, atp8, and N-terminal of cytb showed nucleotide diversity ranges from 0.15 to 0.17 (Fig. 5A), higher than that of Ptinidae (from 0.05 to 0.13; Park et al., under revision). Nucleotide diversity of cox1 was flat from N-terminal to C-terminal (Fig. 5A), whereas atp6, nad5, and cytb showed sharp decrease of nucleotide diversity from N-terminal to C-terminal (Fig. 5A).

Fig. 6.

Nucleotide diversity of 13 PCGs in 21 Lucanidae complete mitochondrial genomes and 88 Scarabaeoidea mitochondrial genomes. X-axis indicates 13 PCGs and Y-axis indicates nucleotide diversity. Orange color indicates nucleotide diversities based on 21 Lucanidae mitochondrial genomes and Yellow color means nucleotide diversities from 88 Scarabaeoidea mitochondrial genomes.

Multiple insertions were found within the alignments (Fig. 5B), which resulted in drastic plummeting of nucleotide diversity (Fig. 5A). Most insertions were species specific (Fig. 5B); however, the insertion between trnGln and trnMet was shared in both species of subfamily Aesalinae (Fig. 5B). The AT-rich control region was the only region to have insertions while other insertions were also found in vicinity of the region (Fig. 5B). Even without the insertions, the control region showed extremely high nucleotide diversity (Fig. 5B).

All 13 PCGs displayed slightly different nucleotide diversities between Lucanidae and Scarabaeoidea groups (Fig. 6). While cox1, cox2, and nad2 were more conserved in Lucanid mitochondrial genomes than those of Scarabaeoidea, rest of the PCGs were more conserved in the superfamily level (Fig. 6). Atp8 and nad4l genes presented the largest differences between the two groups (Fig. 6), which were coincidentally the first and second shortest PCGs. Above half of the 88 available sequences (46 species) of Scarabaeoid mitochondrial genomes were from subfamily Scarabainae of family Scarabaidae, especially from genus Onthophagus (20 species). This unevenness in selected taxa may have reduced the diversities of the Scarabaeoidea group, thus the extreme conserveness of cox1 and cox2 of Lucanidae is a significant phenomenon. These diversity data of each gene will be useful in selecting proper molecular markers.

SSR Identified in F. binodulus Mitochondrial Genome

SSRs were rescued from F. binodulus mitochondrial genome sequences with the pipeline of the SSR database (see Materials and Methods). Among three types of SSRs, classified as SSRs, extended SSRs and potential SSRs, 4 SSRs (three monoSSRs, and one triSSRs), 2 extended SSRs (heptaSSR and octaSSR per each), and 44 potential SSRs were identified (Table 3). Along with unit sequence length, pentaSSRs occupied the most of all types of SSRs (32 SSRs; 64.00%; Fig. 7) and hexaSSRs covered the second largest proportion (12 SSRs; 32.00%; Fig. 7). Only one of each triSSR, heptaSSR, and octaSSR were identified (Fig. 7). Forty SSRs (80.00%) were in genic regions, whereas 10 SSRs (20.00%) were located in the control region.

Table 3.

List of 50 SSRs, extended SSRs, and potential SSRs from Figulus binodulus mitochondrial genomes

No	Name	SSR type	Type	Coordination		Unit sequence	Repeat number	Position	Genes
1	cH0000001	PotentialSSR	HexaSSR	329	340	AAAAAC	2	Genic	nad2
2	cP0000001	PotentialSSR	PentaSSR	781	790	ACTAT	2	Genic	nad2
3	c70000001	ExtenedSSR	HeptaSSR	999	1012	CCTTCTA	2	Genic	nad2
4	cP0000002	PotentialSSR	PentaSSR	2623	2632	AATTC	2	Genic	cox1
5	cP0000003	PotentialSSR	PentaSSR	6208	6217	ATAAA	2	Genic	nad5
6	cP0000004	PotentialSSR	PentaSSR	6638	6647	AAATA	2	Genic	nad5
7	cP0000005	PotentialSSR	PentaSSR	6698	6707	AAAAG	2	Genic	nad5
8	cP0000006	PotentialSSR	PentaSSR	6775	6784	AAATA	2	Genic	nad5
9	cP0000007	PotentialSSR	PentaSSR	6920	6929	CATTA	2	Genic	nad5
10	cP0000008	PotentialSSR	PentaSSR	7575	7584	CCATC	2	Genic	nad5
11	cP0000009	PotentialSSR	PentaSSR	7598	7607	AATTA	2	Genic	nad5
12	cP0000010	PotentialSSR	PentaSSR	7994	8003	AAACA	2	Genic	nad4
13	cP0000011	PotentialSSR	PentaSSR	8562	8571	CTAAT	2	Genic	nad4
14	cP0000012	PotentialSSR	PentaSSR	8770	8779	AAAAT	2	Genic	nad4
15	cH0000002	PotentialSSR	HexaSSR	9197	9208	CCAAAA	2	Genic	nad4
16	cP0000013	PotentialSSR	PentaSSR	9304	9313	AAATA	2	Genic	nad4
17	cH0000003	PotentialSSR	HexaSSR	9462	9473	ATACAA	2	Genic	nad4l
18	cH0000004	PotentialSSR	HexaSSR	9985	9996	TTAACC	2	Genic	nad6
19	cH0000005	PotentialSSR	HexaSSR	10289	10300	ACCTTC	2	Genic	cytb
20	cP0000014	PotentialSSR	PentaSSR	10529	10538	ATTAT	2	Genic	cytb
21	cP0000015	PotentialSSR	PentaSSR	11121	11130	CTTAT	2	Genic	cytb
22	cP0000016	PotentialSSR	PentaSSR	11218	11227	TTATT	2	Genic	cytb
23	cM0000001	SSR	MonoSSR	11471	11481	A	11	Genic	nad1
24	cH0000006	PotentialSSR	HexaSSR	11625	11636	AAAAGA	2	Genic	nad1
25	cH0000007	PotentialSSR	HexaSSR	11851	11862	AAACAT	2	Genic	nad1
26	cH0000008	PotentialSSR	HexaSSR	11934	11945	TTAAAG	2	Genic	nad1
27	cH0000009	PotentialSSR	HexaSSR	12006	12017	AATTAG	2	Genic	nad1
28	cP0000017	PotentialSSR	PentaSSR	12394	12403	AAATA	2	Genic	nad1
29	cP0000018	PotentialSSR	PentaSSR	12474	12483	TTTTC	2	Genic	16S rRNA
30	cP0000019	PotentialSSR	PentaSSR	12851	12860	TTAAA	2	Genic	16S rRNA
31	cP0000020	PotentialSSR	PentaSSR	12920	12929	AAAAT	2	Genic	16S rRNA
32	cP0000021	PotentialSSR	PentaSSR	13190	13199	ATTAC	2	Genic	16S rRNA
33	cP0000022	PotentialSSR	PentaSSR	13217	13226	TTAAT	2	Genic	16S rRNA
34	cH0000010	PotentialSSR	HexaSSR	13316	13327	AAAAAT	2	Genic	16S rRNA
35	cP0000023	PotentialSSR	PentaSSR	13421	13430	AATTA	2	Genic	16S rRNA
36	cP0000024	PotentialSSR	PentaSSR	13597	13606	AAATA	2	Genic	16S rRNA
37	cP0000025	PotentialSSR	PentaSSR	13947	13956	AATTA	2	Genic	12S rRNA
38	cP0000026	PotentialSSR	PentaSSR	14127	14136	ACAGG	2	Genic	12S rRNA
39	cP0000027	PotentialSSR	PentaSSR	14183	14192	AACTA	2	Genic	12S rRNA
40	cP0000028	PotentialSSR	PentaSSR	14280	14289	AATAA	2	Genic	12S rRNA
41	cP0000029	PotentialSSR	PentaSSR	14593	14602	TAGTT	2	Intergenic
42	cP0000030	PotentialSSR	PentaSSR	15021	15030	CTAAA	2	Intergenic
43	cP0000031	PotentialSSR	PentaSSR	15184	15193	CTAAA	2	Intergenic
44	cP0000032	PotentialSSR	PentaSSR	15466	15475	TTTAT	2	Intergenic
45	cM0000002	SSR	MonoSSR	15598	15612	T	15	Intergenic
46	cT0000001	SSR	TriSSR	15674	15691	TAT	6	Intergenic
47	cH0000011	PotentialSSR	HexaSSR	15856	15867	TATTTA	2	Intergenic
48	c80000001	ExtenedSSR	OctaSSR	15868	15883	AATAAATG	2	Intergenic
49	cH0000012	PotentialSSR	HexaSSR	16080	16091	AAATTA	2	Intergenic
50	cM0000003	SSR	MonoSSR	16168	16186	A	19	Intergenic

Fig. 7.

Distribution of number of SSRs along with unit length of SSRs identified on Figulus binodulus mitochondrial genome. X-axis presents SSR types and Y-axis indicates number of SSRs for each type. Numbers on bars means number of SSRs.

Four species of F. binodulus were selected to represent the neighboring lineages of F. binodulus: Ceruchus minor of tribe Ceruchini, and Odontolabis cuvera fallaciosa, Lucanus cervus, and Dorcus parallelipipedus of tribes Aegini, Lucanini, and Dorcini, respectively. SSRs were identified on the selected sequences using the same method (see Materials and Methods). Overall number of SSRs along with 10 types in the 5 species were similar to each other (Fig. 7): PentaSSRs and hexaSSRs explaining the first and second largest proportion of identified SSRs and extended SSRs, respectively. In detail, numbers of diSSRs, heptaSSRs, and octaSSRs in F. binodulus were the lowest among the five mitochondrial genomes (Fig. 7). Number of HexaSSRs showed interesting trend that C. minor was the largest and D. parallellpipedus is the smallest (Fig. 7), congruent to their phylogenetic relation (Fig. 8). In addition, mitochondrial genome of F. binodulus lacked diSSRs, whereas that of C. minor was the largest, which was same to that of L. cervus. It indicates that number of SSRs along with SSR types may not related to their phylogenetic position.

Fig. 8.

Phylogenetic trees of 21 complete mitochondrial genomes of Lucanidae. Neighbor joining (bootstrap repeat is 10,000) and maximum likelihood (bootstrap repeat is 1,000 phylogenetic trees of 21 Lucaninae, Aesalinae, and Syndesinae mitochondrial genomes with one Scarabaeidae complete mitochondrial genome as outgroup: Figulus binodulus (MN180052 in this study), Dorcus seguyi (MF612069), Dorcus seguyi (NC_038212), Dorcus hopei (MF612067), Dorcus parallelipipedus (KT876887), Dorcus hansi (NC_043928), Serrognathus platymelus (NC_044096), Prosopocoilus gracilis (NC_027580), Prosopocoilus Confucius (NC_036038), Prosopocoilus astacoides blanchardi (KF364622), Rhaetus westwoodii (MG159815), Lucanus Mazama (NC_013578), Lucanus cervus (NC_044476), Lucanus sp. (KT876903), Neolucanus maximus (NC_039652), Odontolabis cuvera fallaciosa (MF908524), Ceruchus minor (NC_042613), Sinodendron rugosum (NC_042614), Sinodendron yunnanense (NC_036157), Aesalus sp. (MH120282), Himaloaesalus gaoligongshanus (NC_042922), and Rhopaea magnicornis (NC_013252) as an out group. Phylogenetic tree was drawn based on maximum likelihood tree. The numbers above branches indicate bootstrap support values of maximum likelihood, neighbor joining, and Bayesian inference trees, respectively.

This trend of SSRs and extended SSRs was maintained outside of family Lucanidae. In S. paniceum (Coleoptera: Ptinidae), 21.11% of SSRs (Park et al., under revision) were located in the control region (Table 3), similar to the proportion found in F. binodulus control region of 10 out of 50 SSRs (20.00%), even though the total number of SSRs in F. binodulus was around half of that of S. paniceum (Park et al., under revision).

As SSRs are known to be prone to rapid evolutions (Stolle et al. 2013) as well as three mitochondrial genomes of S. paniceum, a cosmopolitan pest species, showed that 21.69% SSRs had variations among three individuals (Park et al., under revision), the SSRs identified from F. binodulus mitochondrial genomes could be used as molecular markers for distinguishing different populations or individuals enough.

Phylogenetic Interpretation of Lucanidae Mitochondrial Genomes

Thirty-seven genes on mitochondrial genomes covered >14 kb, which could be used for constructing high resolution phylogenetic trees (Simon et al. 1994, 2006; Boore and Brown 1998;Cameron 2014; Lavrov 2014; Smith and Keeling 2015; Yu and Liang 2018; Łukasik et al. 2019), displaying its usefulness to understand their taxonomic classification even though some exceptional cases were reported (Hwang et al. 2001, Park et al. 2020a).

We constructed phylogenetic trees based on 21 Lucanidae mitochondrial genomes covering three out of four subfamilies of Lucanidae except Lampriminae. Phylogenetic trees were constructed with three methods: maximum likelihood (ML), neighbor joining (NJ), and Bayesian inference (BI; see Materials and Methods). Inter-subfamily topology was congruent to the previously known phylogeny (Kim and Farrell 2015) in that Aesalinae is branched out first, then did tribe Sinodendrini and Ceruchini of Syndesinae, respectively (Fig. 8), showing a paraphyletic manner of subfamily Syndesinae. The support values of the phylogenetic trees, however, were low (Fig. 8); therefore, additional mitochondrial genomes from the minor subfamilies would be essential in future studies. Interestingly, the clade covering Ceruchini and all Lucaninae was well supported in all three trees (Fig. 8), which is congruent with that Cheruchus clustered with Lucaninae species not Sinodendron species in the AT ratio/AT skew and GC ratio/GC skew graphs (Fig. 4B and C).

Figulus binodulus was placed sister to all other Lucaninae mitochondrial genomes as previous studies on this subfamily were concentrated to several closely related species such as Dorcus spp. or Prosopocoilus spp. (Kim et al. 2015), not including a large proportion of the subfamily such as the Platycerini clade or the Gondwanan clade found in Kim and Farrell 2015. Lower classifications of the subfamily also turned out to be a mess: Prosopocoilus gracilis was shown to be more related to Dorcus spp. than other Prosopocoilus spp. (Fig. 8), which was a known phenomenon since the publication of the sequence (Wu et al. 2016). The recently revised genus Serrognathus was also not supported well as it did not form a separate clade (Fig. 8).