Specimens of A. psidii (about 50 adults) were collected from photinia (Photinia $\mathrm{\times }$ fraseri Dress) in an open field near Zhoukou City, Henan province, China (geospatial coordinates: 33.778°N, 114.859°E). All specimens were obtained from the same natural population. Morphological identification of the specimens was undertaken with reference to Yan and Bai [26] and Guastella et al. [27]. The specimens were preserved in absolute ethanol and stored at –80 ℃.

Genomic DNA was extracted from a pool of 40 specimens using the TIANamp Micro DNA Kit (Tiangen Biotech Co., ltd., Beijing, China). The DNA quality and concentration were examined with a nucleic acid protein analyzer (Quawell Technology Inc., San Jose, CA, USA).

Genomic DNA was sonicated to 300 bp using a Covaris S220 focused-ultrasonicator (Covaris Inc., Woburn, MA, USA), according to the Illumina protocol. Library generation for Illumina Hiseq sequencing was carried out using the Illumina TruSeqTM DNA Sample Prep Kit (Illumina, San Diego, CA, USA). The resulting fragment library was sequenced on an Illumina HiSeq2500 platform, with a strategy of 150 paired-end sequencing. Data quality control was performed using the NGS QC Toolkit v2.3.3 (New Delhi, India) [28], under the default settings. Adapters, ploy-N, and low-quality reads were removed from raw data. The high-quality reads (avg. Q20 $>$90%, and avg. Q30 $>$80%) were used to assemble the mitochondrial scaffold.

Two assemblers were used for mitogenome assembly, namely GetOrganelle v1.7.5.2 (Kunming, China) [29] and Geneious R11 (Auckland, New Zealand) [30]. For assembly with GetOrganelle [29], we used the GetOrganelle animal database (-F animal_mt) to identify, filter, and assemble target-associated reads. For assembly using Geneious [30], we applied the “map to reference strategy”, with the pre-sequenced mitochondrial cox1 gene sequence as the reference. The parameter settings were selected as described by Yang et al. [31].

The newly sequenced mitogenome of A. psidii as well as mitogenome sequences obtained from GenBank were annotated in the MITOS webserver [32] (http://mitos2.bioinf.uni-leipzig.de/index.py, Leipzig, Germany) to obtain uniform annotations. The RefSeq 63 Metazoa and the Invertebrate genetic code were selected as job settings. Mitogenome structure images for all whitefly species were generated using OGDRAW (Potsdam-Golm, Germany) [33]. Gene boundaries for protein-coding genes were refined by alignment with closely related whiteflies (Table 1). Sequences of protein-coding genes were separately aligned using TranslatorX (London, UK) [34]. A detailed description of parameter settings for TranslatorX [34] is available in the section on Sequence alignment. The secondary structures of 22 tRNA genes were predicted by MITOS [32]. The gene boundaries of two rRNA genes were refined by alignment against published sequences. The rRNA sequences were aligned separately using MAFFT version 7 (Osaka, Japan) [35]. The secondary structures of rRNA genes were inferred by reference to D. yakuba [36]. The annotated mitogenome sequence was submitted to GenBank under the accession number OM362911.

The nucleotide composition of mitogenome sequences was computed with MEGA X (Hachioji, Japan) [37]. AT and GC-skew values were measured using the following formula: AT-skew = (A–T)/(A+T) and GC-skew = (G–C)/(G+C) [38]. Non-synonymous (dN) and synonymous (dS) substitution rates for protein-coding genes were estimated using PAML 4.9 (London, UK) [39].

Comparison of the gene order in whitefly mitogenomes was performed in CREx (Leipzig, Germany) [40] (http://pacosy.informatik.uni-leipzig.de/crex/form), with the common intervals parameter for distance measurement. In addition, mitogenome structures were mapped onto the sequence-based phylogeny to explore the potential phylogenetic information from mitochondrial gene rearrangements.

Each of the 37 mitochondrial genes was aligned separately. Protein-coding genes were aligned using codon-based multiple alignments under the MAFFT algorithm [35] implemented in TranslatorX [34], with the Invertebrate Mitochondrial genetic code. Poorly aligned sites from the protein alignment before back-translate to nucleotides were removed with Gblocks v. 0.91 (Barcelona, Spain) [41]. The options used for less stringent selection were: allowing smaller final blocks, allowing gap positions within the final blocks, and allowing less strict flanking positions. Ribosomal and transfer RNA genes were aligned individually in the MAFFT server [35] and with the “E-INS-i” iterative refinement method. Alignments were concatenated with FASconCAT-G_v1.04 (Bonn, Germany) [42]. The partition option “-l” was selected to simultaneously generate a gene partition output file for the concatenated supermatrix. This partition file can be used directly or with slight adjustment in the subsequent Maximum likelihood analysis.

Three types of concatenated datasets were used in the phylogenetic analysis: (1) PCG_nt, nucleotide alignment including 13 protein-coding genes; (2) PCG_aa, amino acid alignment including 13 protein-coding genes; and (3) PCGRNA, nucleotide alignment including 13 protein-coding genes, two rRNA genes and 22 tRNA genes.

Phylogenetic trees were constructed using Maximum likelihood (ML) and Bayesian inference (BI) approaches. ML analyses were performed using IQ-TREE 1.6.10 (Vienna, Austria) [43]. Models for ML analyses (Supplementary Table 1) were chosen for each concatenated dataset using ModelFinder (Canberra, Australia) [44] implemented in IQ-TREE 1.6.10 [43]. For the nucleotide dataset of protein-coding genes, data blocks were predefined by gene and codon. For the amino acid dataset of protein-coding genes, data blocks were predefined by gene alone. The 22 tRNA genes were defined as a single partition due to very short gene length of each tRNA gene (60~70 bp), while two rRNA genes were defined as separate partitions. We used the option of allowing partitions to have different speeds (-spp). Nodal support values (BS) were evaluated using 1000 ultrafast bootstrap replicates [45]. BI analyses were conducted using PhyloBayes MPI on XSEDE (1.8c) (Montréal, Québec, Canada) [46], as implemented in the CIPRES Science Gateway [47]. The CAT-GTR model was employed for the nucleotide datasets PCG_nt and PCGRNA, while the CAT-mtZOA model was used for the amino acid dataset PCG_aa. All BI analyses involved two chains, with a total length of 10,000 cycles. Constant sites were removed. Convergence between the two chains was evaluated by examining the difference in frequency for all of their bipartitions (maxdiff 0.1). Trees were collected after the initial burn-in (10%) was discarded. A majority-rule consensus tree was computed using the program bpcomp implemented in PhyloBayes [46]. Branch support was assessed by clade posterior probabilities (PP).

A total of 197,079,804 clean read-pairs were obtained from the sequence library. GetOrganelle produced a 15,876 bp contig, while Geneious generated a 14,774 bp contig. For the mapping assembly conducted by Geneious, 28,046 of 197,079,804 reads were assembled, and the “Pairwise % Identity” reached 99.6%. The mean coverage of 14,774 bases was 286.2, with a standard deviation of 63.9. The minimum coverage was 2 and this occurred at the start site of assembly, while the maximum coverage was 448.

Alignment showed that the Geneious assembly was missing a region compared to that obtained from GetOrganelle. After annotation through MITOS, this missing region was found to span the majority of the rrnS gene and the first control region. With the exception of this missing region, Geneious presented an identical sequence to that from GetOrganelle. Therefore, the GetOrganelle assembly was used as the preferred result.

The A. psidii mitogenome contained the complete set of 37 genes usually present in insect mitogenomes. In addition, duplication of the tRNA gene cluster -trnQ_trnI_trnM was found in the control region. The mitochondrial gene order of A. psidii appeared to be unique among whiteflies. Notably, the relative position of nad1, which was nested within the control region, differed from that seen in all other whitefly mitogenomes. Moreover, the adjoining trnI and -trnQ, and the adjoining trnC and trnY exchanged positions with each other. The entire mitogenome had an AT content of 78.6%. On the major strand, the AT- and GC-skew values were –0.239 (higher T than A content) and 0.243 (higher G than C content), respectively.

A. psidii was unique in that additional start codons were used for protein-coding genes: GTG for atp6, GTT for atp8, and TTG for nad2. The remaining protein-coding genes used ATN (ATG, ATT and ATA) as the initiation codon. Eleven of the 13 protein-coding genes terminated with TAA. The abbreviated stop codon T was found for cox1 and nad1. Estimation of the evolutionary rates of protein-coding genes showed that the mitogenome of A. psidii evolved at a low rate (dN = 0.3117) relative to other whitefly species.

The program MITOS was able to locate 20 tRNA genes. trnS1 and trnL1 were identified by alignments with other closely related species. Secondary structures were predicted for all tRNAs (Fig. 1). A standard cloverleaf structure was inferred for 18 tRNA genes. trnC and trnA had T$\mathrm{\Psi }$C arm variable loops instead of the typical T$\mathrm{\Psi }$C arm. trnS1 and trnS2 lacked the DHU arm. Unpaired DHU arms of these tRNA genes were also found in other insects [1, 14, 15, 16, 17, 18].

Fig. 1.

The secondary structures of tRNA genes inferred for Aleuroclava psidii.

The rrnL and rrnS genes were located on the minor strand, with lengths of 1191 bp and 748 bp, respectively. Both these genes had very similar secondary structures to those of D. yakuba. The rrnL gene of A. psidii contained five domains (I–II, and IV–VI) with 40 helices (Fig. 2). The rrnS gene had three domains (I, II, and III) composed of 27 helices (Fig. 3).

Fig. 2.

The secondary structure of rrnL inferred for Aleuroclava psidii. Blue roman numbers denote the domains, while red numbers denote the helices.

Fig. 3.

The secondary structure of rrnS inferred for Aleuroclava psidii. Blue roman numbers denote the domains, while red numbers denote the helices.

The topologies resulting from ML and BI analyses were relatively concordant for the three concatenated datasets. Both Aleurodicinae and Aleyrodinae were monophyletic in all of these topologies with robust support (BS = 100, PP = 1). BI analyses of the datasets PCG_aa and PCG_nt resulted in an identical tree topology (Fig. 4, and Supplementary Figs. 1–3).

Fig. 4.

Phylogenetic tree inferred from the dataset PCG_aa using PhyloBayes under the model CAT-mtZOA (left), and comparisons of the gene order of whitefly mitogenomes (right). In the phylogenetic tree, the numbers around nodes indicate posterior probabilities. The scale bar represents substitutions/site. Bold denotes the newly sequenced species. The branch of Aleurodicus dispersus is depicted as half of its original length. Red indicates the mitochondrial gene rearrangements found in the whitefly species examined.

The main difference between the analyses was the placement of Trialeurodes vaporariorum and Neomaskellia andropogonis. On all BI trees, T. vaporariorum was recovered as the sister to all other Aleyrodinae. N. andropogonis was the sister to a clade comprising Aleurochiton and Pealius. In contrast, the phylogenetic placements of T. vaporariorum and N. andropogonis varied across ML analyses.

A sister relationship between Singhiella and Bemisia was recovered in all BI analyses (PP = 0.99 or 1) and in ML analysis of PCG_aa (BS = 100). At the genus level, the monophyly of Aleurodicus, Pealius, Aleurocanthus and Bemisia was supported in most analyses. The newly sequenced A. psidii was sister to Aleyrodes shizuokensis on the trees from PCG_aa and PCG_nt (BS $>$73, PP = 1). On ML and BI trees from the PCGRNA dataset, A. psidii was closely related to Aleyrodes shizuokensis.

The gene order in the mitogenomes of three psyllid species is identical to that of the ancestral insect mitogenome (Fig. 4). One of these (Arytainilla spartiophila) was selected to make pairwise comparisons with the mitogenomes of whiteflies using CREx analyses. All whitefly mitogenomes investigated in this study had experienced gene rearrangements. The newly sequenced A. psidii had three rearrangements, involving transpositions of trnQ, trnY, and the protein-coding gene nad1. Notably, the nested position of nad1 within the control region was very different to that observed in other whiteflies. At least one gene rearrangement was detected in other whitefly mitogenomes. The gene order in Aleurodicus dugesii was identical to that of Arytainilla spartiophila, except for the altered position of trnY. In contrast, the gene order was highly rearranged among representatives of two clades. One clade included Aleurochiton and Pealius, and the other contained Singhiella and Bemisia. Based on CREx analyses, the distinct gene order involving altered positions of protein-coding genes in Aleyrodidae sp., Aleurocanthus spiniferus, Aleurocanthus camelliae, S. simplex, B. afer, B. tabaci and B. emiliae may have evolved from the ancestral gene order following a tandem duplication random loss event.

Despite the significant economic importance and worldwide distribution of many well-known pests belonging to the family Aleyrodidae, few studies have investigated the phylogeny of this group. Here, we obtained the complete mitogenome sequence of A. psidii, and investigated phylogenetic relationships within Aleyrodidae. These results demonstrate the power of mitogenome sequences for establishing the relationships of whiteflies. The present mitogenome data clearly resolved Aleyrodidae into two monophyletic groups, namely the subfamilies Aleurodicinae and Aleyrodinae. All genera with more than two exemplars included were monophyletic, with the exception of Aleurochiton. The close relationship of Singhiella to Bemisia was confirmed, with both forming a sister group of a clade comprising Tetraleurodes and Aleurocanthus. This arrangement concurs with an earlier study of mitogenomes [48].

The phylogenetic incongruence observed here was mainly due to the unstable positions of T. vaporariorum and N. andropogonis. The long branches leading to T. vaporariorum and N. andropogonis highlight the rapid evolutionary rate observed in both species compared to other whiteflies. Moreover, one of the Aleurodicus species, namely Aleurodicus dispersus displayed the longest branches on all ML trees. Evolutionary rate analysis indicated the substitution rates for Aleurodicus dispersus, T. vaporariorum and N. andropogonis were higher than those of all other whitefly species (Table 2). In addition to estimating the evolutionary rate from the concatenated 13 protein-coding genes, we also calculated the evolutionary rate of each of the protein-coding genes in mitogenomes (Supplementary Table 2). The results showed that Aleurodicus dispersus, T. vaporariorum, N. andropogonis and S. simplex had higher non-synonymous substitution rates. In particular, Aleurodicus dispersus had the highest dN values for the cob, nad1, nad2, nad4, nad4l, nad5 and nad6 genes. Therefore, the different positions of T. vaporariorum and N. andropogonis were likely to be the result of long-branch attraction, especially in the analyses using the homogeneous model.

To further investigate the effect of long-branch attraction, we excluded the long branched Aleurodicus dispersus and reran ML analysis based on the PCGRNA dataset. This resulted in the long-branched N. andropogonis being in a different part of the tree and having a relatively close relationship to the clade including Aleurochiton and Pealius (Supplementary Fig. 4). This analysis confirms the potential attraction between Aleurodicus dispersus and N. andropogonis.

In contrast to ML analyses, the phylogenetic positions of T. vaporariorum and N. andropogonis were stable in the Bayesian inferences under the heterogeneous model. T. vaporariorum was retrieved as the first diverging clade in all BI analyses, while N. andropogonis was consistently placed as sister to a clade comprising Aleurochiton and Pealius. This result demonstrates the power of the heterogeneous model in suppressing long-branch attraction artefacts when using mitogenome sequence data to construct whitefly phylogeny.

Previous studies have suggested that mitochondrial gene rearrangements can be useful for resolving relationships within some taxonomic groups of insects [12, 49, 50, 51, 52, 53, 54, 55]. However, mitochondrial gene rearrangements are rare in Hemiptera. Nevertheless, the mitogenomes of single hemipteran groups often exhibit peculiar evolution that deserves careful examination. In this regard, whiteflies are among the most notable Hemiptera taxa, due to several interesting features. Most strikingly, several prior analyses of mitogenome sequences have reported that whitefly mitogenomes often have higher substitution rates than those of other hemipterans [56, 57]. In fact, the evolutionary rate of whitefly mitogenome sequences is too high to be used in some phylogenetic analyses [58]. In addition, whiteflies show an unpaired degree of mitochondrial gene rearrangements. Previous studies suggested that high rates of mitochondrial DNA substitution and of mitochondrial gene rearrangement were positively correlated in some insect groups [59, 60]. In the present work, we mapped the gene orders of mitogenomes onto the Aleyrodidae phylogeny based on analysis of amino acid sequences (Fig. 4). Comparisons of the gene order of mitogenomes evidenced the clade-specific evolutionary trend in whiteflies.

The mitogenome structures shown in Fig. 4 illustrate the evolutionary divergence patterns of mitochondrial gene rearrangements that frequently occur in the family Aleyrodidae. Originally, the trnW_-trnY_-trnC organization was retained in the ancestral mitogenome. However, in the new mitogenome of A. psidii, the -trnY gene was located between the trnW and -trnC genes. This pattern was shared by most whiteflies included in the present study. In addition, the tRNA genes located within the tRNA gene clusters and originally situated between the nad3 and nad5 genes and between the control region and nad2 genes can often swap positions in some taxa of Aleyrodidae. In particular, rearrangement of the tRNA gene cluster originally situated between the nad3 and nad5 genes can serve as a synapomorphy for several clades within the group. Besides the rearrangements of the tRNA genes mentioned above, the position of nad1 also changed in the mitogenome of A. psidii. This protein-coding gene was originally situated between the trnS2 and -trnL1 genes in the ancestral mitogenome. However, in the A. psidii mitogenome the nad1 gene was located between two control regions. This pattern was distinct from all other whitefly mitogenomes. Assembly from two programs (Geneious and GetOrganelle) resulted in the same position for the nad1 gene. Nevertheless, further sequencing of mitogenomes from the genus Aleuroclava is needed to confirm this result.

Aleurodicus dugesii displayed the closest gene order to the ancestral mitogenome, with only a single transposition of the -trnY gene detected in the Aleurodicus dugesii mitogenome. The trnW_-trnY_-trnC arrangement was generally seen across the entire family Aleyrodidae. Reverse transposition of the -trnL1_-rrnL_-trnV_-rrnS gene cluster to the unique gene arrangement rrnS_trnV_rrnL_trnL1 was found in Aleurodicus dispersus.

Similar to Aleurodicus dugesii, several early-diverging species within Aleyrodinae also had a small cluster of gene rearrangements, namely T. vaporariorum, Aleyrodes shizuokensis and A. psidii. The trnQ_trnI_trnM gene order was shared by T. vaporariorum, Aleyrodes shizuokensis and A. psidii, which may indicate a close relationship for these three species. In comparison, the remaining Aleyrodinae whiteflies had more extensively rearranged mitogenomes. In addition to the frequent rearrangements of tRNA genes, some of the gene arrangements involved protein-coding genes and/or rRNA genes. For N. andropogonis, the gene cluster cox3-trnG-nad3 originally located between atp6 and trnA and encoded on the major strand was translocated to a position between trnP and the control region, as well as being encoded on the minor strand. Furthermore, the rrnS gene was reversed from the typical minor strand to the major strand. Relative to D. yakuba, we identified an inversion in four species of Aleurochiton and Pealius, namely -trnR_-trnA_-nad3_-trnG_-cox3_-trnN_-trnQ located between trnS2 and -nad1. This could serve as a potential synapomorphy for the clade containing the two genera. The reverse position of trnS1 shared by N. andropogonis, Aleurochiton and Pealius may indicate a close relationship of N. andropogonis to Aleurochiton and Pealius.

The gene cluster -nad3_-trnG_-cox3 situated between -rrnS and the control region may be used as a character in support of the clade including T. acaciae, Aleyrodidae sp. and two species of Aleurocanthus. The -trnQ_-trnV cluster located between -rrns and -nad3 together with the -trnR_-trnD cluster between cox3 and the control region in T. acaciae, as well as the -trnN_-trnQ_-trnV_-trnR_-trnA cluster between -rrnS and -nad3 in Aleyrodidae sp. may be autapomorphic for these two species. The pattern of -trnN_-trnQ_-trnV_-trnR_-trnA located between rrnS and nad3 may be a synapomorphy for the genus Aleurocanthus.

The gene cluster -trnS1_-trnE_-trnF lying between atp6 and -nad5 and the gene cluster of -trnV_-trnD_-trnQ_-rrnS_-trnN_-trnR_-trnA_-nad3_-trnG_-cox3 located between -rrnL and the control region were shared by S. simplex and three species of Bemisia (B. afer, B. tabaci and B. emiliae). The transposition of -trnI and -trnW may be autapomorphic for S. simplex.