Assessing the diversity of whiteflies infesting cassava in Brazil

Whitefly and cassava samples

Whiteflies were collected exclusively from cassava (M. esculenta) plants across 12 Brazilian states representative of the five macroregions (North, Northeast, Midwest, Southeast and South; Fig. 1) between March 2016 and February 2019 (Table 1). To gather evidence of whether a given species was colonizing cassava, adults and nymphs from the same field were collected whenever possible (Table 1). Samples were obtained from commercial and non-commercial (subsistence) crops. Whitefly adults were sampled using a hand-held aspirator and nymphs were collected with the aid of a needle. Insects were preserved in 95% ethanol and stored at −20 °C until being used for molecular identification of the species.

To verify the presence of begomoviruses infecting cassava, foliar samples were also collected at some sampled sites (Table S1). The samples were collected randomly regardless of the presence of virus-like symptoms. The leaves were press-dried and stored at room temperature as herbarium-like samples until being used for DNA extraction.

Whitefly species identification

Whitefly species were identified by sequencing of a mtCOI fragment, as previously described (Moraes et al., 2018). When enough adults and nymphs were collected at a given sampled site, ten individuals from each stage were analyzed, and when only one stage was obtained, 20 individuals were tested (Table 1). An initial assessment of whitely diversity was done using PCR-RFLP of the mtCOI fragment. When variation in the RFLP pattern was observed in the first screening, suggesting that more than one species could be present in that site, approximately five additional individuals for each stage were sequenced according to sample availability.

Total DNA was extracted from single individual whiteflies following a Chelex protocol (Walsh, Metzger & Higuchi, 1991). Briefly, adults or nymphs were ground in 30 µl of Chelex buffer (5% Chelex in 1×Tris-EDTA) using a toothpick in a 600 µl tube. Samples were vortexed for 30 s and incubated at 99 °C for 8 min in a PTC-100 thermocycler (MJ Research). Next, the tubes were centrifuged at 14,000 g for 5 min and 20 µl of the supernatant was collected and transferred to a new tube. One microliter of the supernatant was used as a template for PCR amplification of a 800 bp fragment of the mtCOI gene using primers C1-J-2195 and L2-N-3014 (Frohlich et al., 1999; Simon et al., 1994). PCR was performed using 0.2 µM of forward and reverse primers in a final volume of 25 µl using GoTaq Colorless Master Mix (Promega), following the manufacturer’s instructions. The PCR cycles consisted of an initial denaturing step at 95 °C for 5 min, followed by 35 cycles at 95 °C for 30 s, 42 °C for 45 s and 72 °C for 1 min, with a final extension at 72 °C for 10 min. Amplified products were visualized in 0.8% agarose gels stained with ethidium bromide and directly used for RFLP analysis (Bosco et al., 2006).

RFLP analysis of the amplicons consisted of 5 µl of each PCR product digested with 0.1 unit of Taq I (Promega) in a final volume of 20 µl. Reactions were performed at 65 °C for 2 h and visualized in 1.2% agarose gels stained with ethidium bromide. To verify whether the predicted mtCOI restriction pattern corresponded to a given species according to in silico prediction, a subset of PCR products from adults and nymphs representative of distinct patterns from different sampled sites were selected and sequenced. PCR products were precipitated with 100% ethanol and 3 M sodium acetate pH 5.2 (Sambrook & Russel, 2001) and sequenced commercially (Macrogen Inc.) in both directions using primers C1-J-2195/L2-N-3014.

For a small subset of samples that failed to yield a PCR product using primers C1-J-2195 and L2-N-3014, a second screening, using a recently described primer set with improved specificity for species of the B. tabaci complex and B. afer (2195Bt and C012/Bt-sh2), was performed (Mugerwa et al., 2018). Samples that still failed to amplify or had unexpected RFLP pattern were analyzed with specific primers for T. vaporariorum (TvapF and Wfrev) (Scott et al., 2007).

Sequence comparisons and phylogenetic analysis

Nucleotide sequences were first checked for quality and assembled using Geneious v. 8.1 (Kearse et al., 2012). mtCOI sequences were initially analyzed with the BLASTn algorithm (Altschul et al., 1990) to determine the whitefly species with which they shared greatest similarity. Pairwise comparisons between all mtCOI sequences obtained here and those with higher similarities (as determined by the BLASTn search) were performed with the program SDT v. 1.2 (Muhire, Varsani & Martin, 2014) using the MUSCLE alignment option (Edgar, 2004).

For phylogenetic analyses, the final dataset was composed of 142 sequences: 95 obtained in this work and 47 sequences representative of species in the family Aleyrodidae. Sequences were retrieved from GenBank and from the updated mtCOI reference dataset for species of the Bemisia tabaci complex (Boykin, Savill & De Barro, 2017). Multiple sequence alignments were prepared using the MUSCLE option in MEGA7 (Kumar, Stecher & Tamura, 2016). Alignments were checked and manually adjusted when necessary. Phylogenetic trees were constructed using Bayesian inference performed with MrBayes v. 3.0b4 (Ronquist & Huelsenbeck, 2003). The program MrModeltest v. 2.2 (Nylander, 2004) was used to select the nucleotide substitution model with the best fit in the Akaike Information Criterion (AIC). The analyses were carried out running 50,000,000 generations with sampling at every 1,000 generations and a burn-in of 25%. The convergence was assumed when average standard deviation of split frequencies was lower than 0.001. Trees were visualized and edited using FigTree (tree.bio.ed.ac.uk/software/figtree) and CorelDRAW X5, respectively.

Virus detection in foliar samples

Total DNA was extracted as described (Doyle & Doyle, 1987) and used as a template for PCR using the DNA-A universal primer pair PAL1v1978 and PAR1c496 (Rojas et al., 1993). PCR was performed in a final volume of 25 µl using Taq DNA Polymerase (Invitrogen) following the manufacturer’s instructions. The PCR cycles consisted of an initial denaturing step at 95 °C for 5 min, followed by 35 cycles at 95 °C for 1 min, 52 °C for 1 min and 72 °C for 1 min, with a final extension at 72 °C for 10 min. PCR products were visualized in 0.8% agarose gels stained with ethidium bromide. In addition, rolling-circle amplification (Inoue-Nagata et al., 2004) followed by digestion with MspI was performed in a subset of the samples.

Diversity index and statical analysis

Simpson’s index of diversity (1-D) was calculated to verify if there was any difference in whitefly diversity across macroregions. This index represents the probability that two randomly chosen individuals in a given sampled site will belong to distinct species (Morris et al., 2014). Simpson’s index was chosen as its value increases with increasing diversity and assigns more weight to more abundant species in a sample. We assume that species colonizing cassava will be in abundance, whereas rare species that briefly visit the plant without colonizing it will be underrepresented. Simpson’s index was calculated for each sampled site separately and then pooled according to macroregions. To assess the statistical significance of the differences in diversity among regions, the non-parametric Kruskal-Wallis test followed by post hoc multiple comparison test using Fisher’s least significant difference was calculated, using the function kruskal implemented in the Agricolae package in R software (R Core Team, 2017). Non-parametric Spearman’s rank correlation coefficient analysis was performed using the ggpubr package in R software (R Core Team, 2017).

High whitefly species richness in cassava in Brazil

To verify the composition of whitefly communities colonizing cassava in Brazil, sampling was performed across the country, including the main agroecological zones. A total of 66 sites from 12 states were sampled (Fig. 1; Table 1). Out of 1,385 individuals submitted to PCR-RFLP analysis, 58 adults and 37 nymphs from different locations and representing distinct restriction patterns were sequenced (Table S2). The combination of PCR-RFLP followed by sequencing showed reliability and consistency for species identification without misidentification due to incongruence between the two methods.

Based on pairwise comparisons and molecular phylogeny of the partial mtCOI gene, we identified the presence of at least seven species comprising the whitefly community in cassava (Fig. 2; Table 1). Among them, T. acaciae and B. tuberculata, both previously reported in this crop, were the most prevalent, representing over 75% of the analyzed individuals. In addition, based on the criterion of 3.5% divergence to differentiate species within the B. tabaci complex, three B. tabaci species were identified, with BtMED previously reported, and BtMEAM1 and BtNW identified for the first time in cassava fields in Brazil (Fig. 2; Table 1). The species BtMEAM1 represented 18% of the total individuals analyzed, followed by BtMED (1.6%) and BtNW (0.21%).

Furthermore, two putative new species were identified (Fig. 2), provisionally named whitefly new species 1 and 2 (WtNEW1 and WtNEW2). The WtNEW1 mtCOI sequence (KY249522) showed highest identity (80.65%) and clustered close to the T. acaciae clade, comprised of individuals reported here and three other previously reported sequences from cassava in Brazil (Fig. 2). For WtNEW2, two mtCOI sequences obtained from an adult (JX678666) and a nymph (DQ989531) shared 97.81% among them and showed highest identity with B. tabaci (adult: 82.11%; nymph: 81.68%) and clustered as a basal sister clade to the genus Bemisia (Fig. 2). Although whitefly taxonomy is predominantly based on puparial characters (Hodges & Evans, 2005) and there is no taxonomic criterion established based in mtCOI sequences for most of the groups, as has been proposed for the B. tabaci complex, the level of divergence between the two proposed new species with the closest species is similar to the level of divergence observed between species already described within the Aleyrodidae, as demonstrated in pairwise comparisons (Table S3) and phylogenetic analysis (Fig. 2). Nevertheless, further molecular and morphological characterization should be performed. Together, these results indicate the existence of a high whitefly species richness in cassava in Brazil.

Both the prevalence and the capacity to colonize cassava differ among species

Nymphs were collected for samples identified as T. acaciae, B. tuberculata, BtMEAM1 and the two new putative species (Fig. 3A), suggesting that these species may colonize cassava. Nymphs were not obtained at two sites where BtMED was prevalent (SP1 and SP12). Although it could be suggested that this species has the potential to colonize cassava due to the high prevalence of adults at these two sites, the lack of nymphs suggests otherwise. Moreover, at the sites PR4 and MT6, BtMEAM1 predominated among adults but 100% of the nymphs were B. tuberculata, suggesting that the predominance at one stage does not necessarily mean predominance in another stage. Indeed, correlation analysis between the number of adults and nymphs, performed for all sites where both stages were sampled, showed no significant correlation between them (Fig. S1). Further sampling in those sites or free-choice experiments are necessary to confirm the potential of BtMED to colonize cassava. Considering the whole sampling, we detected only three adults of BtNW, suggesting an inability of this specie to colonize cassava.

To verify if prevalence differs among species across distinct developmental stages, the data were separated according to stage and the proportions of individuals were compared for the three most abundant species (Figs. 3B, 3C). Considering the entire data set, T. acaciae was the prevalent species, followed by B. tuberculata and BtMEAM1 ($x_2^2=152.63$, P < 2.2 ×10⁻¹⁶). The same was true according to stage, either adults ($x_2^2=28.61$, P < 6.1 ×10⁻⁰⁷) or nymphs ($x_2^2=169.44$, P < 2.2 × 10⁻¹⁶; Fig. 3B). However, caution is needed to interpret these results as only adults were sampled at some sites where BtMEAM1 and B. tuberculata were prevalent (Fig. 3A), which could bias the analysis, causing an underestimation of the number of nymphs for those species. Therefore, we also analyzed the data considering only those sites where both adults and nymphs were obtained. Again, T. acaciae was the predominant species followed by B. tuberculata and BtMEAM1 considering either the entire data set ($x_2^2=258.61$, P < 2.2 ×10⁻¹⁶) or only nymphs ($x_2^2=164.47$, P < 2.2 ×10⁻¹⁶). When only adults were considered, T. acaciae was still predominant ($x_2^2=113.52$, P < 2.2 ×10⁻¹⁶) but no difference between B. tuberculata and BtMEAM1 was observed ($x_1^2=0.505$, P = 0.477; Fig. 3C). Moreover, it could be argued that samples from Minas Gerais (MG) were overrepresented in our sampling (Fig. 1C), which could also bias the results presented above due to the predominance of T. acaciae in this state (Fig. 3A). To test this possibility, we analyzed the data excluding the samples from MG. In this case, when both stages were considered, B. tuberculata was predominant ($x_2^2=62.09$, P = 3.3 ×10⁻¹⁴) but no difference between T. acaciae and BtMEAM1 was observed ($x_1^2=1.91$, P = 0.166). When each stage was considered separately, B. tuberculata was predominant followed by Bt MEAM1 and T. acaciae (adults: $x_2^2=43.94$, P = 2.9 ×10⁻¹⁰; nymphs: $x_2^2=84.19$, P < 2.2 ×10⁻¹⁶). Together, these results indicate that the potential to colonize cassava differs among species, which could be due either to lower preference for the plant or to differences in the competitive ability among species during cassava colonization. In addition, they reinforce the low efficiency of BtMEAM1 to colonize cassava.

Competitive interference does not explain the differences in prevalence

Interestingly, at least two species were detected co-occurring at 51% of the sampled sites (Fig. 3A). To verify the possibility of competition among T. acaciae, B. tuberculata and BtMEAM1 to explain the observed differences in prevalence (instead of differences in host preference), the competitive capacity of these three species was inferred based on the analysis of predominance at the sites where they occurred together. Initially, we verified if there were any differences in incidence, defined here as the number of sampled sites where at least one individual belonging to one of the three species was detected (Fig. 4A). The results demonstrate that there were no differences in incidence among them ($x_2^2=1.25$, P = 0.537; Fig. 4A). In addition, no differences were observed when the proportion of sites where whitefly species occurred alone or in different combinations was compared ($x_6^2=3.26$, P = 0.776; Fig. 4B). However, when we compared the occurrence between BtMEAM1 and non-B. tabaci species at the sites where they occur alone, the number of sites with non-B. tabaci species was higher ($x_1^2=6.53$, P = 0.011; Fig. 4B). Thus, the competitive capacity was inferred based on the proportion of individuals from each species at the fields where these species were detected co-occurring in different combinations (Fig. 4C). Interestingly, at the sites where BtMEAM1 and B. tuberculata were sampled together, B. tuberculata predominated over Bt MEAM1, suggesting higher competitive potential (Fig. 4C). For all other species combinations, no evidence of differences in competitive capacity were observed (Fig. 4C). Together, these results suggest that, rather than competition, lower host preference by BtMEAM1 explains its non-prevalence compared to T. acaciae and B. tuberculata, resulting in low colonization rate as indicated by the low number of BtMEAM1 nymphs detected in cassava (Fig. 3).

Composition and species diversity of whiteflies differ among Brazilian regions

The predominance of species composing the whitefly community across macroregions varied considerably. While T. acaciae predominated in the North, Southeast and Northeast, it was not detected in the Midwest (Fig. 5A). In addition, B. tuberculata was detected in all regions, and was prevalent in the South and Midwest. BtMEAM1, although not prevalent in any of the regions, was also detected in all regions. Although the number of species detected was higher in the Southeast, where six species out seven were detected, whitefly diversity was significantly higher in fields in the Northeast according to Simpson’s index of diversity (Fig. 5B), with no differences among the other four regions.

No begomoviruses detected infecting cassava

To verify the presence of begomoviruses infecting cassava, we analyzed leaves sampled in some of the fields where whiteflies were collected (Table S1). Based on PCR detection using universal primers for begomoviruses, all plants were negative. Although the possibility of false negatives cannot be completely discarded, it is unlikely since the PCR assay used primers known to detect CMBs (Berrie, Rybicki & Rey, 2001; Rojas et al., 1993; Zhou et al., 1997) and none of the plants displayed symptoms (Table S1). Nevertheless, as a confirmatory step, we performed RCA followed by digestion with MspI in five samples, all with negative results (Table S1).