Comparative evolutionary analyses of eight whitefly Bemisia tabaci sensu lato genomes: cryptic species, agricultural pests and plant-virus vectors

Abstract

Background: The group of > 40 cryptic whitefly species called Bemisia tabaci sensu lato are amongst the world's worst agricultural pests and plant-virus vectors. Outbreaks of B. tabaci s.l. and the associated plant-virus diseases continue to contribute to global food insecurity and social instability, particularly in sub-Saharan Africa and Asia. Published B. tabaci s.l. genomes have limited use for studying African cassava B. tabaci SSA1 species, due to the high genetic divergences between them. Genomic annotations presented here were performed using the 'Ensembl gene annotation system', to ensure that comparative analyses and conclusions reflect biological differences, as opposed to arising from different methodologies underpinning transcript model identification.

Results: We present here six new B. tabaci s.l. genomes from Africa and Asia, and two re-annotated previously published genomes, to provide evolutionary insights into these globally distributed pests. Genome sizes ranged between 616-658 Mb and exhibited some of the highest coverage of transposable elements reported within Arthropoda. Many fewer total protein coding genes (PCG) were recovered compared to the previously published B. tabaci s.l. genomes and structural annotations generated via the uniform methodology strongly supported a repertoire of between 12.8-13.2 × 10³ PCG. An integrative systematics approach incorporating phylogenomic analysis of nuclear and mitochondrial markers supported a monophyletic Aleyrodidae and the basal positioning of B. tabaci Uganda-1 to the sub-Saharan group of species. Reciprocal cross-mating data and the co-cladogenesis pattern of the primary obligate endosymbiont 'Candidatus Portiera aleyrodidarum' from 11 Bemisia genomes further supported the phylogenetic reconstruction to show that African cassava B. tabaci populations consist of just three biological species. We include comparative analyses of gene families related to detoxification, sugar metabolism, vector competency and evaluate the presence and function of horizontally transferred genes, essential for understanding the evolution and unique biology of constituent B. tabaci. s.l species.

Conclusions: These genomic resources have provided new and critical insights into the genetics underlying B. tabaci s.l. biology. They also provide a rich foundation for post-genomic research, including the selection of candidate gene-targets for innovative whitefly and virus-control strategies.

Keywords: Biological species; Cladogenesis; Comparative genomics; Endosymbiont; Genome assembly; Horizontal genes; Phylogenomics; Transposons.

Fig. 1

Comparison of six de novo Bemisia tabaci s.l. genome assemblies. Genome assemblies depicted as circular plots, where the complete plot represents the full genome length. Plots highlight the longest scaffold; scaffold N50/N90; assembly GC and gap coverage (%). Genome assembly completeness (BUSCO v3.0) shown with Insecta (OrthoDB v9; n = 1,658) orthology set. Historical and current B. tabaci s.l. population names summarized in Table 1. Assembly plots generated with assembly-stats [38]

Fig. 2

Bemisia tabaci s.l. functional annotation and enriched GO terms. a Histogram comparison of gene percentages of significant GO terms, identified from eight B. tabaci s.l. genomes. Population names are color coded in the inset and B. argentifolii has been abbreviated to “B. argen”. b Significant differences in GO terms across eight B. tabaci s.l. populations. Y-axis shows log.¹⁰ transformed P-values for each GO term in figure part (a). Results generated using web service WEGO v2.0 (https://wego.genomics.cn/)

Fig. 3

Genomic transposable element content in Bemisia tabaci s.l. genome assemblies. a Summary of major TE classes highlighting copy count (#) and genome repeat coverage (%). b Stacked bar-chart of Kimura sequence divergences of TE classes, expressed as a function of percentage of each genome; Y-axis: Genome percent coverage (%); X-axis: Kimura divergence score. c 3D-Bar graph showing TE copy count of repeat classes: DNA, LINE, LTR, SINE

Fig. 4

A genome wide species level phylogeny with clade specific orthologs. Whole genome comparative analysis computed with Orthofinder. Publicly available whitefly species were B. argentifolii, B. tabaci s.s. and T. vaporariorum (“Greenhouse whitefly”). Phylogenetic relationships estimated with RAxML (maximum likelihood) and MrBayes (Bayesian posterior probability) on a concatenated matrix of protein sequences of 23 species covering 655 OGCs (131,953 amino acids). a Species-level phylogeny with associated node support values (*/*) ⟹ Bayesian PP / BS (bootstrap replicates n = 100); under the best-fitting substitution model LG + G + F + I. b Ortholog set delineation depicted with respect to major Arthropoda clades (Pancrustacea, Hexapoda, Hemiptera and Aleyrodidae); ‘Multi-copy N:N:N’ ortholog sets contain ≥ 1 gene across all species; ‘Patchy’: missing a single species representative. The six new B. tabaci s.l. populations are highlighted in green dashed boxes. c OGC clade sets with relatively low gene counts expanded for clarity

Fig. 5

Reproductive compatibility of eight B. tabaci s.l. populations collected in Uganda and Nigeria. Male parents (top row) and female parents (left column). Symbols represent the degree of reproductive compatibility. The black circle (⚫) represents complete reproductive compatibility between members of the SSA1-SG1 ∪ SG2, the fisheye circle (◉) represents complete reproductive compatibility between members of the SSA2 ∪ SSA3 species, the hexagon (⬢) represents complete reproductive compatibility observed in the SSA1-SG3 population, while the circled-cross ( ⊗) represents complete reproductive incompatibility with no female progeny production in F1 generation. The mating-crosses denoted by double asterisks (**) were carried out by Mugerwa et al. [8]

Fig. 6

Average Nucleotide Identity and genomic synteny among Candidatus Portiera aleyrodidarum from different Bemisia hosts The cladogram on the left summarizes Portiera relationships based on their pairwise Average Nucleotide Identity values (heatmap, middle). On the right, genomic synteny conservation among Portiera strains based on 202 complete Coding Sequences (CDS) (blue) and the CDS presence in the variable region (green). Portiera genomes are represented linearly, the presence of a subcircular conformation of the variable region is represented at the end of the plot (separated by double backslashes). Blue boxes representing syntenic CDS in the direct strand (upwards) or in the complementary strand (downwards), genes from the variable region are denoted in green. Gray lines connect orthologous CDS

Fig. 7

Protein identity across detoxification gene families. A box-plot representation of a curated set of detoxification proteins obtained from B. argentifolii and their putative orthologous protein (each represented by a single dot) in seven analyzed species of Bemisia tabaci s.l: Asia II-5, B. tabaci s.s., SSA1-SG1-Ng, SSA1-SG1-Ug, SSA2-Ng, SSA3-Ng and Uganda-1. A BLAST-combined with manual inspection approach used to check the identity of each protein. All alignments shown include putative orthologous proteins with at least 100 amino acids of the entire sequence aligned (cutoff >  = 85% PID). The number of proteins analyzed between species vary, as it was not always possible to recover a B. argentifolii orthologue in each of the seven analyzed Bemisia species

Fig. 8

Phylogenetic relationships of α-glucosidase (GH-13) genes of thirteen arthropod species. Phylogenetic analysis focused on three Bemisia tabaci s. l. populations: SSA1-SG1-Ug, SSA1-SG1-Ng and B. argentifolii. Non-whitefly taxa T. castaneum, A. pisum, A. gambiae, B. terrestris, B. mori, D. plexippus, D. pulex, D. melanogaster, R. prolixus and T. urticae are uniquely colored. Phylogenetic analysis performed using a Bayesian approach and implemented in Bayesian Evolutionary Analysis Sampling Trees (BEAST version 1.10.2). Clusters (C 1 <—> C 7) are defined based on selection analysis; see Additional file 1: Table S14. The α-glucosidase genes related to sucrose hydrolysis are located in cluster 2, highlighted in purple (**)

Fig. 9

Bemisia tabaci s.l. horizontally transferred genes (HTGs) of bacterial and fungal origin. Total number of HTGs of: (a) bacterial and (b) fungal origin, derived from protein-coding transcripts. Relative proportion of functional predictions (COG categories) for: (c) bacterial and (d) fungal HTGs. Venn diagram of the distribution of HTG orthologous gene clusters (OGCs) among whitefly genomes for: (e) bacteria and (f) fungi