The highly rearranged mitochondrial genomes of three economically important scale insects and the mitochondrial phylogeny of Coccoidea (Hemiptera: Sternorrhyncha)

Sample preparation and DNA extraction

The specimens of U. yanonensis, P. citri and C. rubens were collected from Chengdu, Sichuan Province of China in October of 2019. The specimens were reliably identified by experts of Sichuan Academy of Agricultural Sciences, and were preserved in 100% ethanol. The total genomic DNA of the three scale insects was isolated using the E.Z.N.A.® Tissue DNA Kit (OMEGA, America) and preserved at −20 °C before the sequencing process.

Sequencing, assembly and annotation

The Illumina TruSeq short-insert libraries (insert size = 450 bp) were constructed using 1.0 µg of purified DNA fragments and sequenced by Illumina Hiseq 4000 (Shanghai BIOZERON Co., Ltd). Prior to assembly, raw reads were filtered and high-quality reads were retained and assembled into contigs by SOAPdenovo2.04 (Luo et al., 2012). Then the assembled contigs were aligned to the reference mitogenome of C. japonicus (GenBank accession number MK847519) using BLAST. The aligned contigs (≥80% similarity and query coverage) were arranged according to the reference mitogenome. Finally, the clean reads were mapped to the assembled draft mitogenome to fix the wrong bases; gaps were filled using GapFiller v2.1.1 (https://sourceforge.net/projects/gapfiller/). The mitogenome sequences of U. yanonensis, P. citri and C. rubens were deposited in GenBank under the accession numbers MT611525, MT611526 and MT677923, respectively.

Most tRNA genes were predicted and depicted by MITOS (Bernt et al., 2013); structures of several tRNA genes of C. rubens were predicted manually. PCGs and rRNA genes were identified by homology alignments. Gene boundaries of PCGs were confirmed in ORF finder (https://www.ncbi.nlm.nih.gov/orffinder/). The graphic view of the mitogenomes were computed using CGView Server (http://stothard.afns.ualberta.ca/cgview_server/) (Grant & Stothard, 2008). The probable mitochondrial rearrangement scenarios during the evolution of U. yanonensis, P. citri and C. rubens were predicted by the CREx (Common Interval Rearrangement Explorer) online server (Bernt, 2007) using Drosophila yakuba as a reference (Clary & Wolstenholme, 1985). Nucleotide composition of each gene and codon usage of PCGs were calculated by MEGA v.6.0 (Tamura et al., 2013). The composition skew analysis was conducted by AT-skew = [A−T]/[A +T] and GC-skew = [G−C]/[G +C] formulas (Perna & Kocher, 1995). The software DnaSP v. 5.10 (Librado & Rozas, 2009) was used to calculate the synonymous substitution rate (Ks) and the nonsynonymous substitution rate (Ka). Presumed secondary structures in the control region were predicted by the online tool Tandem Repeats Finder (http://tandem.bu.edu/trf/trf.advanced.submit.html) and DNAMAN v6.0.3.

Phylogenetic analysis

Nucleotide sequences of PCGs derived from four species of Coccoidea, including U. yanonensis, P. citri and C. rubens sequenced in this study, were used in the phylogenetic analysis (Table 1). The species S. coffeae was not included in the dataset due to the unannotated and unreliable status of its sequence as noted in Genbank. The two aphids, Aphis glycines and Diuraphis noxia were used as the outgroups. The 13 PCGs were aligned by MAFFT and concatenated as a combined dataset using SequenceMatrix v1.7.8 (Katoh & Standley, 2013). PartitionFinder v2.1.1 was used to determine the optimal nucleotide substitution models and partitioning schemes by using the Bayesian Information Criterion (BIC) and a greedy search algorithm (Lanfear et al., 2016). Two phylogenetic inferences were conducted with the partition schemes, including Bayesian inferences (BI) and Maximum likelihood (ML) analysis. BI analysis was conducted by MrBayes v3.2.7, with 10 million generations sampling every 1,000 generations, running one cold chain and three hot chains with a burn-in of 25% trees (Ronquist & Huelsenbeck, 2003). Stability of the results of BI analysis was examined by Tracer v.1.5. ML analysis was performed by RAxML v8.2.12 with 1,000 bootstrap replicates (Stamatakis, 2014). Tree files generated by both BI and ML trees were adjusted and visualized in FigTree v1.4.2.

Mitogenome annotation and nucleotide composition

The complete mitogenomes of U. yanonensis, P. citri and C. rubens were all typical double-strand circular molecules with a length of 15,220 bp, 15,549 bp and 15,387 bp, respectively (Fig. 1), which were similar to other mitogenomes of Coccoidea (Deng, Lu & Huang, 2019; Lu, Huang & Deng, 2020). The standard set of 37 genes (13 PCGs, 22 tRNA genes and two rRNA genes) were all found in the mitogenome of U. yanonensis (Table 2), whereas trnV was lost in P. citri (Table 3); C. rubens lacked five tRNA genes, trnC, trnR, trnS2, trnL1 and trnV (Table 4). In U. yanonensis, there were nine overlapping nucleotides located in four pairs of neighboring genes (Table 2); while in P. citri, there were 36 overlapping nucleotides in nine gene boundaries (Table 3). In C. rubens, there were only seven overlapping nucleotides in four gene boundaries (Table 4). The longest overlap was 18-bp long and located between trnS2 and ND1 in P. citri. There were 227 intergenic nucleotides (IGNs) dispersed in 20 locations for U. yanonensis, 126 IGNs in 19 locations for P. citri and 478 IGNs in 19 locations for C. rubens, indicating a loose structure of the three scale insect mitogenomes.

The whole mitogenomes of U. yanonensis, P. citri and C. rubens were strongly biased toward A and T nucleotides (86.6%, 82.7% and 87.5%, respectively). The U. yanonensis mitogenome had negative AT-skew and positive GC-skew, whereas P. citri and C. rubens exhibited positive AT-skew and negative GC-skew. The A+T contents were also rich in the mitochondrial genes, showing the highest in trnF of U. yanonensis and P. citri, and trnG of C. rubens.

Gene rearrangement

The mitochondrial genes of U. yanonensis, P. citri and C. rubens were highly rearranged, being different from the two sequenced scale insects, C. japonicus and S. coffeae (Deng, Lu & Huang, 2019; Lu, Huang & Deng, 2020). When compared with D. yakuba, U. yanonensis and P. citri both showed a conserved gene cluster trnE-trnF- ND5-trnH- ND4- ND4l-trnT-trnP- ND6-CYTB-trnS2- ND1-trnL1-rrnL; C. rubens had three shorter conserved gene clusters, COX1-trnL-COX2-trnK-trnD, COX3-trnG-ND3 and ND5-trnH-ND4-ND4L (Fig. 2). The mitogenome of U. yanonensis exhibited the rearrangement of three cytochrome c oxidase subunit genes (COX1, COX2, COX3), two NADH dehydrogenase subunit genes (ND2 and ND3) and many tRNA genes. Despite the multiple tRNA gene rearrangements, the mitogenome of P. citri also had a reversal of the ancestral gene cluster COX1- COX2-ATP8-ATP6- COX3- ND3. The mitogenome of C. rubens showed fewer rearrangements than U. yanonensis and P. citri, including two PCGs (ND2 and ATP8) and multiple tRNA genes.

The CREx analysis predicted the alternative scenarios how the three scale insect mitogenomes rearranged from the ancestral type of mitogenome of D. yakuba (Figs. 3–5). The mitochondrial gene order of U. yanonensis changed from D. yakuba by nine steps of rearrangement events, including the transposition of trnV and rrnS, the subsequent reverse transposition of trnK and trnD, the reversal of trnS1, and additional three reversal events and three tandem duplication and random loss (TDRL) events (Fig. 3). In P. citri, the first step is the reversal of trnK, followed by two alternative scenarios: the first one contained two reversal events, one TDRL event and one transposition event; the second one included three reversal events, two TDRL events and one transposition event (Fig. 4). Fewer rearrangement events were predicted in C. rubens, including the first step of transposition, the subsequent three reversal events, and final three TDRL events (Fig. 5). Considering the similarly rearranged mitochondrial genes of C. japonicus and S. coffeae, extensive mitochondrial rearrangement events are expected to occur very frequently in other unsequenced scale insects.

Protein-coding genes

The 13 PCGs of U. yanonensis were similar in size to those of P. citri, without truncated or duplicated PCGs (Tables 2 and 3). However, most PCGs of C. rubens were shorter than U. yanonensis and P. citri, especially for ATP8 and ND6 (Fig. 6). Most PCGs of the three mitogenomes utilized the standard ATN start codon (ATA, ATT, ATC and ATG). However, the special start codon TTG was used by COX1 of U. yanonensis (Table 2). Twelve PCGs of each mitogenome had the complete termination codon TAN (TAA or TAG), whereas ND6 of U. yanonensis and P. citri and COX2 of C. rubens ended with an incomplete stop codon T. In the previously sequenced scale insect, C. japonicus, COX2 also ended with an incomplete T (Deng, Lu & Huang, 2019). The relative synonymous codon usage (RSCU) values were calculated for the three mitogenomes (Fig. 7). In U. yanonensis, the most frequently used codon was TTA (Leu) whereas CTG(Leu), TCC(Ser), ACC(Thr), ACG(Thr), GCC(Ala), CAG(Gln), TGC(Cys), CGG(Arg) and AGC(Ser) were not used. In P. citri, the mostly used codon was also TTA (Leu), but CTC (Leu), AGC (Ser) and CGC (Arg) were the least. In C. rubens, TTA (Leu) was the most frequently used codon.

To evaluate the evolutionary rates of the PCGS, the average ratio of Ka/Ks was calculated for each PCG of the three mitogenomes (Fig. 8). The results showed that ND4L had the highest evolutionary rate, followed by ATP8 and ND5, while COX1 and CYTB appeared to be the lowest. The ratios of Ka/Ks were above 1 for most PCGs except for COX1 and CYTB, suggesting that these genes are evolving under positive selection. However, the ratios of Ka/Ks for COX1 and CYTB were below 1, indicating the purifying selection in these genes. The two genes, COX1 and CYTB which with relatively slow evolutionary rates have already been used as efficient phylogenetic markers in insects.

Transfer RNA genes

The typical set of 22 tRNA genes were all detected in the mitogenome of U. yanonensis, but trnV was absent from the mitogenome of P. citri (Figs. 9 and 10). In C. rubens, only 17 tRNA genes were recognized and the three tRNA genes trnA, trnQ and trnW were manually predicted (Fig. 11). Length and A+T content of the tRNA genes were subequal between U. yanonensis and P. citri, whereas the lengths of tRNA genes of C. rubens were generally shorter than U. yanonensis and P. citri. Individual tRNA gene of the three mitogenomes ranged in size from 49 to 75 bp; the longest tRNA gene was trnA in P. citri (Table 3); the shortest tRNA gene was trnY in C. rubens (Table 4). In the mitogenomes of U. yanonensis and P. citri, most of the tRNA genes could fold into cloverleaf secondary structures, but the dihydrouridine (DHU) arms of trnR and trnS1 were consistently lost. In C. rubens, most tRNA genes exhibited reduced DHU arms or T ψC arms. Such reductions of DHU arms were also reported in the tRNA genes of S. coffeae (Lu, Huang & Deng, 2020), suggesting that tRNA gene reduction could be a very common phenomenon in the mitogenomes of scale insects. The anticodons of the tRNA genes were identical among the three scale insects. In the tRNA genes of U. yanonensis and P. citri, a total of 12 and 19 mismatched base pairs were respectively identified and all of them were G-U pairs. In C. rubens, only four mismatched G-U pairs were identified.

Ribosomal RNA genes

There were two rRNA genes identified in in each mitogenome. The length and A+T content of each rRNA gene were subequal between U. yanonensis and P. citri, but the lengths of rRNA genes were much shorter in C. rubens (Tables 2–4). In U. yanonensis, the large ribosomal RNA (rrnL) gene was 1314 bp with an A+T content of 89.6%; the small ribosomal RNA (rrnS) gene was 807 bp with a high A+T content of 90.6%. In P. citri, the rrnL gene was 1396 bp with an A+T content of 86.9%; the rrnS gene was 856 bp with an A+T content of 88.4%. In C. rubens, the rrnL gene was 1,263 bp with a high A+T content of 90.7%; the rrnS gene was 587 bp with an A+T content of 87.9%. Locations of the two rRNA genes were similar to D. yakuba, being neighbored with the CYTB-ND1 PCG cluster (Fig. 2). Instead of the commonly found trnV between the rrnL and rrnS genes in other insects, the intermediate tRNA gene between the two rRNA genes was trnA in U. yanonensis, trnA and trnQ in P. citri, and completely absent in C. rubens.

Control region

Control region (CR), also known as A+T rich region, was the longest and most variable non-coding area in the three mitogenomes (Fig. 12). The CR of U. yanonensis was short with only 260 bp, being located between trnY and ATP8 and with a relatively high A+T content of 81.9% (Table 2). The CR of P. citri was much longer than U. yanonensis (678 bp), being located between trnM and trnI and with an A+T content of 84.4% (Table 3). Two putative CRs were found in the mitogenome of C. rubens: the 830-bp long CR1 between rrnS and trnF and the 800-bp long CR2 between trnF and trnM (Table 4). A+T content of the two CRs was 85.4% and 88.4%, respectively, higher than U. yanonensis and P. citri. The CR of C. japonicus and S. coffeae was 507 bp and 1454 bp, respectively, indicating the highly variable length of CRs in scale insects (Deng, Lu & Huang, 2019; Lu, Huang & Deng, 2020).

The CR of U. yanonensis was composed of 2.9 copies of tandem repeats; the first two copies had a consensus size of 91 bp, whereas the third repeat was 78 bp in length. The CR of P. citri contained three types of secondary structures that might function in regulating the replication and transcription of the mitogenome, including 2.3 copies of 110-bp long tandem repeats, one 40-bp long poly-[TA]n stretch, and a 21-bp long stem-loop (SL) structure. The SL structure was initiated by a “TAA” motif and ended with a “GTA” motif. The longer tandem repeats and extra secondary structures of P. citri resulted in the longer CR than that of U. yanonensis. The CR1 of C. rubens contained 3.6 copies of 33-bp long tandem repeats but had no SL structures. The CR2 of C. rubens included 5 copies of 24-bp long tandem repeats and a combined SL structure. The length, nucleotide composition, number and types of structural elements in CRs of the three mitogenomes were found highly variable, which implied that the scale insect mitogenomes were likely to be regulated in different ways during the mitogenomic replication and transcription processes.