Comparative analysis of the complete mitochondrial genomes of five Achilidae species (Hemiptera: Fulgoroidea) and other Fulgoroidea reveals conserved mitochondrial genome organization

Species collections and taxonomic identification

Adults of these five Achilidae insects were originally collected from Guizhou and Yunnan province in 2017. All samples were preserved in absolute ethanol immediately and stored at −80 °C in the Institute of Entomology, Guizhou University, Guiyang, China before morphological identification and DNA extraction.

The morphological identification of these five insects were performed by Shiyan Xu. Light microscope images were obtained by a Leica M125 stereomicroscope equipped with a Nikon D7000 digital camera, and were digitally processed with Helicon Focus and Adobe Photoshop CS6 software (Fig. 1). Sex of adult individuals was determined according to structures of the external genitalia. The morphological characteristics for comparison and the data in other species of Achilidae were taken from the previous studies (Chen, Yang & Wilson, 1989; Fennah, 1950, 1956; Long, Yang & Chen, 2015; Long et al., 2017; Xu, Long & Chen, 2018).

Nomenclatural acts

The electronic version of this article in portable document format will represent a published work according to the International Commission on Zoological Nomenclature (ICZN), and hence the new names contained in the electronic version are effectively published under that Code from the electronic edition alone. This published work and the nomenclatural acts it contains have been registered in ZooBank, the online registration system for the ICZN. The ZooBank Life Science Identifiers (LSIDs) can be resolved and the associated information viewed through any standard web browser by appending the LSID to the prefix http://zoobank.org/. The LSID for this publication is: urn:lsid:zoobank.org:pub:D1E6581B-15A0-4277-B00E-0564ABD04DE1. The online version of this work is archived and available from the following digital repositories: PeerJ, PubMed Central and CLOCKSS.

DNA extraction and sequencing

Total genomic DNA was extracted from each Achilidae species with TIANGEN Genomic DNA Extraction Kit (TIANGE, Beijing, China) according to the manufacturer’s instructions. Quality of the extracted DNA was checked on 1% agarose gel, sheared to 250–300 bp segments, A- and B- tailed and ligated to Illumina paired-end (PE) adapters. About 350 bp ligated fragments were quality selected on agarose gel and amplified to yield the corresponding short-insert libraries. Subsequently, HiSeq X 10 was used to sequence PE reads and the length of each read was 150 bp. Raw data was adapter clipped and qualified using tools fastx clipper and fastq_quality_filter in FASTAX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/index.html), respectively. The clean data was then used for mitochondrial genome reconstruction by SOAPdenovo (Luo et al., 2012) and GapFiller (Boetzer & Pirovano, 2012) with default parameters under the references of previously published mitochondrial genomes of related species. Finally, we gained complete mitogenomes of these five insects. All the detailed information is shown in Table S1.

Genome annotation and analysis

The initial annotation of the five mitogenomes, including gene prediction and non-coding RNA, were conducted using MITOS WebServer (http://mitos.bioinf.uni-leipzig.de/index.py) (Bernt et al., 2013). PCGs were translated into putative proteins on the base of the invertebrate genetic code using MEGA v5.1 (Tamura et al., 2011).

The identification and structure prediction of tRNAs were performed using MITOS and tRNAScan-SE Search Server available online (http://lowelab.ucsc.edu/tRNAscan-SE) with invertebrate codon predictors (Lowe & Eddy, 1997). While the undefined tRNAs were further compared through alignments with the nucleotide sequences of other species (Negrisolo, Babbucci & Patarnello, 2011; Zhu et al., 2017). The predicted secondary structures of all tRNAs were drawn by Adobe Illustrator CS6. Tandem repeats were identified by the tandem repeats finder online server (Benson, 1999).

The base composition of nucleotide sequences was described by skewness and was measured according to the following formulas (Perna & Kocher, 1995): AT-skew = (A − T)/(A + T) and GC-skew = (G − C)/(G + C). The values of nucleotide composition and relative synonymous codon usage (RSCU) were calculated using MEGA v7.0.26 (Tamura et al., 2011; Liu et al., 2014). The sequence data of the five insect mitogenomes have been deposited to GenBank under the accession numbers MH324927–MH324931 for Betatropis formosana (B. formosana), Magadhaideus luodiana sp. nov (M. luodiana), Peltatavertexalis horizontalis sp. nov (P. horizontalis), Plectoderini sp. (Pl. sp.) and Paracatonidia sp. (Pa. sp.), respectively.

Phylogenetic analyses

Besides five newly sequenced mitogenomes, 21 complete or nearly complete mitogenomes from Fulgoroidea insect species were used in phylogenetic analysis, with Pomponia linearis and Meimuna opalifera in Cicadoidea, and Cosmoscarta sp. and Aeneolamia contigua in Cercopoidea as outgroups (Liu et al., 2014; Song, Cai & Li, 2017). Their accession numbers and information are listed in Table S2. The nucleotide sequences of all PCGs were aligned based on amino acid sequences with Muscle criterion (Edgar, 2004) implemented in MEGA v7.0.26, respectively. The subsequent alignments were then concatenated by Bioedit v7.0.5.3 (Hall, 1999) and were used to reconstruct phylogenetic trees by ML and BI, respectively. The best substitution models and partition schemes for both ML and BI trees (Table S3) were estimated using PartitionFinder v2.1.1 (Lanfear et al., 2017), with the greedy algorithm. For ML analysis, IQ-Tree v1.4.3 (Nguyen et al., 2015) was used with 1,000 replicates of ultrafast likelihood bootstrap (Minh, Nquyen & Von Haeseler, 2013) to obtain bootstrap branch support values. Bayesian analysis was conducted with MrBayes on XSEDE v3.2.6 (Huelsenbeck & Ronquist, 2001) available from the CIPRES Science Gateway (https://www.phylo.org/) (Miller, Pfeiffer & Schwartz, 2010). Four simultaneous Markov chains ran for 20 million generations and sampled every 1,000 generations after discarding the first 25% “burn-in” trees. Nodal supports were assessed by the value of Bayesian posterior probabilities (BPP). The consensus trees were viewed and edited by Figtree v1.4.3.

Genome structure, organization and nucleotide composition

Each of all five newly sequenced mitogenomes was a circular double-stranded DNA molecule (Fig. 2), containing 37 genes (13 PCGs, two rRNAs and 22 tRNAs) and a non-coding control region, as commonly found among metazoans (Negrisolo, Babbucci & Patarnello, 2011; Liu et al., 2014; Wang & Tang, 2017). There were no differences in the gene order of all the five mitogenomes compared with the putative ancestral gene order (Table S4) (Dai et al., 2014, 2016). However, the five tRNA genes (trnW, trnC and trnY; trnT and trnP) order of Delphacidae mitogenomes differed from the putative ancestral gene order observed in Achilidae, Flatidae, Fulgoridae, Issidae and Ricaniidae mitogenomes (Fig. 3). In contrast to the traditional tRNA genes order, trnC and trnW have exchanged position reciprocally. A similar location exchange was also observed between trnT and trnP, and the location of them and nad6 also changed. A total of 23 genes were located on the heavy-strand (H-strand) with the remaining 14 genes on the light-strand (L-strand).

The sizes of the five Achilidae mitogenomes sequenced in this study ranged from 15,214 (Pa. sp.) to 16,216 bp (Pl. sp.) (Table S4), of which the length differences varied from 55 bp between Pl. sp. and B. formosana to 1,002 bp between Pl. sp. and Pa. sp. The differences were congruent with those of control regions.

The nucleotide composition of 26 complete or nearly complete mitogenomes in Fulgoroidea was investigated through the calculation of A + T content, AT-skew and GC-skew in percentages. The variation of AT% ranged from 74.3% to 77.8%, with the average value equal to 76.32% (Fig. 4A; Table S5). The overall A + T content of five newly sequenced mitogenomes ranged from 74.4% in M. luodiana to 77.7% in B. formosana. The high AT-skew values among analyzed mitogenomes (0.091–0.284) indicated the occurrence of more As than Ts, which was also observed in other examined Fulgoroidea mitogenomes (Fig. 4B; Table S5). Moreover, values of the AT-skew of these five mitogenomes were higher than those of Delphacidae insects, which was mainly due to a relatively high AT skewness of PCGs at all codon positions.

Protein-coding genes

Six kinds of triplet initiation codons (ATN, GTG and TTG) have been frequently reported in the mitogenomes of other insects. In the five newly sequenced mitogenomes, except for nad1 of Pl. sp., M. luodiana and Pa. sp. started with GTG and nad5 of M. luodiana started with GTG, the typical ATN codons were used as a start condon. The complete stop codon TAA was more frequently used than TAG, while incomplete termination codon TA/T was found in PCGs of five mitogenomes. This latter evidence well documented in other insect mitogenomes and with functionality of stop codons supposed to be restored through post-transcriptional polyadenylation (Boore, 2004; Cao et al., 2006; Ojala, Montoya & Attardi, 1981).

The average nucleotide composition analysis of the PCGs of Fulgoroidea mitogenomes showed that the PCGs in the H-strand had much higher AT-skew than those in the L-strand (Fig. 5). Moreover, the AT-usage differences between both strands were more significant in Achilidae, Fulgoridae, Flatidae, Issidae and Ricaniidae than in Delphacidae. The nucleotide composition of the PCGs of the five mitogenomes showed an extreme AT nucleotide bias at all codon positions (Fig. 6). Furthermore, the average nucleotide composition of the first and third codon positions of the H-strand contained more A than T, while the second had more T than A. However, in the L-strand, all the three positions had more T than A.

The amino acid composition analysis of the proteins (Fig. 7; Table S6) coded by the mitogenomes also showed that Phe, Ile, Met, Leu2 and Ser2 were the five most common amino acids, making up more than half of the amino acids in sum (Table S6). The same features were also found in the mitogenomes of five Callitettixini species (Hemiptera: Cercopidae) (Liu et al., 2014), while other analyses based on mitogenomes showed significant differences (Wang & Tang, 2017; Yu et al., 2017). Additionally, there were more Phe and Met but less Ile and Leu2 in Achilidae, Fulgoridae, Flatidae, Issidae and Ricaniidae than in Delphacidae. In other words, the amino acid composition of the analyzed mitogenomes have family-level characteristics (CDspT, codons per thousand codons: Phe > Ile > Met > Leu2 > Ser2 in Achilidae, Flatidae and Fulgoridae; Phe > Ile > Leu2 > Met > Ser2 in Ricaniidae and Issidae; Phe > Ile/Leu2 > Met/Ser2 in Delphacidae).

Theoretically, 20 amino acids in the mitochondria corresponding to 22 tRNAs and 62 codons, indicated that each tRNA must recognize at least two types of codon. However, RSCU indicated strong biases of codon usage and amino acid composition (Fig. S1; Table S7). The codon usage comparison of amino acids Leu indicated that the codon UUN was more frequently used than the codon CUN. Some codons (GCG, GGC, GTG and CGC) were rarely used in the PCGs, which was also found in other insects (Chai, Du & Zhai, 2012; Sun et al., 2009; Wang & Tang, 2017).

tRNAs and rRNAs

All the 22 tRNAs of the five mitogenomes were identified with the total length between 1,397 and 1,417 bp. A total of 14 tRNAs were encoded by the H-strand and the other (eight) by the L-strand. The length of tRNAs ranged from 55 to 72 bp. Except for the trnS1, which lacked dihydrouridine (DHU) arm and formed a simple loop, all the tRNAs could be folded into canonical clover-leaf secondary structure with an aminoacyl arm of seven base pairs and anticodon arm of five base pairs (Fig. 8). Moreover, an unpaired A nucleotide was found in the anticodon stem of trnS1, and similar pattern was also detected in the acceptor stem of trnR. The unusual phenomenon had been reported in other hemipterans (Li et al., 2016; Su et al., 2018; Yuan et al., 2015). The Watson–Crick base pairs (A–T and G–C) were frequently observed in stem regions, while G–U wobble and mismatched pairs can also be found.

The percentage of identical nucleotides (%INUC) in the five mtDNAs was calculated (Table S8). The trnK had the highest %INUC (81.7%), and four tRNAs (trnL2, trnG, trnM and trnN) on the H-strand and three tRNAs (trnP, trnQ and trnL1) on the L-strand showed a higher %INUC (>70%) in five mitogenomes, while the %INUC of trnS1 located on the H-strand was only 25%. Moreover, in these conserved regions, six and five G–U wobble base pairs were observed in acceptor and DHU arm, respectively; while three (one AC in trnL2, one GA in trnW and one AA in trnK) mismatched base pairs were found. The above results indicated a high-level conservation of tRNAs encoded on the H-strand.

The large ribosomal RNA (rrnL) and small ribosomal RNA (rrnS) genes were oriented on the L-strand of the five mitogenomes, with a high A + T content (>70%), a negative AT-skew and a positive GC-skew (Table S5).

Gene overlaps and non-coding A + T-rich Region

Different from the nuclear genome, there are short spaces or even overlaps between genes in mitochondrial genomes. The total length of overlapping regions in the five mitogenomes ranged from 42 bp in Pl. sp. mitogenome to 52 bp in the mitogenome of Pa. sp. The longest overlaps (20 bp) occurred between trnF and nad5 in Pa. sp. mitogenome. A seven-bp overlap (ATGATAA) between atp8 and atp6 found in other insect mitogenomes was also observed in the mitogenomes analyzed in this study (Fig. 9A). This fragment was thought to be translated as a bicstron (Stewart & Beckenbach, 2005). However, another seven-bp overlap between nad4l and nad4 generally found in other insect mitogenomes was not observed in this study. An eight-bp overlap (AAGCCTTA) between trnW and trnC was detected in five newly sequenced mitogenomes, which was absent in species of Delphacidae (Fig. 9B). In the mitogenomes of Pa. sp., the total lengths of intergenic spacers (114 bp) was much longer than those in the other four Achilidae species (Table S4).

Compared with the other regions, the A + T-rich regions of the five mitogenomes between the rrnS and trnM exhibited more variation in length, ranging from 908 bp in Pa. sp. to 1,850 bp in Pl. sp. (Table S4). It is obvious that the length differences among them were consistent with the total length differences of their mitogenomes. This region also harbored the highest A + T content (avg. 82.91%) in the mitogenomes, while AT skewness and GC skewness were significantly different, which was also observed in other insect mitogenomes from Fulgoroidea (Table S5).

Two typical repeating sequences scattered throughout the five entire mitogenomes, were identified. Two repeat sequences (135 × 5 and 21 × 13) were found in the A + T-rich region of Pl. sp. mitogenome, comprising half of the whole control region (Table S9). However, the lengths of the repeat sequences in other control regions decreased with the length of corresponding control regions. This feature indicated that characteristics of these regions in Achilidae were taxon-specific, and the different sizes or copy numbers of repeat units had some influences on the size of the control regions and also of the mitogenomes.

Phylogenetic analyses

Because of limited mitogenome sequences of Fulgoroidea, we included in the phylogenetic analyses only 26 species to understand the evolutionary relationships of Achilidae with the other Fulgoroidea families. Phylogenetic trees of ML and BI analyses were constructed based on 13 PCGs nucleotides sequences from 30 species (Table S2; Fig. 10). Overall, both trees generated the same topology, and BI analysis provided more resolution with strong supports than ML analysis. In both trees, Fulgoroidea is divided into two groups: Delphacidae and ((Ricaniidae, (Flatidae, Issidae)), (Fulgoridae, Achilidae)). Four families (Achilidae, Delphacidae, Fulgoridae and Ricaniidae) including more than two species formed separate branches in both BI and ML analyses (bootstrap support values (BS) ≥ 96, and BPP = 1.00).

In both ML and BI topologies, the basal position of Delphacidae was well supported with high node values (BPP/BS = 1/100), which was also supported by many previous studies based on morphological characters and molecular data (Muir, 1930; Asche, 1987; Bourgoin, Steffen-Campbell & Campbell, 1997; Yeh, Yang & Hui, 2005; Gnezdilov, 2016). In both topologies, relationships between Achilidae (Clade C) or Fulgoridae (Clade A) and other families (Clade B) in Fulgoroidea was consistent: a more ancient Achilidae than Fulgoridae got strong supports, and the placement of these two families was congruent with previous studies based on either morphological or molecular data (Bourgoin, Steffen-Campbell & Campbell, 1997; Yeh, Yang & Hui, 2005; Urban & Cryan, 2007). Three families (Flatidae, Issidae and Ricaniidae) formed consistently a clade with strong nodal supports (BPP/BS = 1/100), in concordance with studies based on morphological, mitochondrial and nuclear genes (Muir, 1930; Asche, 1987; Urban & Cryan, 2007).