Rhaetus westwoodi Parry (Coleoptera: Lucanidae) is the type species of monotypic genus Rhaetus Parry, which is distributed only in the region of the Eastern Himalayas. Little is known about the biology, ecology, genetics, and conservation status of this rare beetle. In this study, the mitochondrial genome of R. westwoodi was sequenced for the first time to obtain useful genetic data for this monotypic genus. Like most of known mitogenome in beetles, the mitochondrial genome of R. westwoodi is a closed circular molecule 18,134 bp in length, comprising 13 protein-coding genes, 22 transfer RNA (tRNAs) genes, 2 ribosomal RNAs (rRNAs), and 1 control region. The nucleotide composition is 36.30% adenine (A), 31.20% thymine (T), 21.49% cytosine (C), and 11.01% glycine (G), with a highly biased AT content of 67.50%. The phylogenetic analysis showed that R. westwoodi and Hexarthrius vitalis Tsukamotoi were recovered as one clade with strong nodal support (1.00 Bayesian Phylogenetics and Phylogeography [BPP], 100% maximum likelihood bootstrap [MLB]). The topology also strongly supported that R. westwoodi appeared to be a sister species to Pseudorhaetus sinicus Boileau (1.00 BPP, 100% MLB).
Rhaetus Parry is a monotypic genus with the type species of R. westwoodi Parry from Darjeeling, N. India (Parry 1862). As the monotype of Rhaetus, R. westwoodi is an amazing beetle due to its large size, shiny black body, and striking mandibles in males.
Nagai (2000) published the subspecies R. westwoodi kazumiae from North Myanmar, with the description of a large male. This male had a longer medial tooth, possibly distinguishing it from the typical R. westwoodi male. However, some large specimens of R. westwoodi from N. India also had the same mandibular characteristics of this reported subspecies and are depicted in the illustration of Didier and Séguy (1952). The mandibular “difference” of the subspecies, in our opinion, should be regarded as intraspecific variation among males of R. westwoodi. Therefore, in this study, we treated R. westwoodi as the only valid taxon in Rhaetus.
Little is known about the biology, ecology, and conservation status of R. westwoodi. The species is limited in its distribution to the high mountain forests (elevations averaging over 1,500 m above sea level) of the Eastern Himalayas (N. India, N. Myanmar, Tibet, N. W. Yunnan) (Arrow 1943, 1950; Benesh, 1960; Didier and Séguy 1952; Fujita 2010; Huang and Chen 2013; 2017; Krajcik 2001; Maes 1992; Parry 1862; Wan et al. 2010). Unfortunately, due to their fascinating appearance, specimens of this species have been much in demand as valuable collections, with an asking price of US$450 to 1,999 per large male (>80 mm in length) and US$90 to 350 per small male (approximately 60 mm in length). With this market, individuals are collecting the beetles for financial returns and, thus, may overhunt R. westwoodi in its restricted habitats to the point that the species may become increasingly endangered. In this study, we sequenced its complete mitogenome for the first time to provide genetic data for future conservation work and for establishing the phylogenetic relationships of R. westwoodi with other beetles.
Materials and Methods
Sample collection and DNA extraction
The voucher specimen of R. westwoodi was collected from Motuo, Tibet, in August 2006 by Ming Bai and deposited in the Museum of Anhui University (Hefei, China). Total genomic DNA was extracted from the muscle of R. westwoodi by using the Qiagen DNeasy Kit (Qiagen, Shanghai, China).
Polymerase chain reaction (PCR) amplification and sequencing
Complete mitogenomes were assembled from amplified fragments, with all primers used for amplification listed in Table 1. PCR amplification reactions were conducted in 25-μl volumes containing 10 μM of each primer (forward and reverse), 2 μl template DNA, 12.5 μl 2× EasyTaq SuperMix (+dye), and 8.5 μl sterile double-distilled water to a final volume of 25 μL. PCR amplifications were performed under the following conditions: an initial denaturation at 94°C for 2 min, followed by 35–37 cycles of denaturation at 94°C for 40 s, annealing at 52–58°C for 50 s, elongation at 70°C for 1 min, and then a final extension step at 72°C for 7 min. The temperature of annealing was determined by the length of the fragments. Sequencing was conducted with the Illumina HiSeq 2000 platform. Cluster strands created by bridge amplification were primed, and all four fluorescently labeled and 3-OH blocked nucleotides were added to the flow cell with DNA polymerase. The cluster strands were extended in single nucleotides. Following the incorporation step, the unused nucleotides and DNA polymerase molecules were rinsed away, a scan buffer was added to the flow cell, and then the optics system scanned each lane of the flow cell in imaging units (tiles). Once imaging was completed, chemicals that affect cleavage of the fluorescent labels and the 3-OH blocking groups were added to the flow cell, which prepares the cluster strands for another round of fluorescent nucleotide incorporation. The sequence was submitted to GenBank and assigned accession number MG159815.
Mitogenome assembly, annotation, and analysis
The mitogenomes were assembled using SOAPdenovo (BGI Company, Shenzhen, China), and preliminary annotations were made with the MITOS WebServer (http://mitos.bioinf.uni-leipzig.de/index.py). tRNA genes and their secondary structures were inferred using tRNAscan-SE version 1.21 (http://lowelab.ucsc.edu/tRNAscan-SE/). Those not identified by tRNAscan-SE, in addition to 16S rRNA (rrnL, lrRNA), and 12S rRNA (rrnS, srRNA), were determined according to sequence similarity with related species. The protein-coding genes (PCGs) were determined by ORF Finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html) under the invertebrate mitochondrial genetic code. Nucleotide compositions, codon usage, and relative synonymous codon usage (RSCU) values of PCGs were calculated with MEGA version 6.05 (Tamura et al. 2013). PCGs were translated with DNAMAN version 126.96.36.199 (Lynnon Biosoft, Vaudreuil-Dorion, Canada). Mitogenomes were mapped with CGView (Grant and Stothard 2008). Composition skew analysis was conducted according to the formulas AT skew = [A − T]/[A + T] and GC skew = [G − C]/[G + C] (Perna and Kocher 1995).
In this study, we retrieved 14 mitogenome sequences from GenBank (Table 2) and, to this data set, we added a newly sequenced R. westwoodi, generating a dataset of 15 taxa. Each of the PCGs was aligned individually, based on codon-based multiple alignments by using the MEGA version 7.0.26. Models of nucleotide substitution were selected according to the Akaike Information Criterion with jModelTest version 2.1.4 (Posada 2008). Phylogenetic trees were generated from maximum likelihood (ML) analysis using RAxML (Gillett et al. 2014) and Bayesian inference (BI) with MrBayes version 3.2.5 (Huelsenbeck and Ronquist 2013), both under the GTR + I + G model. Node support in the ML tree was estimated through bootstrap analysis with 1,000 replicates. The BI was conducted with two simultaneous Markov chain Monte Carlo runs of 2 million generations, sampled every 1,000 steps, with the first 25% discarded as burn-in. Phylogenetic trees were viewed and edited in Figtree version 1.4.3. (Rambaut 2016).
Results and Discussion
The complete mitogenome of R. westwoodi (GenBank accession number MG159815) measured 18,134 bp in length as a closed circular molecule (Fig. 1), consisting of 22 tRNA genes, two rRNA genes (rrn-L and rrn-S), 13 PCGs, and 1 control region. Four PCGs, 2 rRNAs, and 8 tRNAs were located on the N-strand, whereas the other 9 PCGs and 14 tRNAs were located on the J-strand (Table 3). Among the 13 PCGs, only 4 genes (nad4, nad4l, nad5, nad1) were located on the N-strand, whereas the other 9 (cox1, cox2, cox3, atp8, atp6, nad2, nad3, nad6, cytb) were located on the J-strand. These 13 PCGs comprised 3,705 codons, with the largest being the nad5 gene (1,806 bp) and the smallest the atp8 gene (155 bp). The average AT content of the PCGs was 65.65%, with the highest in nad4l (70.92%) and the lowest in cox1 (61.33%) (Table 4). All PCG codons started with ATN except cox1, which started with AAT (Table 3). Nine PCGs ended with TAA and TAG codons, with the incomplete stop codons TA and T in the remaining genes. In general, incomplete codon structures signal a halt of protein translation in insects and other invertebrates (Cheng et al. 2016, Masta and Boore 2004, Wu et al. 2014).
The RSCU analysis showed a biased usage of A and T nucleotides in the genome (Fig. 2). The most frequently encountered codons in other insects have been UUA (Leu), UUU (Phe), and AUU (Ile) (Wei et al. 2016). The PCGs coded on the majority strand show more C than G and on the minority strand show less C than G. This result agrees with the skew values of Wei et al. (2010).
All 22 tRNAs ranged from 61 to 71 bp (Table 3). The genes, determined using tRNA scan-SE, were typical metazoan insects mitogenomes (Cao and Du 2014, Sheffield et al. 2009, Wang et al. 2014, Zhang et al. 2014). The tRNA-S1 lacked the cloverleaf-shaped secondary structure owing to the structure of the dihydrouridine arm forming a structural loop (Cameron 2014a). The remaining 21 tRNAs possessed the typical cloverleaf structure.
The mitogenome sequence was found to be the same as that of other similar beetles (Kim et al. 2013, Lin et al. 2017, Liu et al., 2017, Sheffield et al. 2009, Wu et al. 2016). It has two ribosomal genes, the 1,268-bp-long rRNA-L and the 763-bp-long rRNA-S. The rRNA-L gene (AT content, 72.08%) was located between tRNA-L2 and tRNA-V, whereas the rRNA-S gene (AT content, 68.94%) was between tRNA-V and the AT-rich region.
The mitogenome was found to have a 7-bp conserved sequence (TTGTTCA) in the upstream of ND1, which has been observed in other beetles and insects (Cameron and Whiting 2008, Cameron 2014b, Kim et al. 2015, Lin et al. 2017, Sheffield et al. 2009, Yang et al. 2016). Intergenic spacers appeared frequently in the mitogenome, were interspersed throughout the PCGs and RNA genes, and ranged from 1 to 97 bp in length (Table 3).
The AT-rich region was 3,545 bp in length (70.87% AT) and was locked between rRNA-S and tRNA-M. A 16-bp poly-TA was locked in 17,079 bp downstream of the AT-rich region, whereas a 12-bp poly-T, 10-bp poly-A, and 21-bp poly-A were locked in 437 bp, 750 bp, and 696 bp upstream, respectively, of the tRNA-M.
The phylogenetic relationships of Scarabaeoidae were reconstructed based on the concatenated nucleotide sequences of 11 PCGs from 12 beetles as in-groups and other three scarabs as out-groups (14 of them downloaded from GenBank) (Table 2). BPP values ≥0.95 were considered as strong support, whereas BPP values ≤0.75 as weak support, and the ML bootstrap (MLB) values ≥75% as strong support, MLB values between 50% and 75% as moderate support, and ≤50% as weak support (Kim et al. 2015).
The topologies of the BI and ML analyses were highly consistent (Fig. 3). Rhaetus westwoodi and Hexarthrius vitalis Tsukamotoi were recovered as one clade with strong nodal value (1.00 BPP, 100% MLB). The two species appeared to be a sister group to Pseudorhaetus sinicus Boileau, with robust support (1.00 BPP, 100% MLB). The results corresponded to the morphological work of Huang and Chen (2013) and the latest molecular study inferred from the fragments (Wu and Wan 2016). Furthermore, these three species were in the tribe Dorcini first described by Parry (1864) with the same relationships we described herein. In contrast, the resolution did not agree with their relationships that were shown by Bartolozzi and Sprecher-Uebersax (2006), Benesh (1960), and Krajcik (2001), in which R. westwoodi and P. sinicus were in the tribe Rhaetulini and H. vitalis was in the tribe Lucanini.
We are grateful to Ming Bai (Institute of Zoology, Chinese Academy of Sciences) for donation specimens, Hu Li and Fan Song (Department of Entomology, China Agricultural University) for statistical advice, and laboratory colleagues for their kind help. This study was supported by the National Natural Science Foundation of China (no. 31201745, no. 31071954, and no. 31572311). Thanks for the funding provided by the National Natural Science Foundation.