Specimens of Evacanthus bivittatus were collected in Cangyuan County, Yunnan Province, China (23°16′39.67″N, 99°11′21.45″E). Specimens of Carinata ganga and Carinata recurvata were collected from the Lijiang River in Guilin City, Guangxi Zhuang Autonomous Region, China (25°36′28.74″N, 109°57′53.29″E). The specimens were deposited in the College of Life Science, Qufu Normal University, Shandong Province. All specimens were preserved in 100% ethanol and subsequently transferred to the laboratory for storage at –80°C. Adult specimens were identified based on their morphological characteristics. Genomic DNA was extracted from the head and thorax of specimens using the SanPrep DNA Gel Extraction Kit (Sangon, Shanghai, China).

Genomic DNA samples were sent to Personalbio Inc. (Shanghai, China) for library construction and next-generation sequencing. Libraries with an insert size of 400 bp were prepared from each DNA sample using the TruSeq™ DNA Sample Prep Kit (Illumina, USA), then sequenced as 150 bp paired-end reads on the Illumina NovaSeq platform. Species library reads have been deposited in BioProject (PRJNA1059151, PRJNA1081200, PRJNA1081201) (https://www.ncbi.nlm.nih.gov/bioproject/). Paired-end reads were de novo assembled using Novoplasty 3.7 (Dierckxsens et al., 2017). After trimming the adapters and removing short and low-quality reads, more than 4 GB (30–40 million reads) of clean data were obtained for each sample. Circular mitogenome sequences were extracted using DNASTAR 7.1 (Lasergene version 5.0) (Burland, 2000) with the complete mitogenome of Chudania hellerina (GenBank accession: MN227164) as an initial seed. Mitogenomes were then annotated using MITOZ v.1.04 (Meng et al., 2019) and manually checked using Geneious v.8.1.3 (Kearse et al., 2012). Secondary structures of tRNA were drawn manually using Adobe Illustrator CC2017, based on predictions from the MITOS Web Server (Bernt et al., 2013). Mitogenome maps were generated using the Organellar Genome DRAW (OGDRAW) program (Lohse et al., 2013). The complete mitogenome sequences of E. bivittatus, C. ganga, and C. recurvata were submitted to the NCBI sequence archive (GenBank accession no: NC077613.1, NC082036.1, OR979469.1, respectively).Various bioinformatic analyses, such as nucleotide composition, composition skew, codon usage of PCGs, relative synonymous codon usage (RSCU), and mitogenomic organization tables, were conducted using PhyloSuite v.1.2.2 (Zhang et al., 2020).

A total of 37 mitogenomes from the Cicadellidae were used in phylogenetic analyses, with two Fulgoridae family species serving as outgroups. Nucleotide sequences of the 13 PCGs, excluding stop codons, were aligned using the G-INS-i (accurate) strategy and codon alignment mode in MAFFT v.7 (Katoh et al., 2002). The rRNA gene sequences, rrnL and rrnS, were aligned using the Q-INS-I algorithm in MAFFT v.7, taking into account rRNA gene secondary structure.

Regions of ambiguous alignment were removed using Gblocks 0.91b (Talavera and Castresana, 2007). Gene alignments were concatenated using PhyloSuite v.1.2.2. For Maximum Likelihood (ML) and Bayesian Inference (BI) phylogenetic analyses, partitioning scheme and nucleotide substitution models were chosen using ModelFinder (Kalyaanamoorthy et al., 2017) with Bayesian Information Criterion. ML analyses were performed using IQ-TREE v.1.6.3 (von Haeseler et al., 2015) by applying the ultrafast bootstrap approximation approach with 5,000 replicates. BI analysis was conducted using MrBayes v.3.2.7a (Ronquist et al., 2012) in the CIPRES Science Gateway, using four chains (one cold and three hot chains). Two independent runs of 2,000,000 generations were carried out, with sampling every 1,000 generations. The first 25% of trees were removed as burn-in. On reaching an average standard deviation of split frequencies < 0.01, stationarity was assumed. The resulting phylogenetic trees were visualized using the online site: https://www.chiplot.online/ (Xie et al., 2023).

The mitochondrial genomes of E. bivittatus, C. ganga, and C. recurvata consisted of double-stranded circular DNA molecules, containing 37 classic mitochondrial genes, including 13 PCGs, 22 tRNA genes, and 2 rRNA genes, as well as a single control region. Among these genes, 23 were encoded on the heavy strand (H-strand), while the remaining 14 were transcribed from the light strand (L-strand) ( Figure 1 ).

The sizes of the E. bivittatus, C. ganga, and C. recurvata mitochondrial genomes were 15,133, 15,025, and 15,277 bp, respectively, attributable to differences in control region length (Fernández-Silva et al., 2003). The linear arrangements of the 37 genes in these three mitochondrial genomes were similar to those of other Cicadellidae species (Du et al., 2021; Wang et al., 2021a, 2021b). Excluding the control region, these mitochondrial genomes contained 9, 6, and 8 intergenic regions, respectively, as well as 9, 19, and 15 gene overlap regions. The mitochondrial genomes of all three species were highly condensed, with only a small amount of base overlap between neighboring genes, which suggests efficient RNA transcription and protein translation. The composition of the three mitochondrial genomes exhibited a strong bias toward A (adenine) + T (thymine), accounting for more than 70% of bases, consistent with the A-T bias observed in the majority of published Cicadellidae species (Wang et al., 2021a; 2019). This nucleotide composition asymmetry has often been considered an indicator of gene orientation and replication direction during replication and transcription processes (Perna and Kocher, 1995). The base skewness ([A-T]/[A+T]) values of E. bivittatus, C. ganga, and C. recurvata were 0.016, 0.057, and 0.055, respectively, with corresponding G-C skew values of –0.055, –0.105, and –0.12 ( Table 1 ).

The secondary structures of all tRNA molecules were identified using tRNAscan-SE 2.0 (Lowe and Eddy, 1997). The mitochondrial genomes of these three species each contained 22 tRNA genes, distributed throughout the entire mitochondrial genome, ranging in size from 61 to 73 nucleotides. Each mitochondrial genome had 14 and 8 tRNA genes encoded on the H- and L-strands, respectively ( Figure 1 ); this pattern is common in Cicadellidae mitochondrial genomes (Du et al., 2021; Wang et al., 2021a, 2021b). Additionally, two serine tRNA genes (tRNA^Ser(UCN) and tRNA^Ser(AGN)) and two leucine tRNA genes (tRNA^Leu(UUR) and tRNA^Leu(CUN)) were identified. The trnS1 gene lacked the dihydrouridine loop, preventing it from folding into the typical cloverleaf structure, while the other 22 tRNAs could fold into the typical cloverleaf structure. Furthermore, noncanonical base pairs or mismatched base pairs were common in the relatively conserved tRNAs. Secondary structure prediction identified 20, 15, and 21 noncanonical base pairs in E. bivittatus, C. ganga and C. recurvata, respectively. The three mitochondrial genomes had 30, 2, 4, and 17 cases of G-U, A-C, A-A, and U-U base pair mismatches, respectively. The majority of noncanonical nucleotides were G-U pairs, which are known to form weak bonds in tRNA and contribute to the formation of noncanonical base pairs in tRNA secondary structures ( Supplementary Figure S1 ). Furthermore, since the mitochondrial genome is not affected by recombination processes, these base pair mismatches may help eliminate detrimental mutations (Lynch, 1997).

The mitochondrial genomes of E. bivittatus, C. ganga, and C. recurvata all contained genes encoding the 12S and 16S rRNA subunits, which are located on the L-strand and separated by trnV ( Figure 1 ). The total lengths of the rRNAs in these three mitochondrial genomes differed, measuring 1905, 1886, and 1901 bp, respectively ( Supplementary Tables S1–S3 ). The A-T skew values of the two rRNAs in these leafhopper species were negative, while the G-C skew values were positive, indicating a higher AT content. The control region contains the initiation and control elements for mitochondrial genome DNA transcription and replication, and is the largest non-coding region in the mitochondrial genome (Taanman, 1999). The AT-rich region, which is the largest non-coding region in insect mitochondrial genomes, varies in size across diverse taxa, and even among individuals within the same species. Given the high variability of the insect AT control region, it is challenging to analyze its function and role in insect evolution through sequence similarity alignment. Through comparative analysis of insect AT control region sequences, Zhang et al. identified five conserved elements related to mitochondrial genome transcription and replication (Zhang et al., 1995). Taanman further proposed that the control region was the main regulatory region for mitochondrial genome replication and transcription (Taanman, 1999); however, in recent years, with the increasing numbers of sequenced and analyzed insect mitochondrial genomes, the AT control region has been found not to be a conserved structure in insects, and these regions do not contain conserved elements in many species (Dueñas et al., 2006).

The lengths of the 13 PCGs in E. bivittatus, C. ganga, and C. recurvata all totaled 10,929 bp, encoding 3,633, 3,634, and 3,634 amino acids, respectively. Among the 13 PCGs, 9 (ND2, COI, COII, ATP8, ATP6, COIII, ND3, ND6, and CYTB) were encoded on the H-strand, while the remaining 4 (ND5, ND4, ND4L, and ND1) were encoded on the L-strand. As in other leafhopper species, the ATP8 and ND5 genes in E. bivittatus and C. recurvata feature TTG start codons, whereas the other PCGs had traditional ATN (A/T/C/G) start codons. The stop codons of the 13 PCGs in these three species included TAA, TAG, and T-. In E. bivittatus and C. recurvata, incomplete stop codons (T-) were identified in the COI, COII, COIII, and ND1 genes. Furthermore, in E. bivittatus and C. recurvata, the CYTB gene harbored a TAG stop codon. The remaining genes concluded with TAA stop codons ( Supplementary Tables S1–S3 ). The existence of incomplete stop codons (T-) has been implicated in post-transcriptional mRNA regulation, potentially occurring subsequent to polyadenylation following transcription, a phenomenon commonly found in insect mitochondrial genomes (Ojala et al., 1981).

The RSCU was found to be very similar among the sequenced mitochondrial PCGs ( Figure 2 ); codons ending with A and T were more frequently utilized, with the three most prevalent codons being UUA (Leu2), UCU (Ser2), and UCA (Ser2). Collectively, our findings suggest that the length, orientation, and codon usage frequencies of PCGs in Evacanthini leafhopper species are highly conserved.

Next, we analyzed nucleotide polymorphism in the 13 PCGs using a sliding window approach. Among the 13 genes, COI exhibited significantly lower polymorphism ( Figure 3 ). Genetic distance analysis revealed that ATP8 (0.29) and ND2 (0.26) evolved relatively more rapidly, while COI (0.17) showed a slower rate of evolution ( Figure 4 ). Ka/Ks ratio was used to estimate gene evolution rate, which ranged from 0.1 (COI) to 0.44 (ATP8). These findings indicate that all the PCGs underwent purifying selection throughout the lengthy process of evolution. ATP8 exhibited relatively less stringent purifying selection, while COI experienced the most intense selection pressure ( Figure 4 ).

To elucidate the phylogenetic relationships within the Evacanthini of the Cicadellidae, a phylogenetic analysis was conducted on 37 Cicadellidae species, including E. bivittatus, C. ganga, and C. recurvata. Additionally, two species from the Fulgoroidea were selected as outgroups. Bayesian Inference (BI) ( Figure 5 ) and Maximum Likelihood (ML) ( Figure 6 ) phylogenetic trees were constructed using concatenated sequences of the 13 PCGs.

The resulting phylogenetic trees included 27 genera from 11 subfamilies of the Cicadellidae. They exhibited similar topologies and had high support values in both Maximum Likelihood bootstrap values and Bayesian posterior probabilities. The Evacanthini formed a distinct clade. Two species from the Carinata grouped together, as did two Evacanthus species. Two species formed a separate branch within the Nirvanini. The phylogenetic results also suggested that Evacanthini and Nirvanini are sister groups. Previous studies have supported the monophyly of Evacanthinae, which consists of four tribes: Evacanthini, Pagaroniini, Nirvanini, and Balbillini (Dietrich, 2004). Before the morphological phylogenetic analysis conducted by Dietrich (2004), Nirvaninae were regarded as a separate subfamily, while Evacanthini were included in Cicadellinae (Oman et al., 1990). Our phylogenetic studies corroborate Dietrich’ s conclusions, which are based on morphological traits related to the taxonomy of Evacanthinae and its tribes. These findings enhance understanding of the evolutionary relationships within this subfamily and further validate the use of morphological characteristics in insect taxonomy. The findings of this study also support the classification of Evacanthini and Nirvanini as two distinct branches within the Evacanthinae subfamily. A previous study reported a significant difference in the use of the anticodon (AGN) in trnS1 between representatives of Evacanthini and Nirvanini (Du et al., 2021); however, we found that both E. bivittatus and C. ganga used GCT as their anticodon, while C. recurvata had the TCT anticodon ( Supplementary Figure S1 ). Additionally, our phylogenetic analysis showed that C. recurvata belongs to the Evacanthini, suggesting that trnS1 (AGN) anticodon use may not be a major distinguishing factor between Evacanthini and Nirvanini. Our findings also reveal that the Evacanthini tribe is monophyletic.

Although the other subfamilies were clearly clustered, the phylogenetic interconnections among them were inconsistent. The Cicadellinae could be divided into two groups, with Mileewa species forming one group and species from other genera forming another. Previous research has not clarified the relationship between Evacanthinae and Cicadellinae (Dietrich, 2004, 2005; Wang et al., 2017). Further research and additional data are needed to determine the sister groups of Evacanthinae.

Overall, the topology of the phylogenetic trees was largely consistent with the evolutionary relationships among the leafhoppers. Our analysis confirmed the monophyly of the Evacanthini and provided insights into its taxonomic status within the Evacanthinae. Furthermore, we present comprehensive data from three Evacanthini mitochondrial genomes, contributing to clarification of phylogenetic relationships within the Evacanthini.