Earlier than expected introductions of the Bemisia tabaci B mitotype in Brazil reveal an unprecedented, rapid invasion history

Abstract During 1991, in Brazil, the presence of the exotic Bemisia tabaci B mitotype was reported in São Paulo state. However, the duration from the time of initial introduction to population upsurges is not known. To investigate whether the 1991 B mitotype outbreaks in Brazil originated in São Paulo or from migrating populations from neighboring introduction sites, country‐wide field samples of B. tabaci archived from 1989–2005 collections were subjected to analysis of mitochondrial cytochrome oxidase I (mtCOI) and nuclear RNA‐binding protein 15 (RP‐15) sequences. The results of mtCOI sequence analysis identified all B. tabaci as the NAFME 8 haplotype of the B mitotype. Phylogenetic analyses of RP‐15 sequences revealed that the B mitotype was likely a hybrid between a B type parent related to a haplotype Ethiopian endemism (NAFME 1–3), and an unidentified parent from the North Africa‐Middle East (NAF‐ME) region. Results provide the first evidence that this widely invasive B mitotype has evolved from a previously undocumented hybridization event. Samples from Rio de Janeiro (1989) and Ceará state (1990), respectively, are the earliest known B mitotype records in Brazil. A simulated migration for the 1989 introduction predicted a dispersal rate of 200–500 km/year, indicating that the population was unlikely to have reached Ceará by 1990. Results implicated two independent introductions of the B mitotype in Brazil in 1989 and 1990, that together were predicted to have contributed to the complete invasion of Brazil in only 30 generations.

In Brazil, a 1994 report suggested that the introduction of the B mitotype (biotype or referred to by some as MEAM species) occurred in 1991 (Lourenção & Nagai, 1994). Despite the significance of the invasion, the identity of the suspect B mitotype has not been verified beyond a provisional identification, which was based on silverleaf symptoms observed in whitefly-infested squash plants (Cucurbita spp.) (Lourenção & Nagai, 1994). The development of silverleaf symptoms in squash plants within the Americas has been a reliable indicator of B-mitotype invasions Schuster et al., 1991;Yokomi et al., 1990). The silverleaf phenotype was first associated with whitefly-infested squash in Israel during the 1980s (Burger et al., 1983), suggesting that B mitotype was endemic there and possibly to other nearby locales. The proposed endemism of the B mitotype to northeast coastal Africa and the Middle East (NAF-ME) has been corroborated by analyses of protein polymorphisms, molecular markers (RAPDs, microsatellites), and DNA sequence data (mtCOI; nuclear genes) (Brown, 2010;Hadjistylli, 2010;Hadjistylli et al., 2016).
The invasion and colonization of Brazil (2,100,000 mi 2 ; excluding the Amazon rainforest) by the B mitotype of B. tabaci has not been well-studied, with estimates suggesting that it occurred in about 6 years post-introduction in 1991 .
However, there is great potential for the extensive colonization to have occurred earlier than 6 years owing to the arid subtropical climate there, particularly, during the dry seasons, and the widespread cultivation of agricultural crops that are reproductive hosts.
Although the first recognition of the B mitotype in Brazil was in Sao Paulo in 1991, the duration from the arrival time to population upsurges is not known. Characterizing whitefly populations collected before 1991 is crucial to determine sites of introduction, the timeframe from arrival times to explosion into invasiveness, and pathways of invasion in Brazil.
Computer simulations of the demography of invasive insect species can aid in identifying local ecological factors that were critical to the establishment of an invasive insect so that effective mitigation measures can be identified and implemented (Coutts et al., 2018;Excoffier et al., 2009;Lawson Handley et al., 2011;Sandlund et al., 2001). Models that consider the geographical and environmental complexity of landscapes, the geographic positions (sites) of the invasion(s), and life history parameters of the target insect species have been used for modeling short-and long-distance dispersal in temporal and spatial terms. For example, the SPLATCHE software (Currat et al., 2019) can facilitate modeling of the migration of the B mitotype, and is particularly useful for accounting for dispersal constraints in the context of unique environmental and geographic characteristics (Currat et al., 2019).
In this study, adult whitefly from an archived collection of B. tabaci, available in Brazil were analyzed using a mitochondrial and nuclear molecular marker. Whitefly collections from Brazil represented field samples collected immediately before and after the proposed widespread invasion of the B mitotype in 1991 (Lourenção & Nagai, 1994 (Table S1). Three other whitefly samples from the JK Brown laboratory collection were included for pairwise comparisons or to root the tree, a sample belonging to the A mitotype (New World species) from Ecuador, and two from Uganda (courtesy, Dr. James Legg, IITA Tanzania), one belonging to the Q mitotype, and the other identified as the 'Uganda-Sweet potato' outlier.
Total DNA was extracted and purified using an optimized pro- 1.4 mol/L NaCl, 0.2% 2-mercaptoethanol, and 2% hexadecyltrimethylammonium bromide) and ground with a micropestle. The pestle was washed with 600 µl of CTAB lysis buffer to keep all whitefly fragments into the buffer. Proteinase K was added at 0.005 mg/ml, and samples were incubated overnight at 55°C, followed by incubation at 65°C for 15 min to inactivate the Proteinase K. One volume of chloroform was added, and the mix was centrifugated (Eppendorf Model 5415R, Eppendorf, Hamburg, Germany) at 16,000 g for 3 min at 4°C for phase separation. One volume of 100% isopropanol and 40 µg glycogen were added to the supernatant, with incubation for 10 min at 4°C. Total nucleic acids were pelleted by centrifugation at 16,000 g at 4°C for 10 min. The pellet was washed with 70% ethanol, air-dried, and dissolved in 10 mmol/L Tris-HCl, pH 8.0.
The mtCOI-3′ fragment of ~1015 base pairs (bp) consisting of the 3′-end of the mtCOI and contiguous tRNAleu of 66 bp was amplified by PCR from samples collected in Brazil and from two of the eight B mitotype haplotypes of B. tabaci, EN8 (NAFME 8) and Eth3 (NAFME 3). The previously published NAFME 1-2 and 4-7 reference sequences (Paredes-Montero et al., 2021) were downloaded from the GenBank database (Table S1). The mtCOI fragment was ampli- The RP-15 fragment of 773 bp was amplified by PCR from B. tabaci samples collected in Brazil and from two of the eight haplotype references for B. tabaci, EN8 (NAFME 8), and Eth3 (NAFME 3).
The COI and RP-15 amplicons of the expected size were cloned into the TA cloning vector pGEM-T Easy (Promega, Madison, WI) by heat shock-mediated transformation of E. coli DH5α competent cells (Green & Sambrook, 2012). Colony-PCR amplification (Gussow & Clackson, 1989) was carried out using the M13F-5'-TGTAAAACGACGGCCAGT and M13R-5'-AGGAAACAGCTATGACCATG primers (Promega, Madison, WI) to identify clones containing the expected insert size. Cycling conditions were initial denaturation at 94°C for 10 min, 35 cycles at 94°C for 60 s, 53°C for 60 s and 72°C for 1 min, and final extension at 72°C for 10 min. Initially, one COI and two RP-15 clones per whitefly sample were subjected to bi-directional DNA Sanger sequencing.
Each cloned insert was sequenced bidirectionally using an ABI 3700 capillary sequencer at the Genomics Core Facility, The University of Arizona, Tucson, Arizona, USA (http://uagc.arizo na.edu).

| Pairwise distances and phylogenetic analysis
Sequences were aligned using MUSCLE v3.8.31 (Edgar, 2004), implemented in the Align Multiple Sequences tool available in Mesquite v2.75 (Maddison & Maddison, 2018). The mtCOI and RP-15 sequence alignments were trimmed to 1022 and 754 nt in length, respectively, and terminal gaps were treated as missing data. The best-fitting evolutionary models for the mtCOI and RP-15 sequences were estimated based on consensus of the Bayesian information criterion, Akaike information criterion, and maximum likelihood value using jModelTest v2.1.10 (Darriba et al., 2012).
The Hasegawa-Kishino-Yano model (Hasegawa et al., 1985) with gamma-distributed rate variation among sites, and the Kimura-2parameter model (Kimura, 1980) were identified as the best-fit evolutionary model for the mtCOI and nuclear RP-15 data, respectively.
The mtCOI and RP-15 sequences for the A mitotype (New World species) from Ecuador, and Q-like mitotype from Uganda (Table S1), were included for between-mitotype divergence comparisons. The mtCOI and RP-15 phylogenetic trees were rooted with the mtCOI and RP-15 sequence, respectively, of the "Uganda-Sweet potato" Bemisia spp., collected from sweet potato plants in Namulonge, Uganda.

| Dispersal assumptions
Simulations of population expansion were estimated using SPLATCHE v3 (Currat et al., 2019) in a two-dimensional steppingstone model defined on a matrix of ~27.5 × 10 4 demes (pixels or cells) covering Brazil and parts of surrounding South American countries. Each deme represented a surface of ~72 km 2 and exchanged migrants with its four neighbors depending on the pre-defined parameters: (1) carrying capacity, an estimate of the maximum number of whiteflies that could be sustained given predicted resources available in the geographic deme; (2) friction, corresponding to the ease by which whiteflies move from deme to deme; and (3) migration rate, which refers to the percentage of the population likely to spread to a neighboring deme. Maps of land cover and topographic information for Brazil (altitude) were obtained from the geographical information system repository (https://www.diva-gis.org/gdata). The altitude map was coded using distinct friction values for each cell using QGIS v 3.16.3 (http://qgis.osgeo.org). Maps were re-sized to a layer resolution of 0.08 and exported as ASCII raster layers using the conversion tool available in the "Raster" menu, and unnecessary rows in map headers were removed using a text editor. The "land cover" map consisted of one of 22 vegetation types based on the classification by the Global Land Cover 2000 Project (GLC20000) (https://forobs.jrc. ec.europa.eu/produ cts/glc20 00/glc20 00.php).
Carrying capacity and friction values were assigned to demes based on the quality of habitat and barriers to whitefly dispersal (e.g., altitudes of >600 m a.s.l.). For instance, the largest relative carrying capacities at 100%, were assigned to the following vegetation types: agriculture intensive, mosaic agriculture degraded forests, mosaic agriculture degraded vegetation. The lowest carrying capacity of 0% was assigned to micro-niches with a vegetation type unlikely to sustain whitefly populations, which are fresh water flooded forests, mangroves, barren bare soil, water bodies, and permanent snow ice.
The likelihood of whitefly migration rate was given as a friction value that ranged from 0 to 1, where 0 was assigned to demes with no impediments for whitefly migration, that is, migration is strong, and 1 corresponded to extreme friction, or when migration is not possible. Friction type was set to 2 to consider barriers to dispersal given by the vegetation type (land cover) and altitude. For the vegetation map, the same criteria used to define carrying capacities were used to assign friction values, with the lowest friction assigned to irrigated agro-ecosystems (Brown & Bird, 1992;, consistent with the habitat in which the B mitotype is predicted to support rapid expansion through migration and dispersal (Bethke et al., 1991;Butler et al., 1989;Naranjo et al., 2010), and where the highest friction corresponded to ecosystems where B. tabaci was "least likely" to establish, for example, in bodies of water or dense, humid tropical forested areas, and other ecosystems that are inconsistent with known B mitotype habitats. For altitude, extreme friction was set to values above 600 m a.s.l. and lowland areas were set to anticipate unimpeded flight and dispersal.
A migration rate of 6% of the population at each deme was assigned in accordance with the migratory proportion of flight morphs identified for the B. tabaci Arizona-B mitotype (Blackmer & Byrne, 1993;Byrne & Houck, 1990). In the simulation, westward dispersal was predicted for half of the migratory whitefly morphs, based on the predicted direction of flight (dispersal) in Brazil, predicted as, east to west (Gontijo et al., 2013), with the remaining proportion of migratory whiteflies evenly distributed among the other three possible directions of dispersal, excluding ecosystems where they were deemed unlikely to establish (Stansly & Naranjo, 2010 Several studies have measured basic demographic parameters under controlled environments (Carabali et al., 2005;Kakimoto et al., 2007;Ma, 2005;Musa & Ren, 2005;Tsai & Wang, 1996)

| Mitotype identification by mtCOI and RP-15 sequence analyses
Based on comparative mtCOI sequence analysis, the B. tabaci whiteflies collected from 1989 to 2005 in Brazil were identified as haplotype NAFME 8 of the B mitotype, most closely related to the Israel B-like reference (EN8, NAFME 8). Results of the mtCOI pairwise distance analysis indicated that the B. tabaci B-mitotype sequences from Brazil diverged from those of the Ethiopia 3 (Eth3, NAFME 3) and the NAFME 1-7 haplotypes by 1.8% and 0.5-1.9%, respectively ( Figure 1a,b).
Phylogenetic analysis indicated the Brazil mtCOI sequences grouped with haplotypes 6-8 of the NAFME clade, which clustered as a sister clade to the NAFME haplotypes 1-5 (Figure 1a).
The NAFME sister clades were monophyletic with the Q mitotype reference mtCOI sequence from Uganda, to form the NAF-MED-ME lineage (de Moya et al., 2019). As expected, the A mitotype and Uganda-Sweet potato reference mtCOI sequences, representing distinct cryptic species, grouped as outliers.
The oldest B mitotype record in the EMBRAPA whitefly collection studied herein was represented by a sample collected in 1989 in Rio de Janeiro, southeastern Brazil. The second-earliest record was represented by a sample collected in 1990 in Ceará state in northeastern Brazil on the Atlantic coast (Figure 1b). These results suggest that the outbreaks of B mitotype reported during 1991 in São Paulo city (São Paulo state) (Lourenção & Nagai, 1994) were likely associated with invasions associated with westward-dispersing B. tabaci originating from Rio de Janeiro where they were introduced on infested ornamental plants.
Based on phylogenetic analysis of the RP-15 nuclear sequences, females (2n) of the B mitotype analyzed from the 1989-2005 collections were found to harbor two different alleles, which corresponded to the two well-differentiated BZ-B1 and BZ-B2 subclades ( Figure 1c). In addition, one B. tabaci sample, which was represented by a single male (haploid) (40399), yielded one allele, inadvertently serving as an internal control. The B mitotype (female) endemic to Israel (EN8, NAFME 8) also harbored BZ-B1 and two alleles. In contrast, the nuclear RP-15 sequences determined for the reference Eth3 (NAFME 3) (female) from Ethiopia, grouped with the BZ-B2 subclade only.
The 1.5% intra-mitotypic divergence between the BZ-B1 and BZ-B2 haplotypes was similar to the inter-mitotypic divergence between the BZ-B1 and BZ-B2 and the A-mitotype (cryptic species) reference, at 1.5-2.2% pairwise distance, which was 0.1% less than the pairwise distances between the B and Q nuclear types, at 1.6-1.75%.
Collectively, these results suggest that the invasive B haplotype associated with the outbreaks in Brazil, identified as NAFME 8, and the reference EN8 (NAFME 8) from Israel, could possibly represent the offspring resulting from a potentially recent hybridization event.
Further, results indicate that one parent of hybrid NAFME 8 is most closely related to the B mitotype reference from Ethiopia (NAFME 1-3). The other parental allele(s) did not match any of the RP-15 alleles determined here and further exploration was beyond the scope of this study (Figure 1c,d).

| Simulations of demographic expansion of the B mitotype in Brazil
The demographic expansion predictions outward from Rio de Janeiro, the earliest record of the B mitotype presence, indicated that the B mitotype likely dispersed westwardly, unimpeded, at a speed rate of 200-500 km/year, reaching São Paulo in only 15 generations by 1990, but would have been unable in this same year to reach Ceará, the site of second oldest archived B mitotype sample (record) in Brazil (Figure 2). The low-friction migration micro-niches in the tropical semideciduous forest and the "cerrado" (savanna) ecoregions with extensive monocrop acreage appear to have facilitated the rapid spread of B mitotype populations from expansions emanating from Rio de Janeiro. These analyses predicted that outbreaks of the B mitotype reported in São Paulo during 1991 (Lourenção & Nagai, 1994) were linked to the westward-dispersing populations originating in Rio de Janeiro.
Simulations of the demographic expansion from Ceará, the second-earliest record of the B mitotype (this report), indicated that populations invading from this locale migrated at a slower rate than the those from Rio de Janeiro, at 100-300 km/year. The slower rate can likely be attributed to the climatic conditions there, which are Only nonredundant representative nuclear sequences were included in the tree. Colored whitefly cartoons show the position of reference samples EN8 (NAFME 8) and Eth3 (NAFME 3) in the phylogenetic tree. (d) Distribution of putative B. tabaci hybrids in Brazil, lower-left map shows estimated niche ranges of EN8 (NAFME 8) and Eth3 (NAFME 3) whitefly references in their center of origin. Node support corresponded to posterior probabilities. The sequences of a Uganda-Sweet potato B. tabaci mitotype was used to root phylogenetic trees F I G U R E 2 SPLATCHE simulations of demographic expansion of the Bemisia tabaci B mitotype in Brazil. Left-upper, -middle and -lower maps show the simulated expansion of the B mitotype from Rio de Janeiro after 15 generations using minimum, average and maximum demographic parameters, respectively. Right-upper, -middle and -lower maps show the simulated expansion of the B mitotype from Rio de Janeiro after 30 generations using minimum, average and maximum demographic parameters, respectively, and from Ceará after 15 generations

| Plant species associated with early B mitotype collections in Brazil
The plant species binomial names were corroborated using the bo-  (Brown, 2010;. The B mitotype sample archived in 1990 from Ceará state was collected from Spondias mombin, or hog plum, which is native to the tropical Americas, including the Eastern Caribbean region. Because Anacardiaceous species have not been known to serve as reproductive hosts of B. tabaci, it is most likely that the infestations of S. mombin were associated with populations dispersing from densely colonized preferred host plants including nearby senescing crops, even though reproduction on S. mombin was unlikely to occur.

| DISCUSS ION
Although the B-mitotype was first implicated in whitefly B. tabaci outbreaks in São Paulo, Brazil based on the silver-leaf symptoms observed there in pumpkin plants (Lourenção & Nagai, 1994) Further, the year of collection for the whitefly samples analyzed in this study are without question because they represent subsets of specimens archived in the EMBRAPA insect collection before the purported 1991 outbreak, which was not reported until three years later, in 1994 (Lourenção & Nagai, 1994). That the B mitotype was introduced into Brazil for the first time in 1989, is consistent with the earliest observations of silverleaf symptoms in squash throughout Central America and Caribbean region during 1988(Brown, 1994Brown & Bird, 1992;Costa & Brown, 1990), and only slightly later in South America (Viscarret et al., 2003). The two earliest detections of the B mitotype were predicted to be capable of invading and establishing throughout most of tropical-subtropical Brazil, which has an estimated surface area of 2,100,000 mi 2 (excluding the Amazon rainforest), in 2 years or less.
The rapid geographic expansion of B. tabaci from the two independent introductions, traceable to 1989 and 1990, respectively, led to the two migrating populations to coalesce and establish a nearly country-wide infestation by ~1991. By comparison, invasion of the US by the B mitotype, including the suspect initial sites of introduction first in Hawaii and then in Florida, required approximately 4-5 years (1986)(1987)(1988)(1989)(1990)(1991) (Bellows et al., 1994;Brown, Coats, et al., 1995;Gill & Brown, 2009), whereas, complete establishment of the B mitotype in North Queensland, with a surface area of 30,900 mi 2 , is estimated at 3-4 years (Liu et al., 2007). Regular surveillance of B.
tabaci in Brazil and elsewhere in the Tropical Americas would facilitate the early detection of other exotic B. tabaci haplotypes, while also inform pest management strategies to prevent crop loss due to damage caused by B. tabaci as a pest and vector of plant viruses.
Unexpectedly, analysis of the RP-15 nuclear marker revealed that the invasive haplotype of the B mitotype in Brazil and a reference population from Israel, EN8 (NAFME 8), harbored two alleles, BZ-B1 and BZ-B2. Further, only one of the two alleles, BZ-B2, was recovered from the Eth3 reference (NAFME 3) from Ethiopia, suggesting the B haplotype that invaded Brazil (and elsewhere) represents a potentially recent, intra-specific hybridization event. The NAFME region is the predicted center of origin of the B mitotype (Brown, 2010), and recently, the B mitotype has been shown to comprise at least eight phylogeographic haplotypes, several that inhabit overlapping niche ranges with potential as "hybrid zones" (Paredes-Montero et al., 2021). The B-mitotypes, EN8 (NAFME 8) and Eth3 (NAFME 3), are two examples of such haplotypes. Thus, in one hybridization scenario, the Eth3 (NAFME 3) or closely related NAFME 1-2, harboring only one of the two RP-15 alleles, served as one of the two parental haplotypes. The EN8 (NAFME 8) haplotype, which carries two different RP-15 alleles, could feasibly be the founder (hybrid) population of the B mitotype that has been introduced worldwide (Paredes- Montero et al., 2021). Although the other parent haplotype has not been identified among the reference samples analyzed here (albeit, a small sample size), the missing parent could have been displaced by its hybrid offspring (Abbott et al., 2016;Bellows et al., 1994;Chaturvedi et al., 2020;Epifanio & Philipp, 2000). This latter hypothesis offers the most parsimonious explanation, given the substantial knowledge of its invasiveness and extensive damage to food, fiber, and ornamental crops since the earliest introductions, which occurred during the 1980s. In addition, NAFME 8 is known to displace at least two mitotypes native to the Americas (New World A, Jatropha) and/ or exotic or endemic Q haplotype(s) (Brown, 2010;Brown, Coats, et al., 1995;Liu et al., 2007;Paredes-Montero et al., 2020;Zang et al., 2005). Finally, the co-existence of putative hybrids and parental taxa, including Eth3 (NAFME 3) or its closest relative (NAFME 1-2), could be explained based on the diverse climate niches they occupy, spanning the Arabian Peninsula to portions of coastal northeast Africa (Hadjistylli, 2010;Paredes-Montero et al., 2021).
The fixation of one mitochondrial type in the putative intraspecific hybrid, EN8 (NAFME 8), suggests the possibility of unidirectional gene flow as a potential driver of selection that favored the hybridization event. Unidirectional gene flow is well-known to be associated with increased competitive advantage of a dominant parent over a less dominant one, which can result in asymmetric mating (Caballero, 2007;Liu et al., 2007). Alternatively, unidirectional hybridization can be explained by cytoplasmic incompatibility (Nirgianaki et al., 2003;Li et al., 2007;Lv et al., 2021) phenotypes such as male-killing (Jiggins, 2003;Lv et al., 2021;Pan et al., 2012), and feminization of genetic males, which have been attributed to endosymbionts known to be harbored by whiteflies, including Rickettsia, Wolbachia, and Cardinium (Barro & Hart, 2000;Caballero, 2007;Giorgini et al., 2009). An expanded analysis by whole genome sequencing of extant B haplotypes could provide important new insights into the nature of this hybrid, NAFME 8, which has been ostensibly referred to as "superbug" (Kays et al., 1994;Markham et al., 1994;Perring, 2001;Wan et al., 2009).

ACK N OWLED G M ENTS
We thank Dr Norton Polo, Embrapa Genetic Resources and Biotechnology (CENARGEN), Brasília -Brazil, for kindly providing access to the Entomological Museum collection. This study was supported in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brazil (CAPES) -Finance Code 001. Special thanks to The University of Arizona, Federal University of Goiás -Brazil (UFG) and Embrapa for supporting this work.

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

DATA AVA I L A B I L I T Y S TAT E M E N T
The mtCOI and RP-15 sequences were submitted to GenBank, and accession numbers are available in Table S1, which was deposited to Dryad at the following URL: https://doi.org/10.5061/dryad.41ns1 rnfz.