Mitogenomics of five Olidiana leafhoppers (Hemiptera: Cicadellidae: Coelidiinae) and their phylogenetic implications

Sample collection and DNA extraction

Details of sample collection are presented in Table S1. All specimens were preserved in absolute ethanol and stored at −20 °C until analysis. Genomic DNA was extracted from muscle tissues of adult males using the DNeasy® Tissue Kit (Qiagen, Germany). Total genomic DNA was eluted in 70 µL of double-distilled water. The remaining extraction steps were performed according to the manufacturer’s protocol. The obtained genomic DNA was stored at −20 °C until further analysis.

Mitogenome sequencing and assembly

Five Olidiana mitogenomes were sequenced using a next-generation sequencing platform (Illumina HiSeq 4000, Berry Genomic, Beijing, China; 6 GB raw data). Clean sequences were assembled using Geneious Primer version 2019.2.1 (Kearse et al., 2012), with O. ritcheriina (MK738125) (Wang et al., 2019c) as a reference.

Mitogenome annotation and sequence analysis

The locations of 13 protein-coding genes (PCGs) were identified using the ORF Finder tool of the National Center for Biotechnology Information and the invertebrate mitochondrial genetic code. Uncommon start and stop codons were identified by comparing our sequences with those of other Cicadellidae species. The locations and secondary structures of 22 transfer RNA (tRNA) genes were determined using tRNAscan-SE (Schattner, Brooks & Lowe, 2005) and ARWEN version 1.2 (Laslett & Canbäck, 2008). Ribosomal RNA (rRNA) genes were identified based on the loci of adjacent tRNA genes and compared with those of other Cicadellidae species (Wang, Li & Dai, 2017; Wang et al., 2019c). Repeat sequences within the control region were determined using the Tandem Repeats Finder tool (http://tandem.bu.edu/trf/trf.submit.%20options.html) (Benson, 1999). The annotated mitogenome sequences of the five Olidiana species have been deposited in GenBank with the accession numbers MN780581 –MN780585.

Base composition and relative synonymous codon usage (RSCU) in the mitogenomes were analyzed using MEGA version 6.06 (Tamura et al., 2013). Strand asymmetry was calculated using the following formulas: AT skew = (A − T)/(A + T); GC skew = (G − C)/(G + C) (Perna & Kocher, 1995). Intergenic spacers and overlapping regions between genes were manually counted.

Sequence alignment and phylogenetic analysis

To determine the phylogenetic relationships among members of Cicadellidae, 74 species from 12 subfamilies of Cicadellidae as well as 6 treehopper species were included, with two Cercopoidea species (Tettigades auropilosa and Cosmoscarta bispecularis) used as outgroups (Table S2). Phylogenetic analysis was performed by independently aligning the sequences of 13 PCGs and 2 rRNA genes. For each PCG sequence, terminal codons were removed before alignment using MAFFT version 7.0 in the Translator X online server (http://translatorx.co.uk/) with the L-INS-i strategy (Abascal, Zardoya & Telford, 2010; Castresana, 2000). Each rRNA gene was individually aligned using MAFFT with the G-INS-I strategy, and poorly aligned sites were removed using Gblocks 0.91b (Katoh, Rozewicki & Yamada, 2019). The resulting 15 alignments were concatenated using MEGA version 6.

Five datasets were concatenated for phylogenetic analysis: (1) PCGs, all codon positions of the 13 PCGs with 10,044 nucleotides; (2) PCG12, first and second codon positions of the 13 PCGs with 6,696 nucleotides; (3) AA, amino acid sequences of the 13 PCGs with 3,348 amino acids; (4) PCG12R, first and second codon positions of the 13 PCGs and 2 rRNA genes with 8,448 nucleotides; and (5) PCGR, all codon positions of the 13 PCGs and 2 rRNA genes with 11,796 nucleotides. The substitution saturation of four datasets (PCG, PCG12, PCG12R, and PCGR) was tested by plotting the number of transitions and transversions against genetic divergence using DAMBE (Xia, 2013).

A phylogenetic tree was constructed using the ML method with IQ-TREE (Nguyen et al., 2015; Zhang et al., 2019). The best-fit model was selected for each partition under corrected AIC using Partition Finder 2 (Table S3) (Lanfear et al., 2017) and evaluated using the ultrafast bootstrap approximation approach for 10,000 replicates. BI was performed using MrBayes version 3.2.6 (Suchard & Huelsenbeck, 2012; Zhang et al., 2019). Following the partition schemes suggested by Partition Finder, all model parameters were set as unlinked across partitions. Two independent runs with four simultaneous Markov chains (one cold and three incrementally heated at T = 0.2) were performed for 100 million generations, with sampling every 1,000 generations.

Mitogenomic characteristics of Olidiana species

The complete mitogenomes of five Olidiana species, namely, O. alata (MN780581; length 15,205 bp), O. longsticka (MN780582; length 15,993 bp), O. olbliquea (MN780583; length 15,312 bp), O. ritcheri (MN780584; length 15,372 bp), and O. tongmaiensis (MN780585; length 15,363 bp), were sequenced and assembled (Table S2). Their lengths were within the ranges of complete mitogenomes reported for other Cicadellidae species (14,805 bp for Nephotettix cincticeps and 16,811 bp for Idioscopus laurifoliae) (Song, Cai & Li, 2017; Wang et al., 2018). The mitogenomic architecture closely matched that of the inferred insect ancestral mitogenome (Crease, 1999): the newly sequenced mitogenomes had closed, circular DNA, typically comprising 37 genes (13 PCGs, 22 tRNAs, and 2 rRNAs) and a noncoding control region (Fig. 1). Of the 37 genes, most were encoded by the majority strand (J-strand) (9 PCGs and 14 tRNAs), whereas the minority strand (N-strand) encoded 14 genes (4 PCGs, 2 rRNAs, and 8 tRNAs) (Fig. 1, Table S4). However, the lengths of the 37 genes did not significantly differ between the five Olidiana species and other Cicadellidae species. The AT content of the five mitogenomes ranged from 78.0% (for O. alata) to 79.7% (for O. longsticka) and displayed a positive AT skew [0.147 (for O. longsticka) to 0.195 (for O. tongmaiensis)] and a negative GC skew [ −0.269 (for O. tongmaiensis) to −0.202 (for O. longsticka)] (Table 1). Additionally, the five Olidiana mitogenomes comprised 1–4-bp-long intergenic spacers at eight different loci, except for the trnY–COI intergenic spacer, which had 2–10-bp-long intergenic spacers. A total of 12 gene pairs were directly adjacent to each other, whereas the other gene pairs overlapped with each other, with overlap lengths of 1–4 bp, except for trnW–trnC and trnS2–ND1, which had large overlap lengths of 7–15 bp (Table S4). The gaps among the 37 genes in the mitogenomes were relatively smaller than those among genes in most reported Cicadellidae mitogenomes (Wang et al., 2018; Wang et al., 2019a; Wang et al., 2019b; Wang et al., 2019c; Wang et al., 2020a; Wang et al., 2020b).

PCGs and codon usage

Similar to that in other reported leafhopper mitogenomes (Wang et al., 2019b; Wang et al., 2019c; Wang et al., 2020a; Wu et al., 2016; Yang, Mao & Bennett, 2017), in the five Olidiana mitogenomes, the lengths of the 13 PCGs ranged from 150 bp (ATP8) to 1,674 bp (ND5) (Table S4). The AT content of the 13 PCGs ranged from 76.6% to 78.5%. Furthermore, the PCGs displayed positive AT skew (0.157–0.214) and negative GC skew (−0.299 to −0.234) (Table 1). Four PCGs (ND4, ND4L, ND5, and ND1) were coded by the N-strand, whereas the other nine (COI, COII, COIII, ATP8, ATP6, ND2, ND3, ND6, and CYTB) were coded by the J-strand. The Olidiana mitogenomes contained similar start and stop codons, and most PCGs had the typical start codon ATN (ATA/ATT/ATG/ATC) and either TAR (TAA/TAG) or an incomplete (single T) stop codon (Table S4). The presence of incomplete stop codons is a common feature of the mitochondrial genes among other leafhoppers, particularly of ATP8, and these incomplete stop codons are most likely caused by post-transcriptional modifications during mRNA maturation (Wang et al., 2018; Wang et al., 2019a; Wang et al., 2019b; Wang et al., 2019c; Wang et al., 2020a; Wang et al., 2020b; Yuan et al., 2019).

To understand the codon bias of the newly sequenced mitogenomes, RSCU and codon usage were determined. Codon usage was considerably similar between the five Olidiana mitogenomes and other Cicadellidae mitogenomes (Fig. 2). Among the five Olidiana mitogenomes, the most frequently used codon was UUA (for leucine; Leu). Leucine (Leu) 300–347, isoleucine (Ile) 350–404, methionine (Met) 353–375, and phenylalanine (Phe) 282–295 were the most frequently coded amino acids. However, five codons, including UCA, ACC, GUG, CCG, and GCG, were seldom used (Fig. 2). The codon usage pattern of Coelidiinae mitogenomes is highly consistent with that of previously sequenced Cicadellidae mitogenomes (Du, Dai & Dietrich, 2017a; Wang et al., 2020a; Wang et al., 2020b).

tRNAs and rRNAs

Consistent with most reported leafhopper mitogenomes (Li et al., 2017b; Wang, Li & Dai, 2017; Wang et al., 2018; Wang et al., 2019b; Wang et al., 2019c; Wang et al., 2020a; Wang et al., 2020b), the five Olidiana mitogenomes contained 22 tRNA genes, ranging from 57 bp (for trnC; O. longsticka) to 73 bp (for trnK; O. tongmaiensis) in length. The AT content of the tRNA genes ranged from 78.9% to 79.4%, displaying a positive AT skew (0.111–0.145) and negative GC skew (−0.150 to −0.110) (Table 1). In addition, all tRNAs exhibited a highly conserved, canonical cloverleaf secondary structure, except trnS1 (AGN), which lacked the stable dihydrouridine arm commonly found in most hemipterans (Fig. 3) (Cameron, 2014; Li et al., 2017b; Wang, Li & Dai, 2017; Wang et al., 2018; Wang et al., 2019b; Wang et al., 2019c; Wang et al., 2020a; Wang et al., 2020b).

Two rRNA genes (rrnL and rrnS) are highly conserved in Cicadellidae mitogenomes, and each of the five Olidiana mitogenomes contained these two rRNA genes. rrnL ranged from 1,176 bp (for O. ritcheri) to 1,186 bp (for O. longsticka) in length, whereas rrnS ranged from 729 bp (in O. alata) to 788 bp (in O. obliquea) in length (Fig. 1, Table 1). The rRNA genes of Olidiana mitogenomes displayed a positive AT skew (0.171–0.231) and negative GC skew (−0.307 to −0.284) (Table 1). rrnL was located between trnL2 and trnV, and rrnS was located between trnV and the control region (Table S4).

Control region

The control region of Olidiana mitogenomes ranged from 1,075 bp (for O. alata) to 1,804 bp (for O. longsticka). The differences in length in the control region were mainly attributed to the length and number of tandem repeats (R). All variable repeats in Olidiana mitogenomes were identified. Only a short unit (R) with two copies, both 115 bp in length, was present in O. alata. In O. longsticka and O. olbliquea, the first repeat region (R1) was 448 and 226 bp in length, respectively, and both comprised two units. The other two repeat regions, i.e., R2 and R3, were located after R1, and they were 345 and 433 bp (O. longsticka) and 129 and 28 bp (O. olbliquea) in length, respectively; both comprised three copies. O. ritcheri and O. tongmaiensis comprised four types of units: R1, 2 × 135 bp; R2, 2 × 195 bp; R3, 3 × 116 bp; and R4, 4 × 78 bp and R1, 2 × 281 bp; R2, 3 × 191 bp; R3, 3 × 81 bp; and R4, 3 × 4 bp, respectively (Fig. 4). Similar to the long intergenic spacers in other insect species, the repeat regions in leafhoppers may be attributed to an alternative origin of mitogenome replication (Dotson & Beard, 2001). The AT content (84.1%–85.8%) of the control regions was generally higher than that of the other regions. This is in part due to damage or accumulation of mutations in the mitochondrial DNA (Martin, 1995). The control regions of the five Olidiana mitogenomes displayed a slightly positive AT skew [ranging from 0.012 (for O. ritcheri) to 0.065 (for O. longsticka)] and negative GC skew [ranging from −0.231 (for O. ritcheri) to −0.061 (for O. alata)], except the control region of O. longsticka, which displayed a slightly positive GC skew (0.057) (Table 1). Moreover, these control regions were compared with previously reported control region sequences; their differences were very large, and no obvious correlation or similarity was found with existing sequences.

Phylogenetic relationships

No saturation was detected among the four candidate nucleotide sequence datasets (PCGs, PCG12, PCG12R, and PCGR) prepared for ML and BI analyses (all Iss <Iss.cSym or Iss. cAsym; P < 0.05) (Table S5), and the concatenated data were deemed suitable for phylogenetic analysis. Therefore, 10 phylogenetic trees were reconstructed based on five datasets (PCGs, PCG12, AA, PCG12R, and PCGR) using the BI and ML methods, and the main topological structures of the constructed phylogenetic trees were consistent (Figs. 5 and 6). Our results support the view that treehoppers originated from Cicadellidae and further confirm that Cicadellidae is a paraphyletic group, which is also supported by many previous studies (Du, Dai & Dietrich, 2017a; Du, Dietrich & Dai, 2019; Hu et al., 2019; Mao, Yang & Bennett, 2016; Song, Cai & Li, 2017; Wang et al., 2018; Wang et al., 2020a; Wang et al., 2020b; Yu et al., 2017). All analyses clearly support the monophyletic relationship of the 10 subfamilies within Cicadellidae and confirm that Iassinae and Coelidiinae as well as Megophthalminae and treehoppers are sister groups (Figs. 5 and 6, Figs. S1–S7). These results are consistent with those of previous studies (Wang, Li & Dai, 2017; Wang et al., 2018; Wang et al., 2019b; Wang et al., 2019c; Wang et al., 2020a; Wang et al., 2020b; Wu et al., 2016; Yang, Mao & Bennett, 2017).