This study constructed phylogenetic trees using 25 species from the subfamily Iassinae as the ingroup, and Tinobregmus viridescens and Ponana quadralaba as the outgroup (these two groups were downloaded from NCBI; the accession numbers are SRR2496641 and SRR1821957). Specimen identification was completed by LJK through morphological characteristics. After the identification, the genital number was stored in a PCR tube containing glycerol, and the remaining tissues were used for total DNA extraction. Our voucher specimens were stored at the Institute of Entomology, Guizhou University, Guiyang, P. R. China.

Qubit precise quantitative detection and agarose electrophoresis were performed to assess the DNA quality. The concentration of our samples was between 28.83 ng/μL and 230.5 ng/μL by Qubit detection, and the total amount was between 1,614 ng and 13,138 ng. The target band of electrophoresis detection was located at about 10,000 bp, and all samples met the standard of resequencing and database construction. The total genomic DNA of 25 species of Iassinae was re-sequenced using the Illumina HiSeq 6000 sequencing platform (Beijing Berry and Kang Biotechnology) to obtain 150 bp double-end sequencing data. The clean data reported in this paper have been deposited in the Genome Sequence Archive (Chen et al. 2021) in National Genomics Data Center (CNCB-NGDC Members and Partners 2022), China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (GSA: CRA016669) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa.

Using the rapid-assembly pipeline PLWS v1.0.6 to assemble all sequenced genomes (https://github.com/xtmtd/PLWS; Zhang et al. 2019). We used the clumpify.sh script in BBTools-38.91 software for quality control to delete duplicate data (Bushnell 2014), and then the bbduk.sh component of the software for quality control to prune data smaller than Q15 and retain data larger than Q20. When the data were filtered and the length of the retained data was less than 15 bp, the entire sequence needed to be removed; it was also necessary to remove multiple N-repeated data. The bbnorm.sh component in BBTools v.38.91 software was used to reduce the abundance of data and reduce the average depth of high coverage data for downstream analysis (Chen et al. 2018).

The quality-controlled genomic data were assembled using Minia v.3.2.4 software to construct the initial contigs, and the contigs of each locus were grouped to reach the level of contings (Chikhi and Rizk 2013). Redundans v.0.13c software was then used to remove redundant contigs (Pryszcz and Gabaldon 2016). We used Minimap v.2.9 software to generate sequence-extended mapping files, and then SAMtoolsv.1.7 to convert the mapping files into sorted and indexed BAM format (Li et al. 2009; Li 2018). BESST.2.2 software and SOAPdenovo2 component in GapCloser.1.12 software were employed to extend contigs and fill gaps (Luo et al. 2012; Sahlin et al. 2014). Finally, we only kept sequences longer than 1,000 bp.

For matrix generation, we used BUSCO v3.0.2 (Waterhouse et al. 2018) to evaluate the integrity of the Benchmarking Universal Single-Copy Orthologs assembly and extracted USCOs from the genome. The hemiptera_odb10 dataset (n = 2,510) was used as the reference gene set, and the reference built-in species were Acromyrmex echinatior or Rhodnius prolixus. When the proportion of Fragmented USCOs in the evaluation was too large (>20%), then the average USCO length parameter in the control file “length_cutoff” could be modified to restore more “complete” USCOs, so that the BUSCO pipeline could treat more fragmented USCOs as “complete” USCOs sequences. For sequence alignment, amino acid sequences were aligned via MAFFT v7.450 with the L-INS-I method (Katoh et al. 2019). In order to remove unreliable homologous regions, we used trimAl v1.4.1 software for trimming. In continuation, we used the scripts component in PLWS for length filtering and retained the loci with parsimony information sites greater than 95. Based on the GC content strategy and relative composition variability (RCV) strategy, we performed compositional heterogeneity detection and filtering and retained the loci with GC content less than 0.13 and RCV value less than 0.5. Finally, based on the homogeneity detection strategy, the pval value of 0.05 was input after the symmetry test of the stationarity, reversibility and homogeneity (SRH). Based on PhyKIT v1.11.10 (Steenwyk et al. 2021), the filtered sequence was combined with the obtained data using the matrix generation.sh script as written by Zhang et al. (2019) to generate a data matrix USCO80 with 80% integrity (Table 2). For the remaining loci based on alignment marker filtering, IQ-TREE v2.1.3 was used to construct a single gene tree based on the EX-EHO mixed model (Hoang et al. 2018). The bootstrap was calculated, and the average bootstraps (ABS) strategy was selected to delete the loci with weak phylogenetic signals (ABS = 50). The sequences with integrity of 80% were combined to obtain USCO80_abs50 (Table 2).

The extraction of UCEs needed to be compared with the probe set of UCEs. For groups lacking probe sets, the probe design was required before the extraction of UCEs. The probe we used was a hemiptera UCEs probe sequence (Hemiptera-UCE-2.7K-v1) as developed by Faircloth (2017). Get UCE from assembled genomes using the UCE_extraction.sh script in the rapid-assembly pipeline PLWS. “phyluce _ probe _ run _ multiple _ lastzs _ sqlite”, was used to align the genome sequence with the probe sequence. “phyluce _ probe _ slice _ sequence _ from _ genomes” was used to extract the FASTA sequence matching the UCE marker in each genome. “phyluce _ assembly _ match _ contigs _ to _ probes” and “phyluce _ assembly _ get _ match _ count” were used to match the assembled contigs with the probes. “phyluce _ assembly _ get _ fastas _ from _ match _ counts” was used to extract all the successfully matched marker genes into a FASTA file. After extracting the UCE loci, the UCE matrix was constructed in a similar way to the USCO matrix. In short, by MAFFT alignment, trimAl was used to remove unreliable homologous regions, and Scripts components for length filtering to form a heterogeneity test (excluding loci with parsimony informative sites less than 100, GC > 0.5, RCV > 0.25). After the symmetry test of the SRH hypothesis, the pval value of 0.02 was input, and then the sequence was merged to generate a data matrix UCE80 with an integrity of 80%. IQ-TREE was used to construct the single gene tree of each locus based on the EX-EHO mixed model and the bootstrap was calculated. When the gene markers with high phylogenetic signal were obtained, the ABS value of 65 was input, and finally, the matrix UCE80_abs65 with integrity of 80% was generated for phylogenetic analysis. The details are shown in Table 2.

We used a different set of analysis methods for phylogenetic reconstruction to minimize the impact of systematic errors in large-scale genome datasets (Sun et al. 2020). We employed the maximum likelihood (ML) and Multispecies coalescent process (MSC) model methods to construct phylogenetic trees based on different data matrices. For partition ML analysis, we used ModelFinder (Kalyaanamoorthy et al. 2017) to select the best matching surrogate model for each gene partition through the relaxed clustering algorithm ‘-rclusterf 10’ (Lanfear et al. 2014) implemented in IQ-Tree. At the same time, to reduce the computational burden, we used LG ‘-mset LG’ and GTR ‘-mset GTR’ for USCO (amino acid) and UCE (nucleotide) partition ML analysis, respectively. To solve the influence of heterogeneity on phylogenetic inference (Kolaczkowski and Thornton 2004), the general heterogeneous evolution (GHOST) model was used on a single topology to construct the ML tree, in which the selected substitution model was the LG model, and the amino acid substitution frequency model and the rate heterogeneity model between sites were set to FO and H4 models, respectively. To reduce the effects of substitution saturation and compositional heterogeneity, we used Dayhoff6 recoding (6-state amino acid recoding strategies; Hernandez and Ryan 2021): Phylogears v2.2.0 was firstly used to convert the sequence format, and then a phylogeny was inferred with IQ-TREE with “-m GTR+R” model. The rate heterogeneity model between sites was set to the R4 model. All ML analyses ran 1000 SH-aLRTs (Guindon et al. 2010) and 1000 UFBoot2 (Hoang et al. 2018). A single gene tree as generated by UCE and USCO matrices in IQ-TREE, was used to perform MSC supertree analysis on incomplete lineage sorting (ILS) using default parameters (Zhang et al. 2018) in ASTRAL-III v5.6.1.

We use the MCMCtree in PAML v4.9j (Yang 2007) to estimate the divergence time, the Partition ML tree of the matrix USCO80_abs50 as the input tree, and the independent rate of the molecular clock. In order to avoid the falsely specified rate prior and unreasonable narrow posterior confidence interval due to the increase in the number of loci, the loci are merged into a larger partition based on the partition scheme of ML reconstruction estimation in the previous partition (Dos Reis et al. 2014; Jina and Brown 2018). The fossil records of Cicadellidae are extremely rare, and many available fossils are difficult to determine due to the lack of information on important morphological features (such as leg hairs and genitals) used for classification (Feng et al. 2024). Therefore, we calibrated the root age (113–125 million years ago, Mya) of the phylogenetic tree based on the oldest leafhopper fossil (Grimaldi 1990). Two independent runs were performed using the independent rate model, Sampfreq 5, until it collected 300,000 samples, while the first 250,000 samples were discarded as aging values.

We sequenced each sample and obtained raw data of 40 to 100 G, covering a range of 8.79 to 37.54 x. The assembled genome size ranged from 816.3 Mb (Batracomorphus curvatus) to 2,867.2 Mb (Krisna viridula), the number of Scaffolds from 213,833 to 857,534, the length of Scaffold N50 from 48.63 to 216.07 kb, the maximum read length from 59.31 to 639.44 kb, and the GC content (%) from 31.25 to 32.75. Additional statistical information (i.e., average read depth, number of stents, total length, maximum read length, and N50 length) is given in Table 1.

The overall completeness (complete and single-copy/ duplicated + fragmented) of the universal single-copy ortholog extraction was 33.50–70.30% (817–1,711 sites), of which 1.00–2.70% (25–68 sites) were duplicated, and 1.10–2.60% (27–65 sites) were fragmented (Fig. 1 and Table 1). The sequence of the same point was merged into the FASTA file to obtain 2404 USCOs, which were then used for downstream analysis. After sequence alignment, screening and pruning, 2,246 loci were retained. Length filtering and composition heterogeneity detection were performed on the pruned loci, and 1,599 loci were retained. A total of 1,586 loci were obtained by the SRH test, resulting in a USCO matrix, containing 368 loci (Table 2). The average bootstraps (ABS) strategy removed loci with weak phylogenetic signals, and retained 1432 loci. Subsequently, we used 1432 loci to infer a single gene tree and then selected 358 loci with an average bootstrap value greater than 80 to generate a matrix USCO80_abs50 for phylogenetic analysis (Table 2).

Figure 1. 

BUSCO genome completeness assessments: complete (C), complete single-copy (S), complete duplicated (D), fragmented (F), and missing (M).

After alignment and pruning, 2,416 UCE loci were retained, which were reduced to 1,901 loci by length filtering and composition heterogeneity detection. After SRH detection, it was reduced to 1,349 loci and finally reduced to 829 loci after the UCE matrix with 80% integrity (Table 2). Based on the average bootstraps (ABS) strategy, the sites with weak phylogenetic signals were deleted and reduced to 1,324 sites. Finally, 484 sites with an average value greater than 80 were selected to generate the UCE80_abs65 matrix for phylogenetic analysis (Table 2).

In this study, all phylogenetic reconstructions converged on a highly consistent topology. USCO and UCE matrices yielded 5 ML trees (Figs 2, S1–S4, respectively) and 2 species trees (Figs S5, S6, respectively). A comparative analysis of phylogenetic reconstruction shows differences in the topological structure of phylogenetic trees constructed using the same model based on different USCO matrices, while the UCE matrix phylogenetic reconstruction is concentrated on a highly consistent topology. All the results supported that the four genera involved in the study were restored to be composed of two main evolutionary branches, namely ((Batracomorphus + Trocnadella) + (Gessius + Krisna)).

Figure 2. 

ML phylogenomic tree of Iassinae based on the analysis of 358 USCO loci with the Partition model in IQ-TREE. Support values on nodes indicate SH-aLRT/UFBoot2, respectively. SH-aLRT/UFBoot2 approximate values are 100/100 node values. ML, maximum likelihood; USCO, universal single-copy ortholog.

In the branch (Batracomorphus + Trocnadella), Batracomorphus formed a large branch showing monophyly with species aggregation. This large branch divided the genus into three distinct clades with four species (Batracomorphus cornutus, Batracomorphus furcatus, Batracomorphus geminatus and Batracomorphus paradentatus), forming a small branch. In the genus Trocnadella, four species showed significant differences in their phylogenetic relationships as reconstructed by the two molecular markers. The phylogenetic reconstruction using the USCO matrix showed the developmental relationship as: (((Trocnadella arisana + Trocnadella testacea) + Trocnadella furculata) + Trocnadella fasciana))); however, the reconstruction using the UCE matrix showed the developmental relationship as: (((T. testacea + T. furculata) + T. arisana) + T. fasciana))). In the (Gessius + Krisna) clade, based on the phylogenetic relationship constructed by the USCO matrix, the species involved in the genus Krisna were not all clustered together. Among them, Krisna concava and Gessius rufidorsus clustered into a branch to become closely related species, and Krisna rufimarginata formed a separate branch. Krisna furcata and K. viridula always gathered together to form a sister group. In the phylogenetic relationship constructed by the UCE matrix, K. concava and K. rufimarginata were clustered together to form closely related species with G. rufidorsus, K. furcata and K. viridula formed a sister group.

The estimated divergence time (Fig. 3) demonstrated that the recent common ancestor divergence time of the four genera of Iassinae involved in the study was 119 Mya, and the 95% highest posterior density (HPD) interval was 113–125 Mya, with Early Cretaceous. In the Late Cretaceous-Oligocene (26–78 Mya) stage, four genera were differentiated. The differentiation of Batracomorphus, Trocnadella, Krisna and Gessius began at 60 Mya (49–72 Mya, 95% HPD), 44 Mya (33–57 Mya, 95% HPD), 64 Mya (51–78 Mya, 95% HPD) and 26 Mya (17–36 Mya, 95% HPD), respectively. Batracomorphus, Trocnadella, and Krisna originated in the Upper Cretaceous-Eocene, and Gessius originated in the Eocene-Miocene.

Figure 3. 

Divergence time estimation within Iassinae, inferred in MCMCtree on a Partition ML tree on matrix USCO80_abs50. Nodes are labeled with the age range of the most recent common ancestor, blue bar shows the 95% highest probability density interval. Q, Quaternary; Plio., Pliocene.