Structural features and phylogenetic implications of Cicadellidae subfamily and two new mitogenomes leafhoppers

Complete mitochondrial genome sequences facilitate species identification and analyses of phylogenetic relationships. However, the available data are limited to the diverse and widespread insect family Cicadellidae. This study analyzes and summarizes the complete mitochondrial genome structure characteristics of 11 leafhopper subfamilies and two newly sequenced Typhlocybinae species, Empoascanara wengangensis and E. gracilis. Moreover, using 13PCGs and rRNA data to analyze the nucleotide diversity, evolution rate, and the phylogenetic relationship between the subfamilies of 56 species, verifying the taxonomic status analysis of E. wengangensis and E. gracilis. The analysis results show that the genome structures of the subfamilies and the newly sequenced two species are very similar, and the size of the CR region is significantly related to the repeat unit. However, in the entire AT-skews and CG-skews, the AT-skews of other subfamilies are all positive, and CG-skews are negative, while Empoascini of Typhlocybinae and Ledrinae are the opposite. Furthermore, among 13PCGs, the AT-skews of 13 species are all negative while CG-skews are positive, which from Empoascini in Typhlocybinae, Idiocerinae, Cicadellinae, Ledrinae, and Evacanthinae. Phylogenetic analysis shows that ML and PB analysis produce almost consistent topologies between different data sets and models, and some relationships are highly supported and remain unchanged. Mileewinae is a monophyletic group and is a sister group with Typhlocybinae, and the sister group of Evacanthinae is Ledrinae + Cicadellinae. Phylogenetic analysis grouped the two newly sequenced species with other species of Typhlocybinae, which was separated from other subfamilies, and all Erythroneurini insects gathered together. However, E. gracilis grouped into a single group, not grouped with species of the same genus (Empoascanara). This result does not match the traditional classification, and other nuclear genes or transcriptome genes may be needed to verify the result. Nucleotide diversity analysis shows that nad4 and nad5 may be evaluated as potential DNA markers defining the Cicadellidae insect species.


Introduction
The insect mitochondrial genome (mtDNA) is the only extranuclear genetic information carrier in insects. It is usually a closed double-stranded DNA molecule with a measured molecular weight of about 15-20 kb. Usually, it contains 37 genes, including 13 protein coding genes (PCGs), NADH dehydrogenase 1-6 and 4L (nad1-6 and nad4L), cytochrome c oxidase subunits 1-3 (cox1-3), ATPase subunit 6 and 8 (atp6 and atp8), cytochrome b (cytb), two ribosomal RNAs genes (16s and 12s) and 22 transfer RNA (tRNA) genes. A region rich in A + T, the control region, is also present [1,2]. Due to the characteristics of simple structure, small molecular weight, stable composition, conservative arrangement, lack of recombination, maternal inheritance, relatively high evolutionary rate and easy detection, the mitochondrial genome has been widely used for species identification and population genetic research as well as in biogeography and phylogeny [3][4][5]. However, there are repetitive regions and AT-rich regions in mitochondrial genes, and the real mitochondrial sequence and nuclear copy mitochondrial sequence (pseudogene) have great similarities, making it difficult to assemble the mitochondrial gene correctly after sequencing. How to enrich mitochondrial DNA more effectively is worthy of our thinking. In addition, in the application of mitochondrial genes, whether it is a phylogenetic tree constructed or the study of population evolution based on the mitochondrial genome, a single protein coding gene and rRNA gene are commonly used. Compared with the complete mitochondrial genome, A single gene fragment can only reflect part of the effective biological information, and different researchers often get different results based on different genetic data, resulting in the phylogenetic relationship of many species of insects still unresolved. With the continuous development of sequencing technology, it is necessary to use the whole mitochondrial genome as much as possible for phylogenetic analysis., In order to get more accurate results.
The hemipteran insect family Cicadellidae (leafhoppers) includes >2,600 genera and >22,000 species worldwide, including >2,000 species in China [6,7]. Erythroneurini, the larger tribe of the cicadellid subfamily Typhlocybinae, is widely distributed in the six major zoogeographic regions of the world and includes~2,000 species worldwide and >300 species in China [8,9]. All leafhoppers are phytophagous, different species feeding on a wide variety of plants, and the group includes critical agricultural pests and vectors of plant pathogens [10][11][12]. Simultaneously, because of its large number and small individuals, the taxonomic status and phylogenetic relationship between the subfamilies have been controversial, which has been discussed by related researchers. The study of the molecular phylogeny of Cicadellidae began in the 1990s. In 1993, Fang et al. sequenced and analyzed the 16S of 19 genera 21 species of Deltocephalinae, and made a reasonable attempt to study the molecular phylogeny of leafhopper insects [13]. Subsequently, in 1995, Fang et al. combined molecular data and morphological characteristics, and based on cytb to conducted a branch analysis of the Deltocephalinae genera in the New North Territory, and the results still verified the monophyletic of the group [14]. Dietrich et al. conducted a systematic study on the phylogeny of Cicadaceae. Since 1997, phylogenetic studies have been conducted on Flexamia, Dalbulus, Membracoidea, etc., based on mitochondrial and nuclear gene fragments, and the relationship between leafhopper subfamilies, tribes, and genera has been analyzed. Most of the results are similar to those based on morphological characteristics [15][16][17]. Hereafter, Dietrich combined morphology and molecule, divided Cicadellidae into 27 subfamilies, and revised the Cicadellidae classification system proposed by Oman [5,18]. However, these studies have not clearly established the relationship between the subfamilies in the Cicadellidae, and the relationship between some subfamilies and their relative groups remains to be explored [19,20].
Although in recent years, the phylogenetic relationship of the entire Cicadellidae and a few subfamilies within it has been studied based on morphological or molecular biological information [21][22][23][24], but the overall understanding of leafhopper phylogeny is preliminary. These recent studies have partially revised the Cicadellidae classification of high-level elements, but more data are still needed to reconstruct and verify its phylogenetic relationship. Therefore, in this study, we newly sequenced and annotated two species to increase the molecular data of the Cicadellidae, and combined with the mitochondrial gene data of 13 protein-coding genes and two ribosomal RNA genes of 56 Cicadellidae insects from 11 subfamilies (Table 1) to reconstruct the phylogenetic relationship between these subfamilies, and confirm the taxonomic status of Empoascanara wengangensis ) and E. gracilis (Dworakowska, 1992) at the molecular level. Besides, we analyzed the mitochondrial structure of these two species and each subfamily, including genome size and nucleotide composition, codon usage, tRNA secondary structure, A + T control region repeat unit, nucleotide diversity and evolution rate, and compared the similarities and differences between the various subfamilies. It is hoped that this study can provide a reference for future research on leafhopper classification and phylogeny.

Sample collection and DNA extraction
Samples of E. wengangensis and E. gracilis were collected from Duyun (107˚07 0 19 00 -107˚46 0 26 00 E, 25˚51 0 26 00 -26˚25 0 39 00 N) and Anshun (105˚22 0 50 00 -105˚45 0 22 00 E, 25˚33 0 38 00 -25˚55 0 32 00 N), Guizhou province, China, on 17 September 2018 and 13 May 2019. And leafhopper insects are not protected animals and are collected in non-natural reserves. The whole body was preserved in absolute ethanol and then stored at -20˚C in the laboratory. After morphological identification, voucher specimens with male genitalia prepared were deposited in the insect specimen room of Guizhou Normal University. Total DNA was extracted from the entire body without the abdomen.

Genome sequencing, assembly, and annotation
The mitochondrial gene sequence was obtained by sequencing. Primers were designed to amplify the mtDNA sequence in PCR reactions. The PCR reaction was performed using the LA Taq polymerase. The PCR conditions were as follows: initial denaturation 94˚C for 2 min, then 35 cycles of denaturation at 94˚C for 30 s, annealing at 55˚C for 30 s, and extension at 72˚C for 1 min/kb, followed by the final extension at 72˚C for 10 min. The PCR products were sequenced directly, or if needed first cloned into a pMD18-T vector (Takara, JAP) and then sequenced, by the dideoxynucleotide procedure, ABI 3730 automatic sequencer (Sanger sequencing) using the same set of primers. After quality-proofing of the obtained fragments, the complete mt genome sequence was assembled manually using DNAStar [57], and the Blast function in NCBI performed homology search to verify the amplified sequence as the target sequence [58,59]. The nucleotide base composition, codon usage and A + T content values were analyzed with MEGA 6.06 [60]. The secondary structure of tRNA genes was annotated using online tools tRNAscan-SE 1.21 [61] and ARWEN [62]. The tandem repeat sequence in the control area was determined by the online search tool Tandem Repeats Finder  Furthermore, the ratio between the non-synonymous (Ka) and the synonymous substitution rate (Ks) of 13 PCGs was also estimated in DnaSP 5.0.

Phylogenetic analysis
The phylogenetic analysis used the complete mitochondrial genomes of the two newly sequenced erythroneurine species plus 54 Cicadellidae species from our team and GenBank, and three outgroups of Cercopoidea (Table 1). The Gblocks Server online platform was used to eliminate poorly aligned positions and divergent regions of DNA protein alignment, and all alignments were checked and corrected in MEGA 6.06 [60] before phylogenetic analysis. Four datasets were generated: (1) 13 PCGs with 9262 nucleotides (PCGs); (2) the first and second codon positions of the 13 PCGs with 6174 nucleotides (PCG12); (3) 13 PCGs with 9262 nucleotides and two rRNA with 1219 nucleotides (PCGR); (4) and amino acid sequences of the 13 PCGs with 3289 amino acids (PCGAA). The trimmed alignment was used to estimate the phylogeny by Bayesian inference (BI) using MrBayes 3.2.7 [66] and maximum likelihood (ML) using IQ-TREE [67]. BI selected GTR + I + G as the optimal model, running 10 million generations twice, sampling once every 1000 generations, after the average standard deviation of the segmentation frequency drops below 0.01, with the first 25% of the samples are discarded burn-in, and the remaining trees used to generate a consensus tree and calculate the posterior probability (PP) of each branch. ML constructed with the IQ-TREE used an ultrafast bootstrap approximation approach with 10,000 replicates and calculate bootstrap scores for each node (BP).

Genome arrangement, organization and composition
In the Cicadellidae family, the structure and characteristics of the mitochondrial genomes of most leafhoppers are very similar, compared with the common gene rearrangements of Thysanoptera and Hymenoptera, the gene composition and arrangement of the mitochondrial genome of this family are relatively conservative. And the gene rearrangements are relatively rare, only three species have been reported from Deltocephalinae, and all rearrangements occur on the three genes trnW, trnC, and trnY in tRNA. Gene order of other species is the same as the putative ancestral insect (Drosophila yakuba) mitochondrial genome arrangement [1,40,42,43,68]. Among the 56 leafhoppers sequenced, the length of the mitochondrial genome ranges from 14 kbp to 17 kbp. The smallest is Krisna concava of Iassinae, while the largest is Eupteryx minuscula from Typhlocybinae, which is 14,304bp and 16,944bp in length, respectively. The average A+T content of the 11 subfamilies is 77.91% (A, 42.55%; T, 35.36%; C, 12.37%; G, 9.72%), of which Iassinae is the highest at 80.46%, while Ledrinae is the lowest at 75.94%. In all Cicadellidae species that have been sequenced, the content of base A > base T, base C > base G, except the six species from Ledrinae and Empoascini of Typhlocybinae, whose content is opposite. The range of AT-skews between subfamilies is -0.2246~0.1524, and the skew of GCskews is -0.2179~0.1230. All species are positive for AT-skews and negative for GC-skews, except six species from Empoascini in Typhlocybinae and Ledrinae. Iassinae has the highest AT-skews, and Coelidiinae has the lowest GC-skews (Table 2 and S1 Table).
Genome organization and nucleotide composition of the two new mitogenomes sequenced in this study are similar to those of other Erythroneurini reported previously [25][26][27]. The complete mitogenomes of E. wengangensis and E. gracilis are double-stranded plasmids with 14,830 and 14,627bp, respectively (Fig 1). Both contain the usual 13 PCGs, 22 tRNA genes, two rRNA genes, and a control region). Fourteen genes encode in the minority strand (L-strand) while the others encode in the majority strand (H-strand). E. wengangensis has a total of 45bp intergenic space in 12 regions ranging from 1 to 8bp. Eleven genes were found to overlap by a total of 47bp. E. gracilis has a total of 84bp intergenic space in 14 regions ranging from 2 to 15bp, and 11 genes were found to overlap by a total of 32bp ( Table 3).
The AT contents and skew statistics are shown in Table 4. The mitochondrial genomes of E.wengangensis and E. gracilis exhibit heavy AT nucleotide bias, with A + T% for the whole sequence 76.6% and 77.0%, respectively. Similar patterns of nucleotide composition are also found in other leafhopper species [38,47]. The control region (CR) has the strongest A + T% bias, while the PCGs shows the lowest A + T% among whole genes. The whole genome has positive AT-skews (0.015, 0.140) and negative GC-skews (-0.154, -0.157). Analysis of 37 individual genes of the two species shows that AT-skews are mostly positive, while for GC-skews, the genes of E. wengangensis are mostly negative, but E. gracilis are mostly positive (Fig 2 and  S2 Table). Positive AT-skews indicates that the content of base A is higher than that of base T. However, although the AT-skews is negative in a few genes, the difference in absolute value was minimal. For GC-skews, a negative value indicates that the content of base G is lower than that of base C, while a positive value is an opposite. Overall, the base composition of these two species is skewed toward A and C.  (Table 5 and S3 Table). All 62 available codons are used in 11 subfamilies. Synonymous codon usage bias was observed in 56 mitochondrial genomes, and UUA (Leu), UCG (Ser), AUU (Ile), AUA (Met), etc. are the most commonly used codons in many species. Located on the major strand (H-strand), while the other four PCGs are located on the minor strand (L-strand). The average AT content values of PCGs are 75.4% and 75.6% in E. wengangensis and E. gracilis, respectively, and the third codon position (80.8%, 80.9%) has an AT content much higher than that of the first (72.0%, 71.6%) and second (73.4%, 74.5%) positions. AT-skews of all codon positions are positive, while GC-skews are negative. All 13 PCGs have the standard ATN as the start codon, while nad5 and atp8 genes have TTG, a pattern also observed in other leafhopper mitogenomes [25,26]. Conventional stop codons (TAA or TAG) appear in 11 PCGs, except that cox2 and nad5 use an incomplete codon (a single T-) as the stop codon (Tables 3 and 4).

Protein-coding genes and codon usage
The relative synonymous codon usage (RSCU) was calculated and summarized, and All 62 available codons (excluding TAA and TAG) are used in E. wengangensis and E. gracilis. After excluding the stop codons, these two species have 3654 and 3658 PCG codons, respectively (Fig 3 and Table 6). Synonymous codon usage bias was observed in both mitochondrial genomes, and 22 codons are used more frequently than other codons. The four most abundant codons are UUA (Leu), UCG (Ser), GUA (Val) and GAA (Glu), and the least used codons were CGC (Arg) and ACG (Thr). The preferred codons all end with A or U, thus resulting in a strong A + T bias at the third codon position.

Transfer RNA and ribosomal RNA genes
The predicted tRNA length of 56 species is between 60~75bp, all of which are the typical clover-leaf secondary structure while only trnS1 lacks the dihydrouridine (DHU) stem and forms a simple loop. Furthermore, mismatches such as GU, AA, etc. often occur in tRNA. The structures and characteristics of 16S and 12S are similar in each subfamily. The sizes are 1038~1426bp, 707~789bp, with an average of 1188bp and 740bp, respectively. The AT content of 16S is higher than that of 12S. Deltocephalinae and Coelidiinae are the largest of 16S and 12S, respectively. Whole nucleotide compositions, AT-skew and GC-skew in 11 subfamilies of Cicadellidae. All 22 typical tRNA genes are present in the E. wengangensis and E. gracilis mitochondrial genomes, fourteen genes are oriented on the major strand (H-strand), whereas the others are transcribed on the minor strand (L-strand). Their nucleotide lengths are almost identical between species, ranging from 62 (trnG and trnT) bp to 71 bp (trnK and trnM) ( Table 3). The average AT content values of tRNAs are 77.5% and 78.7%, respectively, and the tRNA genes have negligible AT and GC-skews (Table 4).
Based on the secondary structure, a total of 20 and 21 G-U weak base pairs are found in E. wengangensis and E. gracilis of tRNAs respectively (Fig 4 and S1 Fig), forming weak bonds and located in AA stems (11bp), T stems (3 and 2bp) and DHU stems (6 and 8 bp). Most mismatched nucleotides are G-U pairs, which form weak bonds in tRNA and non-classical pairs in tRNA secondary structure, similar to other Cicadellidae [40,46].
Leafhopper ribosomal RNA (rRNA) includes 16S RNA and 12S RNA. These two genes are highly conserved and are encoded on the minor strand (L-strand). Similar to other known insects, the content of A + T% in 16S is higher than that of 12S. The 16s genes of E. wengangensis and E. gracilis are 1188bp and 1202bp in length, with AT content of 81.90% and 82.90%, respectively, and located between trnL2 and trnV. The 12S rRNA genes of both are 727bp and 735 bp in length, with AT contents of 79.80% and 80.30%, respectively, and located after trnV. The rRNA genes showed a positive AT-skew and GC-skew (Table 4). These features are similar to those observed in other insects [49,69,70].

PLOS ONE
Structural feature and phylogenetic implications of Cicadellidae subfamily and two new mitogenomes leafhoppers

Control region
The CR regions of 11 subfamilies range from 54 to 2662 bp, with an average of 1142 bp, the smallest in Iassinae and the largest in Typhlocybinae. The AT content is higher which ranging from 79.35% to 88.08%, with an average of 84.73%, and Mileewinae is the highest while Evacanthinae lowest. Among the 56 species, the repeat units range from 0 to 21. The smallest repeat unit from Typhlocybinae and Iassinae, and the largest repeat unit from Typhlocybinae. Spss 22.0 software was used to analyze the correlation between CR region size and repeat unit in 56 species. The results showed that there was a significant correlation between them (R = 0.672, P < 0.05), and the repeat units of most Cicadellidae insects were positively correlated with the size of CR regions (Fig 5). Like the typical insect mitochondrial genome, the mt genomes of E. wengangensis and E. gracilis have a sizeable non-coding region identified as the control region and located downstream of 12S. Control regions of both species are rich in AT, their lengths are 525bp and 212bp, and the AT contents are 83.1% and 91.1%, respectively ( Table 4). The control regions in the four available Empoascanara mitogenomes are various and not highly conserved, and their lengths range between 212 and 990 bp with variable numbers of repeat sequences (Fig 6). No tandem repeat units were found in E. gracilis; E. sipra includes one type of repeat unit (R); two kinds of repeats (R1, R2) are found in E. dwalata and E. wengangensis with various lengths and copy numbers.

Nucleotide diversity and evolutionary rate analysis
The sliding window analysis shows highly variable nucleotide diversity (Pi values) among 13 PCGs sequences of the 56 mitogenomes (Fig 7). The genes nad2, nad4, nad4L and nad5 have PCGs ranged from 0.251 to 0.655 (0 < ω < 1) (Fig 8), indicating that these genes are under purifying selection [71]. The genes nad4, nad5, nad1 and nad6 exhibit comparatively high Ka/ Ks ratios of 0.655, 0.626, 0.625 and 0.606, while the values of cox3, cox2, cytb and cox1 were relatively low, respectively 0.324, 0.323, 0.310 and 0.251. Nucleotide diversity analysis are primary for identifying the regions with large nucleotide divergence, especially useful for designing species-specific markers [72,73]. These are useful for taxa with highly variable morphological characteristics, especially Cicadellidae species. For a long time, gene cox1 has been considered as a universal barcode for identifying animal species [74], but in these 56 Cicadellidae species, it is the slowest evolving and least changing gene among 13 PCGs. If it is proved that the resolution of cox1 is very low in 13 PCGs, then other genes with sufficient large size, rapid evolution and high Ka/Ks ratio can be used as potential molecular markers in population genetics [73,75]. However, Part of the value of cox1 as a marker is also the availability of conserved primer sites and relatively low AT% bias. In this case, nad4 and nad5 may be evaluated as potential DNA markers that define the Cicadellidae insect species, but the result needs to be verified by multiple parties.

Phylogenetic relationships
Cicadellidae subfamilies' relationship has always attracted related scholars' attention because of the small size, and the morphological characteristics are difficult to distinguish. Early scholars made preliminary explorations of phylogenetic relationships based on morphological characteristics, and then analyzed them with molecular data. However, most studies are based on gene fragments, and data are relatively lacking, so the results do not reflect the relationship between subfamilies very well, and more data are needed for verification [24,28,46,76,77]. In this study, the results are the same as other recent studies of the mitochondrial genome (Figs 9 and 10). Most currently recognized leafhopper subfamilies were recovered as monophyletic but relationships among some subfamilies are not well resolved. And compared with the phylogenetic tree constructed using nuclear genes (28S), the relationship between the subfamilies is partly different, for example, Cicadellinae and Coelidiinae are sister groups, but they are not grouped together in this study. This may be because the study used only one nuclear gene fragment and the amount of data was insufficient [17]. At present, the most researched based on the complete sequence data are the horned leafhoppers of Cicadellidae, and other groups are rarely reported or not. Many studies have shown that Deltocephalinae species constituted one clade and tended to be placed at the tree's basal position as the sister group to the other leafhoppers, as in this study. The relationship between Mileewinae and its related subfamilies has  also been studied. Young transferred Mileewanini from Cicadellinae to Typhlocybinae in 1965, thinking that Mileewanini has more similarities with Typhlocybinae. Mahmood questioned this view and moved the tribe back to Cicadellinae in related work. Later [78], Young suggested Mileewanini as a separate subfamily in 1968 [79], and it has been adopted by many scholars in recent years [6,19,21]. In 2019, He Hongli used PCGs data to conduct a preliminary phylogenetic study on this subfamily's taxonomic status. The results show that Mileewinae is closer to Typhlocybinae than Cicadellinae and similar to Dietrich morphological phylogeny research [80]. Our research also shows that Mileewinae is a monophyletic group and is a sister group with Typhlocybinae. However, this result is only based on a species of Mileewinae, and more molecular data of this subfamily are needed to confirm this conclusion. Similarly, Evacanthinae's kinship issue has yet to be resolved. Although Dietrich and Wang Yang (2017) combined morphological and molecular data to conduct phylogenetic analysis, their relationship with closely related groups has not been clearly reconstructed [6,17,81]. The sister group of Evacanthinae may be Coelidiinae + Neocoelidiinae, Mileewinae + Typhlocybinae, Coelidiinae + Mileewinae + Signoretiinae, or Cicadellinae + Mileewinae. This paper shows that the sister group of Evacanthinae is Ledrinae + Cicadellinae. And phylogenetic analysis grouped the two newly sequenced species (E. wengangensis and E. gracilis) with other species of Typhlocybinae, which was separated from other subfamilies, and all Erythroneurini insects gathered together. However, E. gracilis grouped into a single group, not grouped with species of the same genus (Empoascanara). This result does not match the traditional classification, and other nuclear genes or transcriptome genes may be needed to verify the result. Based on the above research, although mitochondrial genes are now widely used in phylogenetic analysis, but the results obtained by different gene fragment combinations will have certain differences. The amount of data and the model have a great influence on the accuracy of the research results. Perhaps we need to combine other data to analyze the mitochondrial genome in a deeper level.

Conclusion
This paper analyzes and summarizes the complete mitochondrial genome structure characteristics of 11 leafhopper subfamilies and two newly sequenced Typhlocybinae species, E. wengangensis and E. gracilis, and analyzes the elemental composition, location, secondary structure and other characteristics of PCGs, tRNA genes, rRNA genes and control regions. Furthermore, using 13PCGs and rRNA data to analyzes the nucleotide diversity, evolution rate, and the phylogenetic relationship between the subfamilies of 56 species, verifying the taxonomic status analysis of E. wengangensis and E. gracilis. The analysis results show that the genome structures of the subfamilies and the newly sequenced two species are very similar, and the size of the CR region is significantly related to the repeat unit. However, in the entire AT-skew and CG-skew, the other subfamilies AT-skew are all positive, and CG-skew are negative, while Empoascini of Typhlocybinae and Ledrinae are the opposite. Moreover, among 13PCGs, the AT-skews of 13 species are all negative while CG-skews are positive, which from Empoascini in Typhlocybinae, Idiocerinae, Cicadellinae, Ledrinae, and Evacanthinae. This feature is consistent with the phylogenetic relationship displayed by the phylogenetic tree. We speculate that AT-skew and CG-skew have specific indications for the phylogeny of Cicadellidae species.
Phylogenetic analysis shows that ML and PB analysis produce almost consistent topologies between different data sets and models, and some relationships are highly supported and remain unchanged. Mileewinae is a monophyletic group and is a sister group with Typhlocybinae, and the sister group of Evacanthinae is Ledrinae + Cicadellinae. And phylogenetic analysis grouped the two newly sequenced species (E. wengangensis and E. gracilis) with other species of Typhlocybinae, which was separated from other subfamilies, and all Erythroneurini insects gathered together. However, E. gracilis grouped into a single group, not grouped with species of the same genus (Empoascanara). This result does not match the traditional classification, and other nuclear genes or transcriptome genes may be needed to verify the result. Nucleotide diversity analysis shows that nad4 and nad5 may be evaluated as potential DNA markers defining the Cicadellidae insect species. This study confirms the results of previous studies indicating that mitochondrial genome sequences are informative of leafhopper phylogeny. The new data provided here will facilitate future comparative studies of leafhopper mitogenomes and accentuate the need for more comparative data. However, more data is still needed to verify the above research results, especially for the subfamilies, whose taxonomic status is disputed. At present, there are little or no sequencing data.