The complete mitochondrial genome sequence of Oryctes rhinoceros (Coleoptera: Scarabaeidae) based on long-read nanopore sequencing

Sample collection, DNA extraction and ONT sequencing

An adult female O. rhinoceros was collected from a pheromone trap (Oryctalure, P046-Lure, ChemTica Internacional, S. A., Heredia Costa Rica) on Guadalcanal, Solomon Islands in January 2019 and preserved in 95% ethanol. Mrs Helen Tsatsia (Director of Research) and members of the research team at the Ministry of Agriculture and Livestock, Honiara, Solomon Islands Government facilitated the insect collection in Solomon Islands. Initially, the mitochondrial haplotype of the specimen was determined as CRB-G (Marshall et al., 2017) via Sanger sequencing of the partial cox1 gene sequence that was amplified using the universal barcode primers LCO1490 and HCO2198 (Folmer et al., 1994). High-molecular weight DNA was extracted using a customized magnetic (SPRI) bead-based protocol. Specifically, smaller pieces of tissue from four legs and thorax (50 mm³) were each incubated in a 1.7 ml eppendorf tube with 360 μL ATL buffer, 40 μL of proteinase K (Qiagen Blood and Tissue DNA extraction kit) for 3 h at RT, while rotating end-over-end at 1 rpm. A total of 400 μL of AL buffer was added and the reaction was incubated for 10 min, followed by adding 8 μL of RNase A and incubation for 5 minutes. Tissue debris was spun down quickly (1 min at 16,000 rcf) and 600 μL of homogenate was transferred to a fresh tube, where SPRI bead solution was added in 1:1 ratio and incubated for 30 min while rotating at end-over-end at 1 rpm. After two washes with 75% ethanol, DNA was eluted in 50 μL of TE buffer. DNA quality (integrity and concentration) was assessed on the 4,200 Tapestation system (Agilent, Santa Clara, CA, USA) and with the Qubit broad-range DNA kit. To enrich for DNA >10 kb, size selection was done using the Circulomics Short Read Eliminator XS kit. We sequenced a total of four libraries, each prepared with 1 μg of size-selected HMW DNA, following the manufacturer’s guidelines for the Ligation Sequencing Kit SQK-LSK109 (Oxford Nanopore Technologies, Cambridge, UK). Sequencing was done on the MinION sequencing device with the Flow Cell model R9.4.1 (Oxford Nanopore Technologies, Cambridge, UK) and the ONT MinKNOW Software.

An Illumina sequencing library was prepared using a NebNext Ultra DNA II Kit (New England Biolabs, Ipswich, MA, USA) and was sequenced on a HiSeq X10 (150 bp paired end reads) by Novogene (Beijing, China).

Mitogenome assembly, annotation and analysis

The Guppy base caller ONT v.3.2.4 was used for high-accuracy base calling on the raw sequence data, and only high-quality sequences with a Phred score >13 were used for the de novo mitogenome assembly with the program Flye v.2.5 (Kolmogorov et al., 2019) in the metagenome assembly mode. The method recovered the full circular assembly and to verify its accuracy, we first mapped the original reads back to the generated mitogenome assembly using Minimap2 (Li, 2018) with the following parameters: -k15—secondary = no -L -2. Second, we used BWA-MEM (Li, 2013) to map short-read Illumina sequences obtained from the whole-genome sequencing of another O. rhinoceros female collected from the same geographic location as the specimen used for the mitogenome assembly. The read alignment analysis in Pilon (Walker et al., 2014) was used to identify inconsistencies between the draft mitogenome assembly and the aligned short Illumina reads, removing small indels that represent homopolymers (e.g., >4 bp single nucleotide stretches) as an inherent sequencing error of the ONT (Mikheyev & Tin, 2014). Finally, we manually inspected if the Pilon correction occurred only in putative homopolymer regions by comparing the draft assembly with the Pilon-polished version.

The complete mitogenome sequence was initially annotated using the MITOS web server (Bernt et al., 2013), and tRNA genes and their secondary structures were cross-analysed using tRNAscan-SE v2.0 (Chan et al., 2019). To further refine the annotation and to examine mitogenome transcription, we used BWA-MEM to align Illumina reads from a transcriptome study of O. rhinoceros larvae (Shelomi, Lin & Liu, 2019), retrieved from the NCBI (SRR9208133). Finally, we manually inspected and compared our annotation to the complete and near complete mitogenome annotations of other related taxa (Table 1) in Geneious Prime (Biomatters development team, 2020). MEGA X (Kumar et al., 2018) was used to assess the codon usage and nucleotide composition of protein-coding genes. We used Geneious Prime to test if the nucleotide sequence of the cox1 gene recovers the recently reported PCR-RFLP marker (Marshall et al., 2017). This was done by aligning the sequences of the primer pair (LCO1490 and HCO2198) to isolate the amplicon fragment and perform in silico restriction digestion with MseI restriction enzyme. The restriction digestion of the amplicon produces a set of fragment lengths that distinguishes CRB-G from other haplotypes (Marshall et al., 2017). The presence of tandem repeats within the control region was assessed with the Tandem Repeats Finder v.4.0.9 (Benson, 1999) using default parameters. The annotated mitogenome sequence has been deposited in GenBank under accession number MT457815.1.

Phylogenetic analysis

To ascertain if our newly sequenced O. rhinoceros mitogenome can be correctly placed within the Dynastinae subfamily of the Scarabaeidae family, we performed the phylogenetic analyses with 15 additional taxa for which complete or near complete mitogenome sequences could be retrieved from the NCBI. We used thirteen species from five subfamilies of the Scarabaeidae family (Dynastinae, Rutelinae, Cetoniine, Melolonthinae, Scarabaeinae) and members of two other families from Scarabaeoidea as outgroups (Trogidae, Geotrupidae) (Table 1).

Nucleotide sequences of all 13 protein-coding genes (PCGs) were first translated into amino acid sequences under the invertebrate mitochondrial genetic code and aligned using the multiple alignment in Geneious Prime. The aligned amino acid matrix was back-translated into the corresponding nucleotide matrix and the Perl script Degen v1.4 (Zwick, Regier & Zwickl, 2012; Regier et al., 2010) was used to create the degenerated protein-coding sequences in order to reduce the bias effect of synonymous mutations on the phylogenetic analysis. These final alignments from all 13 PCGs were concatenated using Geneious Prime.

We estimated the phylogeny using two methods: the Maximum likelihood (ML) inference implemented in IQ-TREE web server (Trifinopoulos et al., 2016) and the Bayesian inference (BI) in MrBayes (Huelsenbeck & Ronquist, 2001). For the ML analysis, the automatic and FreeRate heterogeneity options were set under optimal evolutionary models, and the branch support values were calculated using the ultrafast bootstrap (Hoang et al., 2018) and the SH-aLRT branch test approximation (Shimodaira & Hasegawa, 1999) with 1,000 replicates. The Akaike information criterion (AIC) in ModelFinder (Kalyaanamoorthy et al., 2017) was used to select the best substitution model, and the BI phylogeny was generated with a total chain length of 1,100,000 (burn-in of 110,000 trees) and sampling every 200 cycles. The consensus trees with branch support were viewed and edited in Figtree v1.4.2 (Rambaut, 2014).

Mitogenome assembly, organization and transcription

The ONT long reads enabled the complete assembly of the circular mitochondrial genome for O. rhinoceros with just one MinION flow cell. This initial assembly was 20,954 bp in length and had a coverage of 3,834x. Based on the coverage ratio between the mitochondrial and draft nuclear genome assembly (unpublished data), we estimated that there were 320 copies of mitogenome for every nuclear genome copy. After including the data from three additional MinION flow cells, we produced a draft mitogenome assembly that was 21,039 bp long and had a coverage of 10,292x. Again, based on the coverage ratio between the mitochondrial and draft nuclear genome assembly, we estimated that there were 294 mitogenome copies per one copy of the nuclear genome in the total dataset. We then used the aligned Illumina reads to identify and correct indels, which are the most common error type produced by the ONT sequencing (Mikheyev & Tin, 2014), resulting in the removal of 141 bp. Of those, 95 bp were removed from the 5- to 9-mers (homopolymers). The final polished assembly was 20,898 bp in length. Our Illumina-only assembly, on the other hand, recovered a 17,665 bp mitogenome length with 3,233 bp missing from the control-region (CR).

The nucleotide composition of the complete mitogenome had high A + T bias (37.7% A, 32.8% T, 19.4% C and 10% G), and the long CR matched this genome-wide pattern (34.3% A, 35.8% T, 20.4% C and 9.5% G). The annotation revealed all 37 genes, with the rearrangement of trnI and trnQ genes, that showed the order: CR-trnQ-trnI-trnM-nad2 instead of CR-trnI-trnQ-trnM- nad2 (Table 2). CR resides within a large non-coding region (6,204 bp long) located between rrnS and trnQ, and the tandem repeats analysis revealed a complex structure of this large region with 11 putative repeats that had a consensus sequence between 7 and 410 bp repeated 2–12 times (Table 3). The length of 22 tRNAs ranged from 63 to 70 bp (Table 2), and their predicted secondary structures exhibited a typical clover-leaf structure. The length of rrnL and rrnS were 1,283 bp and 783 bp, respectively (Table 2).

All protein coding genes (PCGs) started with a standard initiation codon (ATN), 10 of 13 PCGs terminated with the conventional stop codons (TAG or TAA), while three genes (atp6, cox3 and nad5) had an incomplete stop codon T (Table 2). In silico digestion of the cox1 amplicon (delineated with the primer sequences from (Folmer et al., 1994)) produced the fragments 253 bp-, 138 bp-and 92 bp-long (Fig. S1) that are diagnostic for the CRB-G haplotype (Marshall et al., 2017).

The mapping of the transcriptome sequencing reads to the newly assembled mitogenome revealed that all PCGs were transcribed (mean coverage depth per base >23000X, Fig. 1), with cox1 and cox2 showing the highest level of expression when compared to the rRNA genes in the examined larval samples (Fig. 1). We found three domains within the large CR-containing region that also showed detectable transcription levels, with the transcript sizes of 137, 156 and 644 bp respectively. Our attempts to annotate these transcripts were not successful due to the fact that their sequences did not contain any open reading frames, nor did they have any significant BLAST hits within the NCBI’s reference RNAseq or nucleotide databases.

Phylogenetic analysis

The alignment of concatenated PCG sequences of 16 species (degenerated by Degen script (Zwick, 2013)) had 1,449 parsimony-informative sites, 1,818 singleton sites and 7,823 constant sites. ModelFinder identified GTR+F+R3 as the optimal substitution model (AIC: 94067.741, AICc: 94068.053), and given the proportion of invariable sites of 0.469 and the estimated gamma shape alpha of 0.728, we used the GTR+I+G model as the available alternative in MrBayes. This is the most parameter-rich model that has been shown to yield highly robust tree topologies for long alignments (>1,250 bp) in the smaller sample sizes (<17 taxa), irrespective of the model selection results (Abadi et al., 2019).

The topology of the maximum likelihood (IQ-tree) phylogeny was identical to the Bayesian (MrBayes) phylogeny, with O. rhinoceros being grouped with another member of the Dynastinae subfamily (Cyphonistes vallatus) with high confidence (100% SH-aLRT support, 100% ultrafast bootstrap support in ML, posterior probability 100% in BI). Dynastinae and Rutelinae were inferred as sister clades, that together with Cetoniine and Melolonthinae formed a basal split between phytophagus and coprophagus scarab beetles (Scarabaeinae) (Fig. 2).