1. Introduction:
The Bemisia tabaci (Gennadius 1889) (Insecta: Hemiptera: Aleyrodidae) complex of whiteflies has collectively colonized over 600 plant species and includes globally important pests [1,2]. Members of the lineages within this complex have been referred to as races, strains, and biotypes [3,4,5,6], or haplotypes and mitotypes [7,8]. Globally, only a few of these genetic variants are known to cause annual economic losses [9,10,11], whereas most lineages are benign in their native habitats in the absence of irrigated monoculture-agriculture [7,10,12]. Originally described from cultivated tobacco plants in Greece, B. tabaci has been defined to include numerous populations that inhabit mostly tropical/subtropical and climatically moderate temperate regions [1,7,8,13,14]. Despite the accumulated evidence of behavioral differences among selected lineages [3,4,5,15,16,17,18], all but one of these populations have been shown to be morphologically indistinguishable [8,19,20].
Molecular and population studies have demonstrated extensive molecular divergence between different lineages of B. tabaci occurring throughout the world [7,8,21,22,23]. Members of the B. tabaci complex can be divided between biogeographic regions into major clades based on mitochondrial cytochrome oxidase I (COI) gene divergences, which can range as high as 24% among clades [7,24]. The latter observation provides the basis for the hypothesis that B. tabaci comprises a complex of distinct, mostly cryptic species [7,14,17,21,22,25]. The question surrounding the number of species and their nomenclature has been controversial, with some research groups using a small fragment of the 3′-end of the COI gene (~600 bp) to delineate species boundaries that gave rise to 28 species [26,27]. While the genetic structure has been investigated for selected B. tabaci populations [22,28,29], the majority of studies have relied upon comparisons of a partial fragment (658–780 bp) of the COI gene [21,26,30,31,32,33]. Mitochondrial sequences evolve at a relatively high rate and have a low effective population size [34,35], which can be useful for understanding gene flow between populations [36,37,38] and for inferring past demographic events that have shaped the genetic diversity of the complex [39,40]. However, introgression of mitochondrial sequences into different nuclear phylogenetic backgrounds is possible [41], making it important to understand if lineages identified in mitochondrial data are also reflected across the nuclear genome. Only a few studies have extensively sampled nuclear data [22,42,43] for members of the B. tabaci complex, and whole genomic data sets have been lacking, precluding the analyses of important genome-level questions.
Numerous biotypes from the B. tabaci complex of whiteflies have been described that can be analyzed in a phylogenomic context. Although many biotypes have been described, not all of them may be evolutionary significant lineages or cryptic species. A few biotypes are featured prominently in the literature. For example, the “B” and “Q” biotypes likely have a similar Middle East/North African origin and have become highly invasive in a monoculture setting [7]. The “B” and “Q” types are thought to be a responsible for the displacement of the “A” biotype that was once widespread across North America [44]. Another significant population of whiteflies includes individuals from the sub-Saharan African region and have become an intense focus for scientific studies because of their detrimental impacts on cassava crops [45,46]. Attempts have been made to summarize the many biotypes [7]; however new types are often described and it is difficult to summarize these under a single framework. A review suggests that seven major phylogeographic lineages exist that could be distinct species [7] and the analyses performed in this current study samples six of these major phylogeographic clades. Unfortunately, there exist no international rules for zoological nomenclature regarding the description and reference to biotypes in scientific literature, therefore we refer to major lineages for our current study in the context of geographical acronyms suggested by a previous publication [7].
Next generation sequencing technologies provide a cost-effective method to obtain data for thousands of genes across the nuclear genome [47,48]. Extensive sequence data representing much of the genome is expected to provide needed information at the nuclear genome level to complement the insights gained from published mitochondrial COI datasets. A genome of the “B” biotype [4], has an estimated size of 615 Mbp [49,50]. Thus, the size of the B. tabaci genome makes it possible to obtain sequence coverage across the genome using Illumina sequencing technologies and shotgun sequencing at modest cost. This annotated B. tabaci genome was also used as a reference for read-mapping among members of the B. tabaci species complex to obtain sequences of single-copy orthologs.
To investigate the extent of genomic divergence among lineages of B. tabaci, the whole genome was sequenced using the Illumina sequencing platform for nine individuals of the B. tabaci complex. Members of the complex were selected to span a broad biogeographic representation of major clades previously recognized [7]. Read-mapping was used to obtain gene sequences for 2184 single-copy orthologs. The resultant data sets were analyzed using phylogenomic approaches (coalescence and maximum likelihood) to assess the concordance between nuclear and mitochondrial datasets and to evaluate overall divergence among major lineages. To assess the stability of ambiguous nodes, likelihood mapping [51] and quartet sampling [52] were employed to demonstrate the frequency of discordant topologies produced from the concatenated dataset. Distance methods were also employed to compare sequence divergence from the mitochondrial and nuclear data obtained from read-mapping [53].
Automatic barcode gap discovery (ABGD) methods [54] were used to test for the existence of putative species within nuclear and mitochondrial alignments. ABGD methods use barcode gaps that exist between interspecific and intraspecific sequence datasets to assess the existence of unknown species diversity within an alignment [54]. These methods became popular because of the use of COI as a barcode to assess unknown diversity across a wide range of animals. This method has also been used to test for the presence of operation taxonomic units in a phylogenetic dataset containing possible cryptic species [55]. Although this method has traditionally been used on relatively small amounts of sequence data (i.e., COI alignments), the analysis is computationally efficient and has the ability to analyze large amounts of genomic data quickly. The analyses are based on the realistic observation that gaps exist between interspecific and intraspecific sequence divergences [54]. Thus, ABGD analyses have the ability to efficiently quantify these observed gaps and provide an estimate of putative species in genomic datasets.
The goal of this study was to analyze Illumina whole genome data to understand possible species boundaries within the B. tabaci complex of whiteflies. This was done through a series of maximum likelihood, coalescent, likelihood mapping, quartet sampling, and ABGD analyses. The genomic analyses of 2,184 genes provide further evidence to complement previous studies that indicate multiple cryptic species exist within the B. tabaci complex of whiteflies.
2. Materials and Methods
2.1. Sample Preparation and Sequencing
Based on partial mitochondrial COI DNA sequences, members of the B. tabaci complex have been shown to form phylogeographic clades of seven major groups (Table 1) [7]. A range of genetic diversity exists within each group, but divergences among groups are high, potentially marking species boundaries [32]. Nine B. tabaci individuals belonging to six phylogeographic groups, previously identified based mostly on COI divergences (Table 1: mitochondrial major clades I-VI) [7], were selected for sequencing. Three whitefly samples were each collected from the Democratic Republic of Congo (Lulimba), Tanzania (Dar es-Salam), and Uganda (Bukoba), respectively, to represent the Sub-Saharan Africa (SSA) clade. One sample from Sudan (Wad Mani) was selected as the representative of the North Africa–Mediterranean–Middle East (NAF–MED–ME) clade. For the Asian–Pacific–Australia (AS–PAC–AU) region, a representative was collected from Hainan, China. To serve as a representative of the Asia major clade, an individual from Ahmedabad in India was selected. To represent the American Tropics/Caribbean (AM-TROP) clade, one each were selected from Puerto Rico, Arizona, and Ecuador (Table 2). A geographical acronym is applied to each phylogeographical major clade [7] and used throughout to refer to the B. tabaci lineages (Table 1 and Table 2). One adult B. afer individual was sequenced and included as the outgroup.
Whitefly adults were collected alive from plants and placed directly into vials containing either 70 or 95% ethanol. Specimens were shipped to the University of Arizona, and immediately transferred to 95% ethanol for storage at −20 °C. Individual adult whiteflies were ground in 1.5 µL tubes. The DNA was extracted by a Qiagen DNA Micro Kit (Qiagen, Valencia, CA, USA) using a 48-h incubation period, and the DNA was suspended in 52 µl of elution buffer. Extracted genomic DNA was quantified with a Qubit 2.0 fluorometer (Invitrogen, Carlsbad, CA, USA) using the manufacturer’s recommended protocol and reagents.
Paired-end shotgun genomic libraries were prepared with Hyper Library construction kits from Kapa Biosystems (Roche). The B. tabaci libraries were quantified by qPCR and sequenced for 151 cycles from each end of the fragments on an Illumina HiSeq 4000 machine using a HiSeq 4000 sequencing kit version 1 or a NovaSeq 6000 machine. Paired-end reads were 150 nt in length for the B. tabaci specimens. The B. afer library was quantified by qPCR and sequenced on one lane for 161 cycles from each end of the fragments on a HiSeq2500 using a HiSeq SBS sequencing kit version 4. Paired-end reads were 160 nt in length for the B. afer. Multiplexing was performed to achieve between 240 and 400 Gbp of data per sample. Library preparation and sequencing took place at the W.M. Keck Center (University of Illinois, Urbana, IL, USA). The FASTQ files were generated and demultiplexed with the bcl2fastq v 2.17.1.14 or 2.20 Conversion Software (Illumina). Raw reads were deposited in the sequence read archive (SRA) of the National Center for Biotechnology Information (NCBI) database. Adaptors were trimmed and sequencing quality reports were generated with FastQC v 0.11.7 [71].
2.2. Ortholog Prediction and Phylogenomic Analyses
The genome of the “B” biotype has previously been sequenced, assembled, and annotated [49,50]. We used version 1.1 of the genome annotation as a reference sequence for read-mapping of orthologs in newly sequenced individuals [50]. Using the OrthoDB v7 database [72] to subsequently serve as references for read-mapping, 2193 single-copy orthologs were identified. This ortholog set was previously developed for phylogenomics of hemipteroid insects [73]. We read-mapped our newly sequenced B. tabaci genomes against the 2193 reference orthologs using Bowtie2 v 2.3.4.1 [74] and SamTools v 1.7 [75], following a general read-mapping pipeline previously published [55]. The sequence alignment map (SAM) files produced by Bowtie2 were compressed to their binary versions (BAM), which were further sorted and indexed using SamTools. All mapped reads were retained in the BAM files by using the -F 4 option in SamTools. Pileup files were created from the BAM files with the mpileup command in SamTools using a maximum depth of 75, minimum base quality of 28, and minimum mapping quality of 10. The pileup files were then converted to VCF and FASTQ files using BCFtools and vcfutils.pl respectively. The same read-mapping pipeline was used for obtaining the sequences of the mitochondrial COI gene for each individual. The single-copy orthologs and COI data from the “B” reference genome used for read-mapping were included in the final alignments and phylogenetic analyses. The ortholog sequences for each individual genome were aligned with PASTA v 1.8.0 by nucleotide sequence using an increased memory capacity of 2048 MB [76]. Alignments were then masked using a 40% gap threshold with trimAl v 1.4 [77]. A concatenated supermatrix of aligned genes (orthologs) was produced with Sequence Matrix v 100.0 [78]. An uncorrected distance matrix was produced from the final concatenated supermatrix of 2184 genes using PAUP* v 4.0 [53].
The supermatrix of aligned nucleotides was analyzed using several different methods for comparison. A concatenated maximum likelihood (ML) analysis was completed with RAxML v 8.2.11 using a GTR + G model of evolution and 100 rapid bootstrap replicates [79]. RAxML does not implement simpler than GTR models of evolution and invariant sites were not considered due to computational non-identifiability between the alpha parameter and gamma distribution. An optimal partitioning scheme was determined with PartitionFinder v 2.1.1 with branch lengths linked, GTR + G model, BIC model selection, rcluster search, and rcluster-max 100 [80]. The resulting partitioning scheme was implemented in RAxML using a GTR + G model and 100 rapid bootstrap replicates.
Because relationships among species could be confounded by incomplete lineage sorting among different genes, a coalescent approach was also implemented to compute a species tree from the individual gene trees [81,82,83]. All 2184 gene alignments were each individually analyzed in RAxML with a GTR + G model and 100 bootstrap replicates. These derived bipartition trees were summarized in a coalescent analysis with ASTRAL v 5.5.9 using default parameters, and clade support assessed with local posterior probability [81,83]. A coalescent analysis was also performed with ASTRID v 1.4 using multi-locus bootstrapping clade support and input files from concatenated RAxML bipartition files [82].
To provide a comparison with data from the nuclear gene sequences, a 1466 bp fragment of COI from the “Q” biotype (Genbank accession: KJ591614.1) was used as a reference for read-mapping to obtain the COI sequences for each of the genomes that were sequenced. In addition, reads from the B. afer library were mapped against a published B. afer mitochondrial genome to obtain the COI sequence for this species (Genbank accession: KR819174.1). To verify the identity of each genetic variant, the resulting COI sequences obtained from read-mapping were searched with BLAST against NCBI using Megablast, implemented in Geneious v 11.1.3 (Biomatters, Inc., Newark, NJ, USA) [84,85,86]. The BLAST results were evaluated based on best hit by percent sequence identity with ties broken by sequence length. The COI sequence from the “B” reference genome [50] (Zhangjun Fei, Boyce Thompson Institute, Cornell University) was included for subsequent ML analyses. The COI sequences were aligned with ClustalW v 2.1 [87] using the default parameters as set in Geneious [86]. The alignment was trimmed to a length of 1459 bp. A ML analysis of the COI data was completed with RAxML using a GTR + G model and 1000 rapid bootstrap replicates. Uncorrected pairwise distances were calculated in PAUP* for the COI alignment (1459 bp). A second ML analysis was performed utilizing a publicly available dataset of COI data [88]. The available data set was aligned with the mitochondrial sequences obtained from read-mapping. This alignment was then trimmed to a length of 696 bp and a ML analysis was performed with RAxML using parameters as previously described.
To assess frequencies of discordant topologies that are produced from the data set, quartet sampling [52] and likelihood mapping [51] was performed. Quartet sampling methods were implemented using a published available script [52] utilizing the concatenated nucleotide alignment and resulting unrooted nucleotide topology as an input tree for analyses. The number of replicates per branch were set to 200 and the log-likelihood cutoff was set to 2, as suggested by published protocol [52]. Likelihood mapping was performed in IQ-Tree v 1.6.5. [89] using the concatenated nucleotide alignment as the input. All number of quartets possible were drawn, subsequent tree searches were skipped, and a GTR + G model was implemented. The nexus file containing taxon clusters were defined as the following: cluster 1 = B. afer, cluster 2 = reference “B” reference genome, cluster 3 = B. tabaci sampled from Sudan, and cluster 4 = all other ingroup whiteflies sampled.
As a test for putative species an ABGD test [54] was implemented for both the mitochondrial and nuclear data sets. First the test was conducted using the COI alignment under both Jukes-Cantor (JC69) and simple (i.e., uncorrected) distances. The relative gap width (X) was set to 1 because the program failed to run with the default setting (X = 1.5) for the COI data. Remaining settings for the COI analysis were set to the default. Secondly, the ABGD analysis was implemented on the nuclear supermatrix of single-copy orthologs using both JC69 and simple distances with default settings. The B. afer specimen was excluded from the alignment for ABGD analyses because of its high divergence from the ingroup.
3. Results
Estimated sequence coverage of the genome for all samples was 39–65X based on calculations using an estimated draft genome size of 615 Mb [49,50], number of reads, and length of reads (Table 2). Read-mapping produced consensus sequences for nearly all of the 2193 single-copy orthologs for members of B. tabaci complex, but only 1929 genes mapped using the more divergent B. afer genome (Table 2). From the read-mapping results, 2184 genes successfully aligned in PASTA. The final nucleotide supermatrix consisted of 3,673,822 bp of data (109,699 informative). Uncorrected distances from the complete nuclear gene alignment ranged from 0.12%–2.50% between members of the B. tabaci complex (Table 3A).
The results of phylogenomic analyses based on the 2184 gene data set were mostly consistent but with some differences found between concatenated ML and coalescent analyses. Samples from Sub-Saharan Africa (SSA) and the American Tropics (AM-TROP), respectively, each form separate monophyletic clades with high support across all analyses (Figure 1a–d). Only one individual each were sampled from AS–PAC–AU region and from the ASIA region, but together are recovered as a highly supported monophyletic lineage across all analyses. Results of concatenated ML analyses recovered that the “B” reference genome was sister to the remaining ingroup B. tabaci individuals (Figure 1a,b), but the coalescent (ASTRAL and ASTRID) analyses recovered a sister relationship between the “B” reference genome and Sudan individual (Figure 1c and Figure S1). The coalescent analyses of the whole genome data set are the only phylogenomic analyses supporting a monophyletic NAF–MED–ME group which includes the “B” reference genome and the Sudan individual sampled.
The phylogenetic analysis based on COI sequences was largely concordant with the nuclear tree. However, node support for relationships among these major groups was low, as might be expected from a single, rapidly evolving gene (Figure 1d). When our read-mapped COI sequences were analyzed in a ML framework with publicly available data [69], individuals sampled were recovered within clades of similar geographic origin (Figure S2). BLAST results also validate read-mapped COI data obtained and the best BLAST matches are reported (Table 2). The uncorrected distances for COI sequences within the B. tabaci complex ranged from 1.31–17.0% among the SSA, NAF–MED–ME, ASIA, AS–PAC–AU, and AM-TROP lineages (Table 3B). In general, divergence in COI sequences increased with increasing nuclear divergence (Figure 2). However, COI eventually became saturated, with evidence of multiple substitutions at high divergences with respect to nuclear divergence.
Likelihood mapping and quartet sampling on the concatenated nucleotide alignment showed some discordance which may have contributed to the conflicting coalescent and concatenated topologies. The majority of quartets produced from likelihood mapping (75%) support a sister relationship between the Sudan individual and all other ingroup whiteflies as reported for the ML concatenated analyses. However, 25% of the quartets sampled favor a sister relationship between the reference “B” reference genome and the Sudan individual (Figure 3). Quartet sampling methods also showed discordance at this described node (0.19/0/1) and for the ambiguously supported node for a monophyletic SSA + (AS–PAC–AU + ASIA) (−0.20/0/1) (Figure 4).
Results of all the ABGD tests, which uses gaps in mismatch distributions of pairwise distances, suggested the existence of five candidate species in our data set. The candidate species were divided into the following putative species groups: (1) SSA (Uganda, Tanzania, and Democratic Republic of the Congo), (2) NAF–MED–ME (“B” reference genome and Sudan), (3) ASIA (India), (4) AS–PAC–AU (China), and (5) AM-TROP (Arizona, Puerto Rico, Ecuador). Based on these exemplars, the candidate species were consistent, regardless of distance method used for analysis of the COI data and nuclear supermatrix of whole genome data.
4. Discussion
For this study, one B. afer and nine B. tabaci genomes were sequenced using Illumina high throughput sequencing technologies. Coverage for each genome was approximately 39–65X based upon an estimated draft genome size for B. tabaci (615 Mb) [49,50]. This level of coverage allowed for the successful mapping of reads to obtain sequences for over 2184 nuclear single-copy orthologs. Phylogenomic analyses of these sequences revealed evidence of biogeographically partitioned lineages that were also consistent with those identified by the mitochondrial COI gene sequence [7]. Methods to detect gaps (ABGD analysis) in sequence divergence were consistent between mitochondrial and nuclear data, providing additional evidence that these major groups represent multiple species.
The differences among uncorrected pairwise distances of the mitochondrial COI alignment are large relative to nuclear orthologs obtained from whole genome sequencing (Table 3, Figure 2). The mitochondrial divergences between members of the B. tabaci complex are high when compared to divergences between other congeneric species of other insect taxa (e.g., aphids: 0.5–7% [90]; psyllids: 3–10% [91,92]; bark lice: 7–9% [93]). Specifically, up to 17% mitochondrial divergence was observed between members of the B. tabaci cryptic species complex. The extensive observed inter-clade COI divergence was proposed to reflect an ancient evolutionary age, hypothesized to be about 86 million years ago (MYA) [26]. However, within insects, several haplodiploid groups or those with paternal genome elimination [94,95,96,97,98] have been shown to have substantially higher mitochondrial DNA substitution rates. For instance, in thrips [99], scale insects [100], book lice [101], and parasitic lice [102], divergences among congeneric species range from 10–25%. Whiteflies are also haplodiploid [103] and elevated mitochondrial substitution rates may be an underlying biological artifact. Unlike the mitochondrial COI gene, the total divergence across the nuclear alignment was only 2.5% among the B. tabaci exemplars. The highest divergence in COI sequences observed between the B. tabaci samples that were not delineated as separate species by the ABGD analysis was 8.25% (Table 3B). This divergence was between the Democratic Republic of the Congo and the two other Sub-Saharan (Tanzania and Uganda) exemplars. However, the equivalent nuclear divergence among these three samples was only 0.3–0.5% (Table 3A). Similarly, the “B” reference genome and the Sudan individual differed by 6.37% for the COI gene, but only 0.57% in nuclear genes. In contrast to sequence divergence measured, lineages that were demarcated as different species by the ABGD analyses consistently showed greater than 13% and 1% divergence in the COI and nuclear data, respectively.
Phylogenetic relationships obtained from the analyses of whole genome and COI data showed some variation in regard to the placement of certain members. A monophyletic relationship between the “B” reference genome and Sudan individual was recovered by coalescent analyses and the ML COI analysis (Figure 1c–d). However, concatenated analyses of nucleotide data favor an alternate topology with the “B” reference genome sister to all remaining Bemisia members. The sister taxon of the ASIA and AS–PAC–AU groups also varied between the concatenated and coalescent analyses. The concatenated ML analyses recovered a sister relationship between the two lineages that overlap in the Asian continent (ASIA and AS–PAC–AU) and the SSA lineage, albeit with relatively weak support (Figure 1a,b). However, the coalescent analyses resulted in a sister relationship between the AM-TROP and the ASIA + AS–PAC–AU lineages (Figure 1c).
Likelihood mapping and quartet sampling methods demonstrate the discordance in topologies produced from the concatenated nucleotide alignment, which helps explain the difference in results between coalescent and ML concatenated results. Likelihood mapping showed 75% of quartets support the topology produced from the ML concatenated analyses, but the remaining 25% support a sister relationship between the “B” biotype reference genome and the Sudan individual (Figure 3). Quartet sampling showed similar discordance for this node with a weak majority of quartets supporting the topology tested (QC = 0.19), but no skew in discordant topologies were detected (QD = 0) (Figure 4), suggesting only one discordant topology is favored similar to the result of likelihood mapping (Figure 3). Quartet sampling methods also display ambiguous support for the monophyly among ASIA, AS–PAC–AU, and SSA individuals with a weak majority of quartets supporting the alternate topology recovered by coalescent analyses (QC = −0.20) (Figure 4). The discordant frequencies detected from quartet methods demonstrate the instability of this node as suggested by coalescent analyses. The conflicting topologies reported between coalescent and concatenated ML analyses can be explained by incomplete lineage sorting and non-random introgression [81,82,83]. This would not be unexpected given recent analyses based on single nucleotide polymorphism data that suggest historical introgression patterns may have played a role in diversification of the complex [42].
The results of the analyses performed corroborate recent findings that suggest that the Sub-Saharan African group of whiteflies that occur on cassava form a single distinct lineage [45,46]. All analyses recover the three sampled Sub-Saharan individuals together in a monophyletic group, but there is discordance between analyses within the Sub-Saharan clade suggesting incomplete lineage sorting among Sub-Saharan populations. These results are consistent with the recent finding that introgression exists within this lineage [45]. A recent dating analysis based on genomic data estimates a divergence of the Sub-Saharan African population from the “B” (reference genome) to be approximately 5.26 million years ago [46], much younger than previous estimates based on COI data [26]. The nuclear divergences produced from this study may be an accurate reflection of this recent divergence within the complex (5 MYA). Therefore, dating analyses based solely on COI data may overestimate species divergences.
This study used whole genome sequences from shotgun sequencing methods to analyze members of the B. tabaci complex within a phylogenetic context. The read-mapping pipeline employed in this study was highly efficient in the recovery of orthologous sequences for members of the B. tabaci complex. Of the 2193 genes used as a genomic reference, nearly all were recovered with the read-mapping pipeline (Table 2). The read-mapping pipeline identified 1938 orthologs from the B. afer specimen, but this individual was meant to serve as an outgroup for phylogenetic analyses. However, all recovered COI sequences were validated using the described BLAST search protocol including the divergent B. afer individual (Table 2). The ML COI analysis on publicly available data [88] also validate read-mapped mitochondrial sequences obtained and are recovered within clades which correspond to geographic origin or biotype (Figure S2).
The results demonstrate that there is nuclear divergence among populations, however divergences among the single-copy nuclear orthologs analyzed are lower relative to COI sequences. The five putative species suggested by ABGD analyses complement phylogenetic analyses, suggesting population structure that divides members of the B. tabaci complex. High mitochondrial divergences may be a biological artifact that persists within the complex; thus, analyses of only COI data may overestimate potential diversity when analyzed in a phylogenetic framework. By exploring the nuclear divergences among orthologs identified with read-mapping and comparing them to results produced from analyses of COI data, evidence is growing for future recognition of multiple species within the complex. Although the analyses completed within this study are no panacea for the discord among the perception of B. tabaci taxonomy, the published genomic data set will no doubt provide an important resource to complement existing lines of evidence that putative, but undescribed species exist within the B. tabaci complex of whiteflies.
5. Conclusions
The B. tabaci complex of whiteflies has long been thought to comprise multiple undescribed species, however much disagreement exists on how many species should be recognized. The ABGD analyses suggest that at least five species exist both from the analyses of nuclear orthologs and COI data obtained. However, many populations remain under sampled and there is potential for the existence of unknown B. tabaci diversity within the complex. To avoid confusion regarding the recognition of species, caution should be taken regarding official taxonomic description of species, until a reasonable consensus is reached considering present and future evidence.
Supplementary Materials
The following are available online at
Author Contributions
R.S.d.M. contributed to methodology, formal analysis, visualization, original draft preparation, review and editing. J.K.B. and J.R.P.-M. contributed to conceptualization, resources, review and editing. A.D.S. contributed to software, methodology, review and editing. R.M.W. contributed to methodology, funding acquisition, review and editing. K.P.J. contributed to funding acquisition, supervision, review and editing. K.K.O.W. contributed to methodology, data curation, review and editing.
Funding
Support for this work was provided to K.P.J. by the National Science Foundation (grant number DEB-1342604). Support for R.M.W. was provided by the Swiss National Science Foundation (grant number PP00P3_170664).
Acknowledgments
We thank P. Gillespie for providing the sample of B. afer, and Ali Idris, Everlyne Wosula, Clerisse Casinga, and Peter Sseruwagi for providing the B. tabaci from Sudan, Tanzania, DR Congo, and Uganda, respectively. We thank Ray Gill (previously, California Department of Food and Agriculture, CA), for providing morphological identification of all whitefly samples. We also thank the three anonymous reviewers for their helpful comments.
Conflicts of Interest
The authors declare no conflict of interest.
Figures and Tables
Enlarge this image.Figure 1. Colors correspond results suggested by automatic barcode gap discovery (ABGD) tests for putative species. (A) Result of the concatenated maximum likelihood (ML) analysis of all nucleotide sites from the nuclear supermatrix (lnLikelihood = –6919264.24). Scale bar represents substitutions per site. Clade values are depicted as bootstrap support. (B) Result of the partitioned ML analysis of the nuclear supermatrix (lnLikelihood = –6907791.99). Scale bar represents substitutions per site. Clade values are depicted as bootstrap support. (C) Result of ASTRAL coalescent analysis of 2184 gene trees. Scale bar represents coalescent units. Clade values are depicted as local posterior probabilities. (D) Result of ML analysis of cytochrome oxidase I (COI) sequences obtained from read-mapping (lnLikelihood = −6174.94). Scale bar represents substitutions per site. Clade values are depicted as bootstrap support.
Enlarge this image.Figure 2. A scatter plot of uncorrected pairwise distances of COI versus uncorrected pairwise nuclear divergence for each whitefly specimen analyzed.
Enlarge this image.Figure 3. Likelihood map produced from IQ-tree showing the distribution of quartets that support a sister relationship between the reference “B” genome and the Sudan sample produced from the analysis of the concatenated alignment. Results show that 75% of quartets sampled support the Sudan sample together with the remaining samples (not including the reference) while 25% support the Sudan sample together with the reference, when all nucleotide sites are analyzed. These quartets are also relatively decisive in that no quartets provide ambiguous information (i.e., no points in the center or edges of the triangle).
Enlarge this image.Figure 4. The B. tabaci complex of whiteflies showing the results of quartet sampling (QC: quartet concordance (frequency a concordant topology is observed)/QD: quartet differential (scaled fraction of discordant topologies)/QI: quartet informativeness (proportion of replicates that are informative)). Red asterisk indicate nodes of instability as discussed in the text. Nodes represented by 1/NA/1 are most stable and no discordant topologies are detected across quartets sampled (QD = NA).
Classification of the B. tabaci sibling species group based on nuclear (lineages), mitochondrial (major clades), and esterase/RAPD-PCR patterns (biotype).
| Nuclear Lineages (Current) | Mitochondrial Major Clades a | Biotypes d | |||
|---|---|---|---|---|---|
| Name | Acronym | Number | Name | ||
| Sub-Saharan Africa | SSA | I | Sub-Saharan Africa—East and South | S | |
| II | Sub-Saharan Africa—West | ||||
| North Africa–Mediterranean–Middle East | III | North Africa–Mediterranean–Middle East | B, B1, B2 | ||
| NAF–MED–ME | J, L, Q | ||||
| Non-Ug/MS | |||||
| Asia 2 | ASIA | IV | Asia 2 | K, P, ZHJ2 | |
| ZHJ1 | |||||
| G | |||||
| Cv, I | |||||
| Asia–Pacific–Australia | V | Asia–Pacific Islands–Australia | H, M, Na (Nauru), | ||
| AS–PAC–AU | An | ||||
| ZHJ3 | |||||
| American Tropics: North-Central America and Caribbean | AM-TROP | VI | American Tropics: North and Central/ Caribbean | A, A1, A2, C, D, F, O, N, R | |
| VII | American Tropics: South America b | - | |||
| Unassigned groups c: | |||||
| Uganda | - | ||||
| Italy | T | ||||
| Benin | E | ||||
a Major clade nomenclature based on published classification of the B. tabaci species complex [7]. b A representative whitefly from this clade was removed due to low-quality alignments. c Unstudied groups that represent potential new lineages. d Biotypes A–Q [5,20,56]. Biotype O described [57], but later referred to as biotype N [56]. Biotype R identified by Ian Bedford [58]. Biotype S [59]. Biotype T identified in Nebrodi, Italy [60] and further characterized [61]. The MS ‘biotype’ was identified based on RAPD-PCR patters and COI analysis [62]. Biotypes An and Na from Australia and the Pacific Islands were designated based on RAPD profiles [63]. The biotype Cv [64,65,66]. The ZHJ1 and ZHJ2 [67,68,69]. The ZHJ1-3 biotypes [70].
Table 2Summary of collection information and genomic coverage for each individual sequenced. The last column (% Pairwise Identity) corresponds to the best BLAST hit match from the Megablast search of COI data obtained from read-mapping. Plant host indicates the scientific name of the plant the sample was collected. The color of rows corresponds to the putative species suggested by ABGD analyses.
| Geographical Acronym | Country | State/Province | City | Plant Host | Reads | Gbp | Coverage (X) | Mapped Genes | Aligned Genes | Aligned bp | Best Blast Match | % Pairwise Identity |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| AM-TROP | Puerto Rico | - | - | Jatropha gossypifolia | 106,701,526 | 320 | 52 | 2192 | 2183 | 3,662,977 | KX397317 | 99.7 |
| AM-TROP | Ecuador | Guayas | Isla Puná | Datura stramonium | 80,297,211 | 241 | 39 | 2192 | 2183 | 3,669,262 | EU427728 | 98.6 |
| AM-TROP | United States | Arizona | Tucson | Gossypium hirsutum | 86,382,401 | 259 | 42 | 2192 | 2183 | 3,672,677 | AY521259 | 100 |
| NAF–MED–ME | Sudan | Al Jazirah | Wad Mani | Phaseolus vulgaris | 94,837,626 | 285 | 46 | 2193 | 2184 | 3,672,242 | FJ188567 | 99.0 |
| SSA | DR Congo | - | Lulimba | Manihot esculenta | 92,877,327 | 279 | 45 | 2193 | 2184 | 3,668,485 | AF344246 | 99.7 |
| SSA | Uganda | Kajera | Bukoba | Manihot esculenta | 133,240,668 | 400 | 65 | 2193 | 2184 | 3,668,433 | AM040604 | 99.7 |
| SSA | Tanzania | Dar es-Salam | Dar es-Salam | Manihot esculenta | 103,951,478 | 312 | 51 | 2193 | 2184 | 3,671,055 | JQ286457 | 99.7 |
| AS–PAC–AU | China | Hainan | Hainan | Gossypium hirsutum | 97,765,484 | 293 | 48 | 2193 | 2184 | 3,672,435 | HG918196 | 99.8 |
| ASIA | India | Ahmedabad | New Delhi | Nicotiana tabacum | 103,114,416 | 309 | 50 | 2193 | 2184 | 3,671,735 | KF751570 | 100 |
| B. afer | Australia | New South Wales | Sydney | Hardenbergia sp. | 86,662,872 | 277 | 45 | 1938 | 1929 | 3,002,554 | AJ784260 | 87.9 |
Colors of rows correspond to putative species suggested by ABGD analyses. A) Uncorrected pairwise distances among B. tabaci individuals based on the nuclear supermatrix of 2184 genes. B) Uncorrected pairwise distances among individuals based on COI data.
| (A) | Geographical Acronym | Specimen | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| 1 | B. afer | - | ||||||||||
| 2 | NAF-MED-ME | SRS1302125 | 0.0755 | - | ||||||||
| 3 | NAF-MED-ME | Sudan | 0.0757 | 0.0057 | - | |||||||
| 4 | AM-TROP | Arizona | 0.0805 | 0.0164 | 0.0138 | - | ||||||
| 5 | AM-TROP | Ecuador | 0.0806 | 0.0163 | 0.0138 | 0.0012 | - | |||||
| 6 | AM-TROP | Puerto Rico | 0.0808 | 0.0159 | 0.0136 | 0.0023 | 0.0023 | - | ||||
| 7 | SSA | DR Congo | 0.0808 | 0.0223 | 0.0199 | 0.0221 | 0.0221 | 0.0222 | - | |||
| 8 | SSA | Tanzania | 0.0830 | 0.0250 | 0.0225 | 0.0246 | 0.0247 | 0.0249 | 0.0051 | - | ||
| 9 | SSA | Uganda | 0.0808 | 0.0221 | 0.0198 | 0.0221 | 0.0221 | 0.0221 | 0.0031 | 0.0052 | - | |
| 10 | AS–PAC–AU | China | 0.0799 | 0.0161 | 0.0135 | 0.0152 | 0.0153 | 0.0154 | 0.0215 | 0.0241 | 0.0216 | - |
| 11 | ASIA | India | 0.0791 | 0.0150 | 0.0125 | 0.0143 | 0.0143 | 0.0145 | 0.0205 | 0.0231 | 0.0206 | 0.0093 |
| (B) | Geographical Acronym | Specimen | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| 1 | B. afer | - | ||||||||||
| 2 | NAF–MED–ME | SRS1302125 | 0.2200 | - | ||||||||
| 3 | NAF–MED–ME | Sudan | 0.2180 | 0.0637 | - | |||||||
| 4 | AM-TROP | Arizona | 0.2229 | 0.1439 | 0.1508 | - | ||||||
| 5 | AM-TROP | Ecuador | 0.2180 | 0.1412 | 0.1487 | 0.0171 | - | |||||
| 6 | AM-TROP | Puerto Rico | 0.2229 | 0.1426 | 0.1501 | 0.0178 | 0.0199 | - | ||||
| 7 | SSA | DR Congo | 0.2282 | 0.1594 | 0.1552 | 0.1588 | 0.1602 | 0.1581 | - | |||
| 8 | SSA | Tanzania | 0.2193 | 0.1489 | 0.1579 | 0.1531 | 0.1545 | 0.1511 | 0.0825 | - | ||
| 9 | SSA | Uganda | 0.2249 | 0.1524 | 0.1627 | 0.1592 | 0.1606 | 0.1572 | 0.0824 | 0.0131 | - | |
| 10 | AS–PAC–AU | China | 0.2312 | 0.1377 | 0.1412 | 0.1410 | 0.1348 | 0.1403 | 0.1594 | 0.1646 | 0.1681 | - |
| 11 | ASIA | India | 0.2260 | 0.1419 | 0.1467 | 0.1597 | 0.1576 | 0.1604 | 0.1695 | 0.1681 | 0.1736 | 0.1317 |
© 2019 by the authors.
Details
; Walden, Kimberly K O 4 ; Paredes-Montero, Jorge R 5 ; Waterhouse, Robert M 6
; Johnson, Kevin P 7 1 Department of Entomology, University of Illinois Urbana-Champaign, 505 S. Goodwin Ave., Urbana, IL 61801, USA; Illinois Natural History Survey, Prairie Research Institute, University of Illinois, Champaign, IL 61820, USA
2 School of Plant Sciences, University of Arizona, 1140 E. South Campus Dr., Tucson, AZ 85721, USA
3 Department of Entomology, Purdue University, 901 W. State St., West Lafayette, IN 47907, USA
4 Department of Entomology, University of Illinois Urbana-Champaign, 505 S. Goodwin Ave., Urbana, IL 61801, USA
5 School of Plant Sciences, University of Arizona, 1140 E. South Campus Dr., Tucson, AZ 85721, USA; Facultad de Ciencias de la Vida, Escuela Superior Politécnica del Litoral, ESPOL, Campus Gustavo Galindo Km 30.5 Vía Perimetral, P.O. Box 09-01-5863, Guayaquil, Ecuador
6 Department of Ecology and Evolution, University of Lausanne, and Swiss Institute of Bioinformatics, 1015 Lausanne, Switzerland
7 Illinois Natural History Survey, Prairie Research Institute, University of Illinois, Champaign, IL 61820, USA