Mitogenomic phylogeny of Typhlocybinae (Hemiptera: Cicadellidae) reveals homoplasy in tribal diagnostic morphological traits

Abstract The subfamily Typhlocybinae is a ubiquitous, highly diverse group of mostly tiny, delicate leafhoppers. The tribal classification has long been controversial and phylogenetic methods have only recently begun to test the phylogenetic status and relationships of tribes. To shed light on the evolution of Typhlocybinae, we performed phylogenetic analyses based on 28 newly sequenced and 19 previously sequenced mitochondrial genomes representing all currently recognized tribes. The results support the monophyly of the subfamily and its sister‐group relationship to Mileewinae. The tribe Zyginellini is polyphyletic with some included genera derived independently within Typhlocybini. Ancestral character state reconstruction suggests that some morphological characters traditionally considered important for diagnosing tribes (presence/absence of ocelli, development of hind wing submarginal vein) are homoplastic. Divergence time estimates indicate that the subfamily arose during the Middle Cretaceous and that the extant tribes arose during the Late Cretaceous. Phylogenetic results support establishment of a new genus, Subtilissimia Yan & Yang gen. nov., with two new species, Subtilissimia fulva Yan & Yang sp. nov. and Subtilissimia pellicula Yan & Yang sp. nov.; but indicate that two previously recognized species of Farynala distinguished only by the direction of curvature of the processes of the aedeagus are synonyms, that is, Farynala dextra Yan & Yang, 2017 equals Farynala sinistra Yan & Yang, 2017 syn. nov. A key to tribes of Typhlocybinae is provided.


| INTRODUC TI ON
The subfamily Typhlocybinae (Hemiptera, Membracoidea, Cicadellidae) is a large group of mostly tiny, delicate leafhoppers that feed primarily on leaf parenchymal cell contents of their host plants, thus differing from the phloem-or xylem-feeding (Cicadellinae) preferences exhibited by most other leafhoppers (Bartlett et al., 2018;Dietrich, 2013a). This group is distributed worldwide and comprises numerous agricultural pests (Nault & Ammar, 1989;Sun, 2004;Wearing et al., 2011). Examples include the potato leafhopper, Empoasca fabae, and tea green leafhopper, Matsumurasca onukii (Backus et al., 2005;Chasen et al., 2014;Qin et al., 2015). However, the vast majority of species, which feed on a wide variety of plants, appear to be of little or no economic importance. The group remains poorly studied, partly because their delicate nature makes them difficult to preserve for morphological studies. Based on the large number of described species (approximately 5000), Typhlocybinae is the second largest cicadellid subfamily (after Deltocephalinae; Bartlett et al., 2018) but the actual diversity of the group is probably much higher (Dietrich, 2013a). Typhlocybines are readily distinguished from other leafhoppers by the following morphological characters: forewing fully developed without closed anteapical cells; hind tarsomere Ⅰ acuminate, without transverse row of blunt setae (Dietrich, 2005).
The current tribal classification of Typhlocybinae is based almost entirely on a few characters of the wing venation. Compared with other cicadellid tribes, typhlocybines have the venation of the fore-and hind wings relatively reduced. Thus, particular patterns of reduction and loss or consolidation of veins have traditionally been used to define the different tribes. Alebrini, the tribe traditionally considered to be the most "primitive," is the only tribe with an appendix on the forewing, a trait shared with non-typhlocybine leafhoppers. Empoascini lack the section of the hind wing submarginal vein that extends along the costal margin. Erythroneurini and most Typhlocybini have the hind wing submarginal vein completely lacking at the wing apex. Erythroneurini and some Dikraneurini have hind wing vannal vein unbranched but Dikraneurini retain a complete submarginal vein. Young (1965) transferred Mileewini (including genera Amahuaka and Ujna) into Typhlocybinae based on intuitive morphological criteria, but he later suggested that Mileewinae should be treated as a separate subfamily (Young, 1968), a classification that has been followed by most subsequent authors. Phylogenetic analyses of Membracoidea based on morphology and DNA sequence data yielded inconsistent results. The concatenated maximum likelihood (ML) analysis of anchored-hybrid data by Dietrich et al. (2017) placed Eurymelinae (sensu lato) as sister to Typhlocybinae but with only 54% bootstrap support. An earlier morphology-based analysis of Cicadellidae (Dietrich, 1999) placed Mileewinae (in part, Mileewini) as sister to Typhlocybinae but the first molecular phylogeny of this family (Dietrich et al., 2001) did not consistently recover Typhlocybinae as monophyletic and its relationships to other subfamilies were also inconsistently resolved. A subsequent combined analysis of morphological and 28S rDNA sequence data supported the monophyly of Typhlocybinae and its sister relationship to Mileewini . The sister-group relationship receiving moderate bootstrap support (lower than 85%) was also recovered by Chen et al. (2021) and Lu et al. (2021).
As for most highly diverse and poorly studied insect groups, research on typhlocybine systematics continues to focus mainly on discovery and description of new taxa. Relatively few efforts have been made to assess the phylogeny of the group and test the monophyly of recognized taxa. Previous analyses have supported the monophyly of Typhlocybinae but its sister-group has remained uncertain (Balme, 2007;Dietrich et al., 2005;Lu et al., 2021). Zhang (1990) proposed an intuitive morphology-based hypothesis of tribal relationships within Typhlocybinae, considering Alebrini to be the earliest diverging lineage based on the retention of an appendix (shared with other leafhoppers) in the forewing, with Dikraneurini sister to Empoascini based on the relatively well-developed submarginal vein (smv) of the hind wing. In an unpublished dissertation, Balme (2007) conducted the first explicit cladistic analysis of the group, using 73 discrete morphological characters and two molecular markers (16S rRNA and Histone 3) and proposed a classification including four tribes, Alebrini + (Typhlocybini + Empoascini) + Dikraneurini, with Erythroneurini treated as a synonym of the latter tribe. An anchored hybrid enrichment-based phylogenomic analysis of Membracoidea based on 388 genetic loci and more than 99,000 aligned nucleotides  included only 1-2 representatives of each tribe but recovered Alebrini as sister to Empoascini and this clade sister to a clade comprising Typhlocybini and Erythroneurini + Dikraneurini with strong support. This dataset did not include a taxon sample large enough to test the monophyly of individual tribes or examine relationships among genera within tribes. Most recently, Lu et al. (2021) analyzed a much larger sample of typhlocybine taxa using data from fragments of three nuclear and two mitochondrial genes, recovering the same tribal relationships found by Dietrich et al. (2017) and also recovering Zyginellini as sister to Typhlocybini but with low support.
Here, we use sequence data from complete mitochondrial genomes in an attempt to improve resolution of phylogenetic relationships within Typhlocybinae and examine the evolution and stability of wing characters traditionally used for the classification of typhlocybine tribes. A total of 110 leafhopper mitochondrial genome assemblies have been previously deposited in GenBank. Among them, only 19 typhlocybine species are included representing the tribes Typhlocybini, Empoascini, Erythroneurini, and Zyginellini. Data for species of Alebrini and Dikraneurini were not previously available.
Therefore, prior to our study, mitogenome data for Typhlocybinae remained insufficient to facilitate a comprehensive phylogenetic analysis of the subfamily.
For this study, we assembled and annotated 28 new typhlocybine mitogenomes using next-generation sequencing (NGS) data, and conducted a comprehensive phylogenetic analysis including 19 previously available mitogenomes, representing all currently recognized tribes, to examine relationships among major lineages of the subfamily. We performed ancestral character state reconstruction (ACSR) to examine the evolution of key morphological characters.
We also used molecular divergence time methods to estimate the times of origin of various typhlocybine clades.

| Taxon sampling and DNA extraction
In this study, single individuals of 28 typhlocybine species were newly sequenced, including one representative of Alebrini, three from Dikraneurini, three from Zyginellini, and 21 from Typhlocybini.
Vouchers are deposited at the Institute of Entomology, Guizhou University, Guiyang, China (GUGC). Nineteen additional mitogenome sequences were obtained from GenBank. The 47 mitogenomes cover all currently recognized typhlocybine tribes (sensu Dworakowska, 1979). Based on results of prior analyses of Membracoidea, we selected representatives of Mileewinae (three species), Cicadellinae (five species) and Evacanthinae (two species) as outgroups (Balme, 2007;Dietrich, 2013b;Dietrich et al., 2010;Takiya, 2007). Included taxa including voucher number, GenBank accession number, and collection locality are listed in Table S2. All specimens were identified to species before DNA extraction. Genomic DNA was extracted from the legs and from the thoracic muscle tissue using the DNeasy Blood and Tissue kit (Qiagen, Germany) following the animal tissue protocol.
Finally, we generated three matrices for the tree inference: (1) amino acids sequence with the 13 protein-coding genes (PCGs_faa); (2) nucleotide sequence of 13 protein code genes with the third codon excluded (PCG12_fna); (3) nucleotide sequence of PCG12_fna plus the two ribosomal RNAs (PCG12_fna plus two rRNAs). Third codon positions were excluded from the nucleotide-based analyses to reduce the possibility of bias or long-branch attraction due to substitution saturation among species belonging to different genera (Leebens-Mack et al., 2005;Stefanović et al., 2004). We used both partitioned and nonpartitioned approaches for phylogenetic inference.
Partitioned maximum likelihood reconstructions were performed using IQ-TREE v1.6.3 (Nguyen et al., 2015) with 1000 ultrafast bootstrap (UFBoot; Hoang et al., 2018) and 1000 SH-aLRT replicates (Guindon et al., 2010). The best partitioning scheme and substitution models were selected based on MODELFINDER (Kalyaanamoorthy et al., 2017) implemented in IQ-TREE. Nonpartitioned reconstructions were made using site heterogeneous models in both maximum likelihood (ML) and Bayesian inference (BI). Posterior mean site frequency (PMSF) model (Wang et al., 2018) was used for the PCGs_ faa matrix by specifying a profile mixture model with the option "-mtInv+C60+FO+R" in IQ-TREE. The corresponding partitioned tree (PCGs_faa matrix) was treated as an initial guide tree. Bayesian inference using PhyloBayes MPI v1.8b (Lartillot et al., 2013) was performed for the PCGs_faa matrix as well. Two separate chains were independently run for 10,000 generations under the CAT+GTR model (Lartillot & Philippe, 2004) using a starting tree derived from PMSF ML analyses. We used the program bpcomp (maxdiff value) and tracecomp (minimum effective size) to check for convergence, that is, when the maxdiff value is smaller than 0.3 and minimum effective sizes are larger than 50. In discussing branch support, we consider values greater than 98% (SH-aLRT, UFBoot2) and 0.99 (posterior probability) represent "high" support; values of 80%-98% for SH-aLRT, 95%-98% for UFBoot2 and 0.95-0.99 for posterior probability indicate "moderate" support; and values lower than 79% for SH-aLRT, 95% for UFBoot2 and 0.95 for posterior probability are "low" support. Pairwise distances (p-distance) based on the nucleotide sequence of 13 PCGs are shown in Table S3.

| Ancestral character state reconstruction
Ancestral states were reconstructed for the following eight morphological characters that have been used previously to infer evolutionary trends within the group and to define tribes: (a) ocelli (absent  Table S4. The condition of character f is uncertain for species previously placed in Zyglinellini because the curved vein connected below the main stem of CuA could be interpreted either as a branch of CuA or as part of the submarginal vein, so this character was scored as uncertain (?) for the 6 included Zyglinellini species (Table S4). The description of each character and its states followed Dietrich (2013a) and Dworakowska (1993). Ancestral character state reconstruction (ACSR) was performed with a maximum likelihood approach using a single-rate Mk1 model in MESQUITE v3.20 (Maddison & Maddison, 2017). To account for phylogenetic uncertainty, we traced character history on 10,000 sampled posterior Bayesian trees and summarized results on the BI consensus tree. In addition, we also traced character history and summarized results on the ML tree based on the matrix of the PCG12_fna.

Because fossil leafhoppers have not yet been incorporated into ex-
plicit morphology-based phylogenetic analyses, we used available fossils to calibrate the root nodes of their respective tribes in order to provide conservative estimates of the minimum ages of these groups. Based on the previous molecular timetree of Membracoidea  and information on the oldest Membracoidea included in the Paleobiology Database (PDBD) Navigator website (http://paleo biodb.org/navig ator/), the maximum age of the root node was constrained with a relaxed lower bound of 174.1 million years ago (MYA). The fossil species Jassoqualus hispaniolensis from

Oligo-Miocene Dominican amber and Youngeawea bicolorata from
Eocene Baltic amber were used to provide minimum age calibration points for Evacanthinae and Mileewinae, respectively. Unidentified fossil species of the tribes of Cicadellini sp. and Dikraneurini sp. from Dominican amber were used to provide minimum age calibration points for Cicadellini and Dikraneurini, respectively ( Table 1). The input tree was retrieved from the PMSF topology. Two priors for the rgene gamma and sigma 2 gamma were set to (2, 20, 1) and (1, 10, 1), respectively. We used an independent rates model and alpha of 0.5.
The first 10 4 iterations were discarded as burn-in, and 10,000 samples were sampled every 5 iterations.  Dietrich and Vega (1995) Note: The "N/A" indicates not applicable.

| Morphology
The length of the body reported in the descriptions includes the forewings at rest. Male genitalia were prepared using the techniques described by Oman (1949). Morphological terminology follows Dietrich (2005Dietrich ( , 2011 and Dworakowska (1993). The genital segments of the specimens examined were macerated in 10% NaOH for 7-10 h (or boiled for 1-3 min), then transferred to glycerin after rinsing in distilled water to remove traces of NaOH (10%) for further research. Male genitalia were drawn using a LEICA m125 microscope and drawings were processed using Adobe Illustrator CS6. Habitus photographs were taken with a KEYENCE VHX-6000 digital camera and optimized with Adobe Photoshop CS6. All specimens studied are housed in the Institute of Entomology, Guizhou University, Guiyang, China (GUGC).

| Mitogenome features
Newly sequenced mitochondrial genomes of 28 species representing 21 genera in the subfamily Typhlocybinae ranged from 14,396 bp to 16,931 bp in length, including 37 typical insect mitochondrial genes and a control region (Table S2). Mitochondrial gene arrangement is highly conservative and consistent with the hypothetical ancestral insect (Figure 1a), except in the newly sequenced Shaddai sp. (GenBank, MW284820), which has the gene trnW translocated behind trnY, with 691 bp and 460 bp intergenic spacers between trnI and trnQ, trnY and trnW, respectively (dashed line box in Figure 1b).
This is the first report of a mitochondrial gene rearrangement within this subfamily.

| Phylogeny
The aligned amino acid dataset (PCGs_faa) included 3,467 sites. The aligned nucleotide sequence datasets included 6934 sites (PCG12_ fna) and 9218 sites (PCG12_fna+2rRNAs), respectively. The character partitions and models used in the partitioned analyses are shown in Table S5. A total of five phylogenetic trees were generated, including one Bayesian inference tree (maxdiff value = 0.11, minimum effective sizes >58; Figure S1) and four maximum likelihood trees ( Figures S2-S5). Overall, the trees resulting from different analyses are similar, with areas of incongruence limited to branches that received only low to moderate support in one or more analysis. The BI tree was selected as the preferred topology for subsequent analy-   (Figure 3). Many branches pertaining to relationships within tribes were consistently well supported in all analyses but some were unresolved, particularly within Typhlocybini; for example, relationships among four major lineages of Typhlocybini were resolved inconsistently across analyses and form a polytomy in the Bayesian consensus tree (Figure 3).
The phylogenetic results, combined with evidence from morphological characters (see below) and genetic distances (  (Qin et al., 2015;Yu et al., 2015) and analyses of phylogeny and genetic distance, we also show that the species previously identified as "Empoasca vitis (GenBank, NC_024838)" probably equals Matsumurasca onukii (COX1 p-distance = .3%). Other recent studies have shown that the latter species, a major pest of tea, has been widely misidentified in China (Qin et al., 2014(Qin et al., , 2015Yu et al., 2015). Unfortunately, we were unable to check voucher specimens from previous studies in order to confirm the identifications suggested by our results.

| Ancestral character state reconstruction
Eight typhlocybine morphological characters were traced over the BI consensus tree and ML tree (PCG12_fna). Overall, these reconstructions indicate that some characters used previously to diagnose tribes of Typhlocybinae exhibit considerable homoplasy (Figure 4;

| Divergence time estimates
A chronogram for Typhlocybinae divergence dates and topologies based on whole mitochondrial protein-coding genes is shown in

| Mitogenome features
Most leafhopper mitochondrial genomes are highly conserved with an arrangement of the 37 genes consistent with that of the inferred ancestral insect, Drosophila yakuba (Clary & Wolstenholme, 1985). However, more than ten species in Cicadellinae, Deltocephalinae, Iassinae,  Note: The value of deltaL indicates that logL differs from the maximal log1 in the comparison. bp-RELL, bootstrap proportion using RELL method (Kalyaanamoorthy et al., 2017;Kishino et al., 1990); p-KH, p-value of one sided Kishino-Hasegawa test (Nguyen et al., 2015); p-SH, p-value of Shimodaira-Hasegawa test (Kishino et al., 1990); c-ELW, Expected Likelihood Weight (Hoang et al., 2018). p-AU, p-value of approximately unbiased (AU) test (Wang et al., 2018). Plus signs denote the 95% confidence sets. Minus signs denote significant exclusion. All test performed on 1000 replicates using the RELL method.  Details of relationships within individual tribes are difficult to compare between our results and those of Lu et al. (2021) given that the taxon samples of the two studies overlap only partially.

F I G U R E 4 Ancestral character states reconstructions analysis based on
Nevertheless, we note that several branches pertaining to relationships among genera of Typhlocybini, as well as a few branches within the other tribes, were extremely short and received low bootstrap support in both studies. This suggests that neither complete mitogenome data nor data from a few nuclear and mitochondrial genes will suffice to completely resolve typhlocybine phylogeny with high confidence, although increased taxon sampling may also help improve phylogenetic resolution in this group.
Considering the morphological similarities shared by Typhlocybinae and Mileewinae, the latter hypothesis seems to be more plausible.
The main areas of uncertainty are the relationships between Typhlocybini and other tribes, and relationships among several deep nodes within Typhlocybini. Different analyses recovered Typhlocybini either as sister to Alebrini +Empoascini (Figure 2a,b) or as sister to Erythroneurini +Dikraneurini (Figure 2c, Ancestral state reconstructions of key morphological characters previously used to define and diagnose tribes within Typhlocybinae (Dworakowska, 1979(Dworakowska, , 1993Evans, 1963Evans, , 1971Hamilton, 1983;Zhang, 1990; Table 3) indicate that some wing characters traditionally used to diagnose tribes are highly stable. Our analyses generally support the monophyly of tribes Empoascini, Dikraneurini, Erythroneurini, and Typhlocybini as defined morphologically by most authors. Thus, our analyses suggest that the wing vein characters traditionally used to define these groups are reliable.
Nevertheless, our results also indicate that some morphological characters have undergone homoplastic changes during the evolution of the group. Ocelli may be present or absent within Empoascini and Typhlocybini. In Empoascini, all genera have ocelli, except Beamerana and Paulomanus. The latter genera, which have hind wing venation identical to that of many Empoascini, nevertheless grouped with Typhlocybini in a recent morphology-based phylogeny and were tentatively placed in the latter tribe (Xu et al., 2021).
Few genera in Typhlocybini have ocelli, for example, Hiratettix and Caknesia, but we were able to include only the former in our analyses and it occupies a relatively derived position within the tribe, suggesting that ocelli were lost and regained at least once in this tribe ( Figure 4a). Loss of the hind wing submarginal vein is also apparently homoplastic with independent losses apparently occurring in Erythroneurini and Typhlocybini (Figure 4e). Hind wing veins RP and MA have also apparently become confluent in the common ancestor of Erythroneurini and Dikraneurini and independently in Empoascini and Typhlocybini ( Table 3). Loss of the branched hind wing anal vein occurs only in the Dikraneurini +Erythroneurini lineage but further analyses with a larger sample of taxa will be needed to determine the extent of homoplasy in this character.

The molecular divergence time estimates indicate that
Typhlocybinae and extant tribes originated during the Middle and Late Cretaceous, respectively. Thus, diversification of major lineages in this group seems to have roughly coincided with the diversification of angiosperms during the Cretaceous (Foster et al., 2016;Hamilton, 1992Hamilton, , 1994Ledyard, 1974 Song et al. (2017) analyzed relationships within Cicadomorpha using mitochondrial genome data, including a species identified as "Typhlocyba sp." (GenBank accession number, KY039138). In our phylogenetic results this sample did not group with Typhlocyba choui, but instead grouped with the genus Aguriahana with high support in all phylogenetic trees (PP = 1, SH-aLRT = 100, UFB = 100; Figure 1). In addition, the COX1 p-distance between "Typhlocyba sp."

| CON CLUS IONS
This study sequenced 28 typhlocybine mitogenomes, representing all currently recognized tribes. We report the first known mitochondrial gene rearrangement within Typhlocybinae (Alebrini).

Remarks.
The new genus is similar to Typhlocyba in coloration and male genitalia but can be distinguished from the latter by the absence of basal macrosetae on the subgenital plate and presence of dense microsetae at the midlength of the paramere.
Etymology. This generic name is feminine, formed by the Latin word "Subtilissimus" which means "slender or fine," referring to its aedeagal shaft in lateral view.

Distribution. China (Yunnan)
Subtilissimia fulva  Etymology. The species name was derived from the Latin adjective "fulvus" = reddish-yellow, tawny, amber-colored, which refers to the coloration of the longitudinal patch on the crown.

Remarks.
The new species is very similar to Subtilissimia pellicula sp. nov. in size and coloration, but differs in having the aedeagal processes arising near the shaft mid length, longer, and not forked near the apex in lateral view (Figure 7g,h).
Etymology. The species name was derived from the Latin noun "pellicula" = "small skin", which refers to a membranous structure on the shaft. It should be treated as a noun in apposition.

Remarks.
The new species is very similar to Subtilissimia fulva sp. nov. in size and coloration, but differs in having the aedeagal processes arising near the shaft apex, curved dorsally, and forked near the apex in lateral view (Figure 7o,p).
Note. Yan and Yang (2017) described the species of Farynala sinistra and Farynala dextra, distinguishing them according to numbers of macrosetae on the pygofer side and the orientation of the aedeagal processes being opposite such that they appear to be mirror images of one another. Our current phylogenetic results (Figures S1 -S5) and the small genetic distance between these two taxa based on COX1 p-distance = .51% strongly indicates that these taxa are conspecific. Thus, the two species names are here treated as synonyms.

CO N FLI C T S O F I NTE R E S T
The authors declare no conflict of interest to disclose.