# A landscape genetic analysis of important agricultural pest species in Tunisia: The whitefly Bemisia tabaci

Ahmed Ben Abdelkrim, Tarek Hattab, Hatem Fakhfakh, Mohamed Sadok Belkadhi, Faten Gorsane

► To cite this version:

Ahmed Ben Abdelkrim, Tarek Hattab, Hatem Fakhfakh, Mohamed Sadok Belkadhi, Faten Gorsane. A landscape genetic analysis of important agricultural pest species in Tunisia: The whitefly Bemisia tabaci. PLoS ONE, 2017, 12 (10), pp.e0185724. 10.1371/journal.pone.0185724. hal-01928434

HAL Id: hal-01928434

https://hal.umontpellier.fr/hal-01928434

Submitted on 26 May 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# A landscape genetic analysis of important agricultural pest species in Tunisia: The whitefly *Bemisia tabaci*

Ahmed Ben Abdelkrim<sup>1,2*</sup>, Tarek Hattab<sup>3,°</sup>, Hatem Fakhfakh<sup>1,4</sup>, Mohamed Sadok Belkadhi<sup>5</sup>, Faten Gorsane<sup>1,4*</sup>

¹ Laboratoire de Génétique Moléculaire, Immunologie et Biotechnologie. Faculté des Sciences de Tunis, Université Tunis El Manar, Tunis, Tunisie, ² Institut Jacques Monod, CNRS UMR 7592, Université Paris Diderot, Sorbonne Paris Cité, Paris, France, ³ Institut Français de Recherche pour l'Exploitation de la Mer, IFREMER, UMR 248 MARBEC, Avenue Jean Monnet CS, Sète, France, ⁴ Faculté des Sciences de Bizerte, Zarzouna, Université de Carthage, Bizerte, Tunisie, ⁵ Centre Technique des cultures protégées et géothermiques de Gabes, Gabes, Tunisie

© These authors contributed equally to this work.

* fatengorsane@gmail.com, faten.gorsane@fss.rmu.tn

<figure><img src="image_2.png" /></figure>

## Abstract

Combining landscape ecology and genetics provides an excellent framework to appreciate pest population dynamics and dispersal. The genetic architectures of many species are always shaped by environmental constraints. Because little is known about the ecological and genetic traits of Tunisian whitefly populations, the main objective of this work is to highlight patterns of biodiversity, genetic structure and migration routes of this pest. We used nuclear microsatellite loci to analyze *B. tabaci* populations collected from various agricultural areas across the country and we determine their biotype status. Molecular data were subsequently interpreted in an ecological context supplied from a species distribution model to infer habitat suitability and hereafter the potential connection paths between sampling localities. An analysis of landscape resistance to *B. tabaci* genetic flow was thus applied to take into account habitat suitability, genetic relatedness and functional connectivity of habitats within a varied landscape matrix. We shed light on the occurrence of three geographically delineated genetic groups with high levels of genetic differentiation within each of them. Potential migration corridors of this pest were then established providing significant advances toward the understanding of genetic features and the dynamic dispersal of this pest. This study supports the hypothesis of a long-distance dispersal of *B. tabaci* followed by infrequent long-term isolations. The Inference of population sources and colonization routes is critical for the design and implementation of accurate management strategies against this pest.

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Ben Abdelkrim A, Hattab T, Fakhfakh H, Belkadhi MS, Gorsane F (2017) A landscape genetic analysis of important agricultural pest species in Tunisia: The whitefly *Bemisia tabaci*. PLoS ONE 12(10): e0185724. https://doi.org/10.1371/journal.pone.0185724

**Editor:** Maohua Chen, Northwest A&F University, CHINA

**Received:** June 27, 2017

**Accepted:** September 18, 2017

**Published:** October 3, 2017

**Copyright:** © 2017 Ben Abdelkrim et al. This is an open access article distributed under the terms of the [Creative Commons Attribution License](https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** This study was partially supported by the Ministry of Higher Education and Scientific Research of Tunisia (LR 99 ES 12) and the USA MERC (Middle East Research Coorporation) project TA-MOU-08-M28-048.

**Competing interests:** The authors have declared that no competing interests exist.
Introduction

The sweet potato whitefly, *Bemisia tabaci*, is an economically important agriculture pest con-
sisting of at least 34 cryptic members [1]. It is highly polyphagous and invasive, colonizing
more than 1000 different plant species and causing significant losses by feeding and through
acting as a vector for more than 300 plant viruses [2,3]. Based on mitochondrial cytochrome
oxidase I (COI) DNA retrieved from sequences of worldwide *B. tabaci*, 34 genetic groups have
been discriminated [4]. They belong to four major clusters: Sub-Saharan Africa; Asia; the New
World and the latest cluster that includes North Africa, the Middle East and Asia Minor [1].
Currently, the global status of this pest refers to two distinct species corresponding to Middle
East-Asia Minor 1 (MEAM1, formerly biotype B) and Mediterranean (MED, formerly biotype
Q) [1,5]. *B. tabaci* is now considered as a complex of well-defined groups that are referred to
biotypes and which are distinguishable according to host specialization, reproductive compati-
bility, differential resistance to different classes of insecticide and efficiency in transmitting
phytoviruses [4,6].

Both direct observation and model predictions show that *B. tabaci* is expanding its distribu-
tion in the Mediterranean basin [7]. Temperature is the major environmental variable that
influences the development, survival and reproduction of *B. tabaci* populations, and so their
potential geographic range [7,8]. Host plant availability and therefore agricultural and land use
practices, are also important factors that modulate the population distributions and structure
of *B. tabaci* across the landscape [9,10]. Understanding how climate and landscape heterogene-
ity shape the distribution of pest populations and gene flow between landscape patches is cru-
cial in pest management science [11,12]. Answer to this question requires the use of the
landscape genetics approach [13–15], which combines concepts and tools from population
genetics, landscape ecology, geography and spatial statistics in order to assess how landscape
resistance modulates animal species move across the environment [12]. Thus, landscape resis-
tance or its inverse landscape permeability can be modeled to (i) identify landscape and envi-
ronmental features that constrain genetic connectivity, (ii) and to perform a spatially explicit
exploration of genetic diversity and population structure, mainly to inform resource manage-
ment and conservation [12,16,17].

The population genetics of many pest species are shaped by environmental gradients
[18,19]. A combination of genetic diversity information based on molecular markers and envi-
ronmental approaches therefore has the potential to provide a powerful framework for the
study *B. tabaci* population dynamics and dispersal in Tunisia. Previously, it was reported that
the high divergence between the mitochondrial cytochrome oxidase I (mtCOI) sequences
among *B. tabaci* clades suggested that ecological factors, particularly climatic and tectonic
events, have contributed to the extreme diversification of *B. tabaci* [4,6,20]. In Tunisia, the
identification and biotype distribution of naturally occurring *B. tabaci* population were
assessed a few years ago [21]. In depth, analysis of the genetic affiliation of these populations
by PCR-RFLP (TaqI) of the mtCOI gene and sequencing indicated that the sampled whiteflies
clustered into the B and Q biotypes. Genetic investigation driven by several approaches includ-
ing molecular variance (AMOVA), selective neutrality and genetic haplotype network tests
supported the hypothesis that Tunisian *B. tabaci* populations underwent a potential expansion
followed by gene flow restriction [21].

Although mtCOI markers have been extensively used to provide insights into *B. tabaci*
groups and thus, could be used as “integrated barcodes”, they do not answer the same ques-
tions as integrated nuclear, morphological, and ecological data sets [22,23]. Such data sets
could be useful for confirming hypotheses based on other sources of data. mtDNA markers
also have limitation such as maternal inheritance, recombination, inconsistent mutation rates,
heteroplasmy and compounding evolutionary processes [22,24]. In contrast, nuclear molecular markers such as simple sequence repeats (SSRs) can be considered as more valuable tools for assessing the genetic structure and traits of natural *B. tabaci* populations. Out of the 50 described loci [9,25–29], a panel of seven microsatellite loci, widely used in *B. tabaci* genetic studies [28], were selected in this work and scored in naturally occurring populations located in different geographical parts of Tunisia. Results gave further insights into population structure and differential patterns of genetic diversity across defined geographical groups.

Molecular data are also important and useful in conjunction with other sources of information such as ecology data [22,30]. By merging SSR molecular markers with the biotype status of *B. tacaci* population and species distribution modeling (SDM) [31], we were also able to analyze potential migration corridors among *B. tabaci* populations in Tunisia. We used large-scale SDM to infer the climatic suitability of the landscape for *B. tabaci* and pathways based on least-cost distance between individuals with a high level of genetic relatedness. Using this approach, we mapped landscape resistance to *B. tabaci* genetic flow and we identified the landscape elements, such as natural barriers, that isolate *B. tabaci* populations and the potential connection corridors between sampling localities. These corridors reflected the dual sources of the connection: both genetic distance and environmentally suitable paths of migration inferred from the SDM. Thus, modeling landscape resistance and potential migration corridors represents a major advance for the prioritizing *B. tabaci* elimination campaigns, at the national level.

# Materials and methods

## *Bemisia tabaci* samples

A total of 1475 female *B. tabaci* adults corresponding to 59 distinct populations were collected during 2014–2016 from both open fields and greenhouses (Fig 1). Populations were sampled from 16 localities within the main vegetable-growing areas across Tunisia (Table 1). Whitefly samples were collected alive and immediately preserved in 96% ethanol at 20°C until DNA extraction was performed.

## DNA extraction

Twenty-five adult insects corresponding to one population were frozen in liquid nitrogen before homogenization in an extraction buffer (100 g/l Proteinase K; 0.45% Triton; 0.45% Tween; 1M Tris-HCl, pH: 8). The mixture was incubated for 60 min at 55°C, then at 100°C for 10 min, and finally at 0°C for 5 min. A centrifugation step was performed at 1,000 rpm for 5 min to pellet debris. DNA quantification was performed with ND-1000 spectrophotometer (Nanodrop Technologies, USA).

## Identification of *B. tabaci* biotypes

The Biotype status of *B. tabaci* samples was determined through (i) The PCR-RFLP (*TaqI*) of a fragment of the mitochondrial gene (mtCOI) [32] and (ii) The PCR amplification of a microsatellite fragment BEM23-(GAA)₃₁ imp [33]. For the first method, PCR reaction was performed in a final volume of 25 µl containing 100 ng of each primer [34], 5 µl of 10X buffer, 1 µl of dNTP (10 mM each), 1.5 µl MgCl2 (1.5 mM), 1U of *Taq* DNA polymerase (MP BIOMEDICALS), 2 µl of extracted DNA, and sterile water. PCR amplifications were performed on a T professional trio system thermalcycler (Biometra, Germany) using the following cycling program: One first denaturation step of 4 min at 94°C followed by 35 amplification cycles (30 s at 94°C, 50 s at 55°C, 1 min at 72°C). A final step of extension is performed at 72°C for 10 min.
<figure><img src="image_4.png" /><figcaption>**Fig 1.** Climatic suitability map of *Bemisia tabaci*. (A) across the geographical extent used for model calibration and (B) across Tunisia. Blue dots represent recorded occurrences.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0185724.g001

PCR products were purified (PureLink PCR Purification Kit, Invitrogen, Paris, France) and digested with 1 U of *TaqI* restriction enzyme with 2 µl of the corresponding buffer (10X) in a final volume of 20 µl, for 2 h at 65°C. PCR products or cleaved fragments were separated in 2% agarose gel and visualized by staining with ethidium bromide. For the second approach, microsatellite marker amplifications were carried out using primers developed by [33] and validated by [35, 36]. PCR-amplified alleles were scored as described below.

## Simple Sequence Repeat (SSR) genotyping

Seven microsatellite markers were selected based on their previous use [9]. Amplifications were carried out using the PCR conditions described by [28]. A green microfluidic chip was used following the DNA 1K analysis kit protocol. For each microsatellite marker, PCR-amplified bands were identified and scored by the Experion System (Automated Electrophoresis Station). Electrophoregrams were analyzed using the Experion software, System Operation and Data Analysis Tools (version 1.0).

## Genetic analysis and population structure

To evaluate the genetic diversity of *B. tabaci* populations, the basic genetic population parameters such as the number of alleles detected per locus (Na), the polymorphism information content (PIC) values, the observed (Ho) and the expected heterozygosity (He) as well as significance values for deviations from the Hardy-Weinberg equilibrium (HWE) were determined.

Multiple co-inertia analysis (MCOA) approach was performed to identify alleles characterizing geographic groups and to evaluate the efficiency of each marker to build population structure since this technique allows each locus to be analyzed separately [37]. Then, the proportion of variance attributable to genetic differences between and among the main geographical regions (North, Sahel and South) was estimated, for each biotype, using hierarchical analysis of molecular variance (AMOVA). To cluster sampled populations into genetic groups we used a discriminant analysis of principal components (DAPC). To achieve that, principal component analysis (PCA) was used to transform the raw SSR data into a set of linearly uncorrelated variables. Afterwards, we identified the optimal number of clusters that minimize the
<table><caption>Table 1. Collections' details for the 59 samples of *Bemisia tabaci* from Tunisia.</caption><thead><tr><td>Sample</td><td>Locatity</td><td>GPS coordinates</td><td>Geographic groups</td><td>Host plant</td><td>Biotypes</td></tr></thead><tbody><tr><td>Pop1</td><td>Teboulba</td><td>35°38'N, 10°56'E</td><td>Sahel</td><td>Tomato</td><td>B</td></tr><tr><td>Pop2</td><td>Mansoura</td><td>36°51'N, 11°7'E</td><td>North</td><td>Melon</td><td>Q</td></tr><tr><td>Pop3</td><td>Mansoura</td><td>36°51'N, 11°7'E</td><td>North</td><td>Cucumber</td><td>Q</td></tr><tr><td>Pop4</td><td>SidiNeji</td><td>36°51'N, 11°7'E</td><td>North</td><td>Cucumber</td><td>Q</td></tr><tr><td>Pop5</td><td>Manouba</td><td>36°48'N, 10°6'E</td><td>North</td><td>Melon</td><td>Q</td></tr><tr><td>Pop6</td><td>Bouficha</td><td>36°17'N, 10°27'E</td><td>North</td><td>Tomato</td><td>Q</td></tr><tr><td>Pop7</td><td>Bouficha</td><td>36°18'N, 10°27'E</td><td>North</td><td>Melon</td><td>Q</td></tr><tr><td>Pop8</td><td>RasDjbal</td><td>37°12'N, 10°7'E</td><td>North</td><td>Zucchini</td><td>Q</td></tr><tr><td>Pop9</td><td>Ghar El Melh</td><td>37°10'N, 10°10'E</td><td>North</td><td>Zucchini</td><td>Q</td></tr><tr><td>Pop10</td><td>Ghar El Melh</td><td>37°10'N, 10°10'E</td><td>North</td><td>Melon</td><td>Q</td></tr><tr><td>Pop11</td><td>Ghar El Melh</td><td>37°10'N, 10°11'E</td><td>North</td><td>Melon</td><td>Q</td></tr><tr><td>Pop12</td><td>El Sahline</td><td>35°44'N, 10°43'E</td><td>Sahel</td><td>Pepper</td><td>B</td></tr><tr><td>Pop13</td><td>Slimane</td><td>36°41'N, 10°28'E</td><td>North</td><td>Cauliflower</td><td>B</td></tr><tr><td>Pop14</td><td>Slimane</td><td>36°41'N, 10°28'E</td><td>North</td><td>Potato</td><td>B</td></tr><tr><td>Pop15</td><td>Slimane</td><td>36°41'N, 10°28'E</td><td>North</td><td>Tomato</td><td>B</td></tr><tr><td>Pop16</td><td>Slimane</td><td>36°42'N, 10°29'E</td><td>North</td><td>Zucchini</td><td>B</td></tr><tr><td>Pop17</td><td>Slimane</td><td>36°42'N, 10°28'E</td><td>North</td><td>Cauliflower</td><td>B</td></tr><tr><td>Pop18</td><td>Monastir</td><td>35°40'N, 10°45'E</td><td>Sahel</td><td>Tomato</td><td>B</td></tr><tr><td>Pop19</td><td>Monastir</td><td>35°46'N, 10°49'E</td><td>Sahel</td><td>Zucchini</td><td>B</td></tr><tr><td>Pop20</td><td>Kébili</td><td>33°42'N, 8°57'E</td><td>South</td><td>Lantana</td><td>Q</td></tr><tr><td>Pop21</td><td>Gabes</td><td>33°53'N, 9°47'E</td><td>South</td><td>Tomato</td><td>B</td></tr><tr><td>Pop22</td><td>Rjim-Maatoug</td><td>33°19'N, 8°0'E</td><td>South</td><td>Black Nightshade</td><td>Q</td></tr><tr><td>Pop23</td><td>Gabes</td><td>33°53'N, 9°47'E</td><td>South</td><td>Zucchini</td><td>B</td></tr><tr><td>Pop24</td><td>Kébili</td><td>33°42'N, 8°57'E</td><td>South</td><td>Tomato</td><td>Q</td></tr><tr><td>Pop25</td><td>Teboulba</td><td>35°38'N, 10°56'E</td><td>Sahel</td><td>Tomato</td><td>B</td></tr><tr><td>Pop26</td><td>SidiNeji</td><td>36°51'N, 11°7'E</td><td>North</td><td>Cucumber</td><td>B</td></tr><tr><td>Pop27</td><td>Bouficha</td><td>36°17'N, 10°27'E</td><td>North</td><td>Melon</td><td>Q</td></tr><tr><td>Pop28</td><td>RasDjbal</td><td>37°12'N, 10°7'E</td><td>North</td><td>Zucchini</td><td>Q</td></tr><tr><td>Pop30</td><td>El Sahline</td><td>35°44'N, 10°43'E</td><td>Sahel</td><td>Pepper</td><td>B</td></tr><tr><td>Pop31</td><td>Slimane</td><td>36°41'N, 10°28'E</td><td>North</td><td>Cauliflower</td><td>B</td></tr><tr><td>Pop32</td><td>Slimane</td><td>36°41'N, 10°28'E</td><td>North</td><td>Potato</td><td>B</td></tr><tr><td>Pop33</td><td>Slimane</td><td>36°41'N, 10°28'E</td><td>North</td><td>Zucchini</td><td>B</td></tr><tr><td>Pop34</td><td>Slimane</td><td>36°41'N, 10°28'E</td><td>North</td><td>Cauliflower</td><td>B</td></tr><tr><td>Pop35</td><td>Monastir</td><td>35°46'N, 10°49'E</td><td>Sahel</td><td>Tomato</td><td>B</td></tr><tr><td>Pop36</td><td>Monastir</td><td>35°46'N, 10°49'E</td><td>Sahel</td><td>Zucchini</td><td>B</td></tr><tr><td>Pop37</td><td>Kébili</td><td>33°42'N, 8°57'E</td><td>South</td><td>Lantana</td><td>Q</td></tr><tr><td>Pop38</td><td>Gabes</td><td>33°53'N, 9°47'E</td><td>South</td><td>Tomato</td><td>B</td></tr><tr><td>Pop39</td><td>Rjim-Maatoug</td><td>33°19'N, 8°0'E</td><td>South</td><td>Black Nightshade</td><td>Q</td></tr><tr><td>Pop40</td><td>Gabes</td><td>33°53'N, 9°47'E</td><td>South</td><td>Zucchini</td><td>B</td></tr><tr><td>Pop41</td><td>Kébili</td><td>33°42'N, 8°57'E</td><td>South</td><td>Tomato</td><td>Q</td></tr><tr><td>Pop42</td><td>Teboulba</td><td>35°38'N, 10°56'E</td><td>Sahel</td><td>Tomato</td><td>B</td></tr><tr><td>Pop43</td><td>SidiNeji</td><td>36°51'N, 11°7'E</td><td>North</td><td>Cucumber</td><td>B</td></tr><tr><td>Pop44</td><td>Bouficha</td><td>36°17'N, 10°27'E</td><td>North</td><td>Melon</td><td>Q</td></tr><tr><td>Pop45</td><td>RasDjbal</td><td>37°12'N, 10°7'E</td><td>North</td><td>Zucchini</td><td>Q</td></tr><tr><td>Pop46</td><td>Ghar El Melh</td><td>37°10'N, 10°11'E</td><td>North</td><td>Zucchini</td><td>Q</td></tr><tr><td>Pop47</td><td>El Sahline</td><td>35°44'N, 10°43'E</td><td>Sahel</td><td>Pepper</td><td>B</td></tr><tr><td>Pop48</td><td>Slimane</td><td>36°41'N, 10°28'E</td><td>North</td><td>Cauliflower</td><td>B</td></tr></tbody></table>

(Continued)
<table><caption>**Table 1.** (Continued)</caption><thead><tr><th>Sample</th><th>Locatity</th><th>GPS coordinates</th><th>Geographic groups</th><th>Host plant</th><th>Biotypes</th></tr></thead><tbody><tr><td>Pop49</td><td>Slimane</td><td>36°41'N, 10°28'E</td><td>North</td><td>Potato</td><td>B</td></tr><tr><td>Pop50</td><td>Mansoura</td><td>36°51'N, 11° 7'E</td><td>North</td><td>Melon</td><td>B</td></tr><tr><td>Pop51</td><td>Kébili</td><td>33°20'N, 8°18'E</td><td>South</td><td>Melon</td><td>Q</td></tr><tr><td>Pop52</td><td>Kébili</td><td>33°34'N, 9° 1'E</td><td>South</td><td>Melon</td><td>Q</td></tr><tr><td>Pop53</td><td>Kébili</td><td>33°42'N, 8°58'E</td><td>South</td><td>Lantana</td><td>Q</td></tr><tr><td>Pop54</td><td>Gabes</td><td>33°52'N, 10° 5'E</td><td>South</td><td>Zucchini</td><td>B</td></tr><tr><td>Pop55</td><td>Kébili</td><td>33°44'N, 9°14'E</td><td>South</td><td>Melon</td><td>Q</td></tr><tr><td>Pop56</td><td>Kébili</td><td>33°44'N, 9°14'E</td><td>South</td><td>Tomato</td><td>Q</td></tr><tr><td>Pop57</td><td>Kébili</td><td>33°39'N, 9° 0'E</td><td>South</td><td>Melon</td><td>Q</td></tr><tr><td>Pop58</td><td>Gabes</td><td>33°53'N, 9°47'E</td><td>South</td><td>Tomato</td><td>B</td></tr><tr><td>Pop59</td><td>Gabes</td><td>33°52'N, 10° 5'E</td><td>South</td><td>Tomato</td><td>B</td></tr><tr><td>Pop60</td><td>Kébili</td><td>33°42'N, 8°58'E</td><td>South</td><td>Melon</td><td>Q</td></tr></tbody></table>

Pop: population

**Solanaceae:** Tomato (*Solanum lycopersicum*.), Eggplant (*Solanum melongena*), pepper (*Capsicum annuum*), Potato (*Solanum tuberosum*), Black Nightshade (*Solanum nigrum*); **Cucurbitaceae:** Melon (*Cucumismelo*), Zucchini (*Cucurbita pepo*), **Brassicaceae:** Cauliflower (*Brassica oleracea*) and **Verbenacea** (*Lantana camara*).

https://doi.org/10.1371/journal.pone.0185724.t001

variation within clusters using the Bayesian Information Criterion (BIC). A discriminant analysis was then applied to 14 principal components, explaining 59.4% of the total variance of the SSR data. Based on the retained discriminant functions, a membership probability of each population to each cluster was calculated, which can be interpreted in order to assess how clear-cut or admixed the clusters are [38]. Visualization of the relationships between population and clusters in the DAPC space were realized using a scatterplot. For the visualization of relationships between genetics groups, we calculated Nei's genetic distances [39] between groups and we constructed a neighbor-joining tree.

## Species distribution modeling

Ecological niche factor analysis (ENFA) [40] was used to create a climatic suitability map that depicts areas where *B. tabaci* is unlikely to occur. ENFA is a specific multivariate ordination technique that relies on presence-only data to compare a species' environmental niche and the environmental characteristics of the studied area. To this end, ENFA compares the distribution of the environmental predictors between the occurrence locations and the whole studied area and assigns a degree of suitability to each point on a map (typically from 0 to 100). We used as environmental predictors 19 bioclimatic variable data provided by WorldClim database (version 1.4) [41]. This bioclimatic dataset is derived from monthly temperature and rainfall value data for the period 1950–2000 and represents annual trends, seasonality and extreme or limiting environmental factors. In particular, we used bioclimatic surfaces with pixel size of 30 arc seconds corresponding to a resolution of approximately 1 km and covering North Africa, the Middle East and part of Asia (Fig 1). This geographical extent was chosen for two reasons: (1) it matches the range of distribution of one of the high-level genetic groups of *B. tabaci* [1] (2) it is advised to develop SDMs across the range of climatic conditions in which a given species occurs [42,43].

To calibrate the ENFA model, we merged our spatial genetic dataset with 207 geographical occurrences located within the bioclimatic surface areas and obtained from the Global Biodiversity Information Facility (GBIF; http://data.gbif.org, 2016) (Fig 1). To evaluate the accuracy
of the ENFA, we performed a ten-fold cross validation based on the Boyce index, which measures how far model predictions differ from random distribution of the observed presences across prediction gradients [44]. It is continuous and varies between -1 and +1. Positive values indicate a model whose predictions are consistent with the presence distribution in the evaluation dataset. Values close to zero mean that the model is no different from a chance model, while negative values indicate an incorrect model, which predicts poor quality areas where presences are frequent [45].

## Predicting population connection corridors

Climatic suitability cells located in Tunisia were first taken from the ENFA prediction map. Given the strong link between *B. tabaci* and the great number of host species for this phytopathous insect [46], we assigned zero values of habitat suitability to cells located outside vegetation areas in Tunisia. Land cover data were obtained from the GLCN Land Cover datasets (http://www.glcn.org/databases/lc_gc-africa_en.jsp) at 300-m resolution. The resulting map, reflecting both abiotic (i.e. temperature and rainfall) and biotic (host and/or prey availability) constraints, was used to create a resistance surface. To this end, a transition matrix was computed from the habitat suitability values to define the connectedness between adjacent pixels. We used King's graph as a neighborhood function, in which a given pixel is considered to be connected to the eight adjacent pixels. Each cell on a resistance surface is assigned a cost value, with high costs given to unsuitable cells that restrict dispersal and low costs to suitable cells that facilitate dispersal [47]. Thus, higher resistance values represent specific factors that limit the ability of an individual *B. tabaci* to traverse a cell, which could be in this case a function of mortality risk, energetic costs, food availability and interactions between each of these factors.

Based on a Wright's *F*<sub>ST</sub> paire-wise genetic distance matrix, we extracted likely genetically identical samples: pairs with Wright's *F*<sub>ST</sub> index [48] lower than percentile 5% of all pair-wise comparisons. Then, using the resistance surface we calculated accumulated cost surfaces between each pair of these genetically related samples using Dijkstra's algorithm[49,50]. The accumulated cost surfaces are a measure of nearness that optimizes the geographical distance travelled and the costs traversed [50]. The resulting accumulated cost surfaces were then normalized between 0 and 1 to ensure comparability between the different resistance sets. Finally, all corridors were merged in one map representing the network of potential connection corridors between sampled populations.

All data processing and statistical analysis were performed using R 3.1 [51].

## Results

### Genetic analysis

We analyzed 59 *B. tabaci* populations collected on different hosts throughout the country and clustered them into three major geographic groups: North, Sahel and South (Table 1). *B. tabaci* populations were assigned to either biotype B or Q according to PCR-RFLP (*TaqI*) patterns of the mtCOI gene fragment in conjunction with a diagnostic microsatellite marker. Two restriction patterns (414, 162, 144 bp) and (631,144 bp) corresponding to biotype B and Q respectively, were revealed (S1 Fig). Results were corroborated by Bem23 microsatellite analysis since the alleles > 371 bp were scored only within Q biotype. Biotype B was detected in the three geographic groups, all of which were on vegetables. Biotype Q detection was restricted to the North and the South of the country, on both vegetables and ornamentals (Table 1)

In addition, seven loci were selected to analyse populations' structure. These gave clear patterns of discrimination between the sampled geographical groups. A total of 113 alleles were scored across the studied populations. The number of alleles (Na) per locus ranged from 6
(detected in locus 177 in the Sahel) to 18 (detected in locus BT4 in the North). The North displayed the highest average Na value (13.42) while the Sahel showed the lowest (8.00). The observed (Ho) and expected (He) heterozygoties ranged from 0.001 to 1.00 and from 0.74 to 0.92 respectively. Regarding polymorphism information content (PIC), all considered loci displayed high values, ranging from 0.82 to 0.92, attesting to the effectiveness of these SSR markers for detecting variability (Table 2).

Out of the 113 revealed alleles, 19, 57 and 29 were not detected in the North, the Sahel and the South, respectively. Subsequently, we outlined the most group-specific alleles for each region corresponding to 23 alleles in the North, 11 alleles in the South and only 3 specific alleles in the Sahel (Table 3).

## Genetic relationships and population structure analysis

We next addressed the capacity and efficiency of marker panels to exhibit genetic structure and to assess the contribution of each specific marker by MCOA according to [37]. This approach permits the identification of alleles of interest that contribute most to the genetic topology of *B. tabaci* populations, based on allelic frequencies and single-marker analysis (Fig 2). For example, with regard to the BT4 locus, the first axis accounts for 68.84% of total variance while the second axis accounts for 31.06%. The three geographic groups seem to be mainly structured by seven alleles corresponding to 312, 323, 326, 328, 336, 342 and 362. Regarding the locus Bem23, the percentages of variance derived from the two first axes were 97% and the three geographic groups were mainly structured by seven alleles (219, 213, 223, 211, 205, 211, 202 and 216). Each SSR locus was interpreted separately allowing us to spotlight alleles that contributed the most to discriminating the three geographic groups. Among revealed alleles, several are located close to the center of the MCOA two-dimensional space and are considered as rare and group-specific alleles (Fig 2).

Based on the three geographic groups of *B. tabaci* populations, the AMOVA analysis was performed for each of the two invasive biotypes separately. When considering the B and Q biotypes, the largest part of variation was clearly observed within populations (88.241% and 89.956%, respectively) whereas variation among localities within geographic groups accounted for 0.141% and 2.063%, respectively. Results also showed that genetic differentiation among the three geographic groups was very limited (Table 4).

DAPC analysis was performed to evaluate the overall pattern of variation among natural populations and thus infer the relationships between them. DAPC results supported the existence of three genetic groups (Fig 3A). DAPC was able to assign more than 60% of all populations where they were sampled into the corresponding group. The first group (group 1) was composed of 76.86% of the North population, 17.24% of the South population and 6.89% of the Sahel population. The second group (group 2) was composed of 90% of the Sahel population and 10% of the North population. The last group was composed of 60% of the South population, 35% of the North population and 5% of the Sahel population (Fig 3B). With regard to the genotypes assigned to each group, the North group presented cross-links with both the Sahel and the South (Fig 3C).

## Species distribution modeling

An ENFA analysis was performed and used to create climatic suitability map that depicts broad-scale climatically suitable areas for *B. tabaci* (Fig 1). The ten-fold cross validation procedure indicated a good ENFA model quality according to the calculated Boyce index mean and standard deviations of 0.77 ± 0.014. These high values imply that the habitat occupied by *B. tabaci* clearly differ from the average environmental conditions found in the broader study
<table><caption>Table 2. Genetic diversity parameters calculated for tested microsatellite loci.</caption><thead><tr><th></th><th></th><th>North</th><th>Sahel</th><th>South</th><th>Overall</th></tr></thead><tbody><tr><td rowspan="4">Locus.53<sup>a</sup></td><td>Na</td><td>9.000</td><td>7.000</td><td>10.000</td><td>12.000</td></tr><tr><td>He</td><td>0.863</td><td>0.835</td><td>0.842</td><td>0.877</td></tr><tr><td>Ho</td><td>0.897</td><td>1.000</td><td>0.850</td><td>0.898</td></tr><tr><td><i>P-value</i></td><td></td><td></td><td></td><td>0.106</td></tr><tr><td rowspan="4">Locus.68<sup>b</sup></td><td>Pic</td><td></td><td></td><td></td><td>0.855</td></tr><tr><td>Na</td><td>17.000</td><td>9.000</td><td>13.000</td><td>19.000</td></tr><tr><td>He</td><td>0.914</td><td>0.860</td><td>0.905</td><td>0.924</td></tr><tr><td>Ho</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><td rowspan="4">Locus.145<sup>b</sup></td><td><i>P-value</i></td><td></td><td></td><td></td><td>0.555</td></tr><tr><td>Pic</td><td></td><td></td><td></td><td>0.910</td></tr><tr><td>Na</td><td>16.000</td><td>7.000</td><td>9.000</td><td>17.000</td></tr><tr><td>He</td><td>0.877</td><td>0.785</td><td>0.747</td><td>0.865</td></tr><tr><td rowspan="4">Locus.177<sup>b</sup></td><td>Ho</td><td>1.000</td><td>1.000</td><td>1.000</td><td>1.000</td></tr><tr><td><i>P-value</i></td><td></td><td></td><td></td><td>0.001*</td></tr><tr><td>Pic</td><td></td><td></td><td></td><td>0.845</td></tr><tr><td>Na</td><td>10.000</td><td>6.000</td><td>11.000</td><td>11.000</td></tr><tr><td rowspan="4">Bem23.A1<sup>c</sup></td><td>He</td><td>0.777</td><td>0.760</td><td>0.885</td><td>0.866</td></tr><tr><td>Ho</td><td>0.552</td><td>0.500</td><td>0.950</td><td>0.678</td></tr><tr><td><i>P-value</i></td><td></td><td></td><td></td><td>0.001*</td></tr><tr><td>Pic</td><td></td><td></td><td></td><td>0.844</td></tr><tr><td rowspan="4">BT4<sup>d</sup></td><td>Na</td><td>10.000</td><td>6.000</td><td>10.000</td><td>14.000</td></tr><tr><td>He</td><td>0.843</td><td>0.755</td><td>0.816</td><td>0.849</td></tr><tr><td>Ho</td><td>0.379</td><td>0.100</td><td>0.350</td><td>0.322</td></tr><tr><td><i>P-value</i></td><td></td><td></td><td></td><td>0.014*</td></tr><tr><td rowspan="4">BT159<sup>d</sup></td><td>Pic</td><td></td><td></td><td></td><td>0.825</td></tr><tr><td>Na</td><td>18.000</td><td>11.000</td><td>15.000</td><td>21.000</td></tr><tr><td>He</td><td>0.920</td><td>0.885</td><td>0.901</td><td>0.932</td></tr><tr><td>Ho</td><td>0.897</td><td>0.900</td><td>0.850</td><td>0.881</td></tr><tr><td rowspan="4">Mean NA</td><td><i>P-value</i></td><td></td><td></td><td></td><td>0.308</td></tr><tr><td>Pic</td><td></td><td></td><td></td><td>0.920</td></tr><tr><td>Na</td><td>14.000</td><td>10.000</td><td>16.000</td><td>19.000</td></tr><tr><td>He</td><td>0.870</td><td>0.885</td><td>0.914</td><td>0.912</td></tr><tr><td rowspan="4">Mean HE</td><td>Ho</td><td>0.793</td><td>0.800</td><td>0.750</td><td>0.700</td></tr><tr><td><i>P-value</i></td><td></td><td></td><td></td><td>0.217</td></tr><tr><td>Pic</td><td></td><td></td><td></td><td>0.897</td></tr><tr><td>Mean HO</td><td>13.429</td><td>8.000</td><td>12.000</td><td></td></tr></tbody></table>

NA: Number of alleles detected; HE: Expected heterozygosity; HO: Observed heterozygosity

*: significant P values; Non-significant P values are not indicated. PIC: Polymorphic index content.

<footnote><sup>a</sup>[26]</footnote>

<footnote><sup>b</sup>[28]</footnote>

<footnote><sup>c</sup>[33]</footnote>

<footnote><sup>d</sup>[25].</footnote>
<table><caption>Table 3. Alleles inferred from loci that were specifically revealed within each geographic group.</caption><thead><tr><th></th><th>Locus 53</th><th>Locus 68</th><th>Locus 145</th><th>Locus 177</th><th>Bem23</th><th>BT4</th><th>BT159</th></tr></thead><tbody><tr><th>North</th><td>132</td><td>185-193-209-216-233-246</td><td>167-176-179-210-223-237-243</td><td></td><td>205-271-386</td><td>334-356-362-365-368</td><td>272-291</td></tr><tr><th>Sahel</th><td></td><td></td><td></td><td></td><td>202</td><td>330</td><td>251</td></tr><tr><th>South</th><td>167-170</td><td>230</td><td></td><td></td><td>207-371-379</td><td>304</td><td>211-218-255-274</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0185724.t003

mean temperature of 10–36°C (Fig 4A). The broad-scale distribution shows that high latitude regions, due to their low temperature, can be considered as beyond the climatic limits of the pest. These include the eastern and western Ghats and Himalayan Mountains in India, the Zagros Mountains in Iran, Iraq and south-eastern Turkey and the Atlas Mountains in North

<figure><img src="image_5.png" /><figcaption>**Fig 2.** Multiple co-inertia analysis (MCOA) analysis. Single markers coordinated for the first two axes of the PCA. Corresponding plots are drawn on the same scale for the seven markers involved. The first two axes of the % PCA are shown. The three groups of populations are labeled within confidence ellipses (P = 0.95), with an envelope formed by the most discriminating alleles that are joined by lines. The shown barplot of Eigen-values indicates the relative magnitude of each axis with respect to total variance. Distribution of groups based on the common congruence values (Cos²) for the two components of the first axis.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0185724.g002
<table><caption>Table 4. Analysis of molecular variance (AMOVA) for biotype B and biotype Q populations of *B. tabaci*.</caption><thead><tr><th></th><th>Source of variation</th><th>Sum of squares</th><th>Variance components</th><th>Percentage variation</th></tr></thead><tbody><tr><td rowspan="5">Biotype B</td><td>Among geographical groups</td><td>7.577</td><td>0.018</td><td>0.576%</td></tr><tr><td>Among localities within groups</td><td>13.663</td><td>0.004</td><td>0.141%</td></tr><tr><td>Among populations within localities</td><td>81.308</td><td>0.339</td><td>11.042%</td></tr><tr><td>Within populations</td><td>84.000</td><td>2.710</td><td>88.241%</td></tr><tr><td>Total</td><td>186.548</td><td>3.071</td><td>100.000%</td></tr><tr><td rowspan="5">Biotype Q</td><td>Among geographical groups</td><td>6.542</td><td>0.084</td><td>2.647%</td></tr><tr><td>Among localities within groups</td><td>11.289</td><td>0.065</td><td>2.063%</td></tr><tr><td>Among populations within localities</td><td>73.508</td><td>0.169</td><td>5.334%</td></tr><tr><td>Within populations</td><td>80.000</td><td>2.857</td><td>89.956%</td></tr><tr><td>Total</td><td>171.339</td><td>3.176</td><td>100.000%</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0185724.t004

Africa (Fig 1A). In Tunisia, climatically suitable areas display a discontinuity in the West-Central region. This is likely due to the lower temperatures of the region, which is an extension of the Atlas Mountains (Figs 1B and 4B).

## Identification of dispersal paths

A maximum $F_{ST}$ value of 0.27 was set as thresholds to select genetically related samples used for corridors' construction, thus pairs of localities with pair-wise $F_{ST} > 0.27$ were filtered out. Accordingly, a total of 77 pair-locations from the 1711 possible pair-wise combinations were retained to calculate accumulated cost surfaces (Fig 5A). These results outlined potential *B. tabaci* dispersal corridors in Tunisia (Fig 5B). Regions located in the North displayed strong connectivity with the South and Sahel geographic groups as 65 out of 77 least-cost corridors pass through these regions. The gene flow between the northern and southern populations is therefore predicted to take place through indirect paths crossing the Sahel region and to a lesser extent through direct paths crossing the central part of the country (Fig 5B).

## Discussion

Analysis of the genetic structure of a large representative sample of *B. tabaci* populations throughout the agricultural areas of Tunisia was carried out to generate new insights into the traits and population dynamics of this pest.

The typical short-generation time of insects makes their population dynamics highly sensitive to climate variability and landscape configuration [52]. In this work, we developed a genetic landscape analysis of *B. tabaci* populations whose purpose was to predict the potential migration corridors between populations throughout Tunisia. By using SSR molecular markers and geographical data, our approach combined the analysis of genetic population diversity and structure, reflecting the amount of genetic exchanges between different *B. tabaci* populations and SDM predictions. The effect of both temperature rainfall and host/prey availability on the spatial distribution of populations was then considered to take into account the functional connectivity of habitats within a varied landscape matrix to subsequently identify potential connection corridors between groups of *B. tabaci* populations.

*B. tabaci* was previously described in Tunisia [21,53], but no exhaustive data were available regarding its features. This first investigation of genetic diversity of Tunisian populations of *B. tabaci* was performed using PCR-sequencing of mitochondrial cytochrome oxidase I (mtCOI). It was assumed that the largest part of variability could be attributed to differences between
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Discriminated analysis of principal components (DAPC) using SSR-based genotypes. (a) DAPC first and second ordination axes. (b) Neighbor-joining tree constructed on the DAPC distances representing the genetic relationship between the inferred 3 clusters as well as the distribution of geographic groups within each of them. (c) Location of the populations (points) and their probabilities for membership in each geographic groups. The colors correspond to the genetic clusters defined by the DAPC analysis.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0185724.g003

localities devoted to each geographic group rather than between groups themselves [21,53]. Beyond cytoplasmic DNA markers, the present work is also based on nuclear SSR markers which are more suitable and accurate for the study of the genetic structure of whitefly *B. tabaci* populations [54]. SSR markers are often used for the genetic analysis of insect populations and are effective at differentiating populations [55]. Indeed, SSR markers have been useful in revealing the distribution patterns and genetic structure of *B. tabaci* populations on the island of La Réunion [26], around the Mediterranean basin [56, 57, 58] and have also enabled the analysis of the variability of these markers to infer the geographical structure within the Asia-Pacific genetic groups [9].

Based on allelic frequencies of single-markers, the MCOA approach allowed the identification of those alleles that contribute most to the genetic topology of *B. tabaci* populations. The
<figure><img src="image_7.png" /><figcaption>**Fig 4.** (a) Gaussian kernel density estimates of annual mean temperature values within occurrence data. Red lines represent the result of 1000 permutations based on random samples of the initial occurrence data (70%). (b) Annual mean temperature map. (c) Development rate and (d) survival rate curves as a function of temperature (in°C) for immature and adult whiteflies estimated by [7].</figcaption></figure>

https://doi.org/10.1371/journal.pone.0185724.g004

MCOA two-dimensional space is a useful graphical tool that compares the usefulness of each marker for discriminating populations. This method provides a systemic view of the whitefly population structure focusing on the part of each allele in building it. Table 3 and the Fig 2 show that the frequency of 4 specific alleles allowed the separation of the population from the North (allele 362 of the locus BT4, allele 205 of the locus Bem23, allele 246 of locus 68 and allele 176 of locus 145). Besides, only the alleles 251 of locus BT159 and 202 of locus Bem23 allowed the separation of the Sahel population from the other two populations. Considering each SSR marker separately, the MCOA analysis showed that group-specific alleles are located close to the centre of the plot. Indeed, the Bem23 alleles (386, 379 and 371) that are specific to the Q
<figure><img src="image_8.png" /><figcaption>**Fig 5.** (a) Habitat suitability map of *Bemisia tabaci* reflecting climatic and host/prey constraints. The color gradient represent an habitat suitability index (high values indicate high habitat suitability). White lines represent the shortest path between each pair of genetically related samples. (b) Map of potential migration corridors between sampled populations. The color gradient represent a connectivity index (high values facilitate gene flow and low values indicate low populations connectivity).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0185724.g005

biotype are placed in the centre of the North and the South populations across the MCOA ordination space. In spite of this, the exhibited structure seems to be more devoted to the efficiency of common alleles to cluster populations into three distinct groups. When all the SSR markers are taken into account, their usefulness to fit in a common and similar structure is achieved even if their efficiency differs. Out of the seven SSR markers, loci 145 and 177 are good structuring markers that displayed three well-defined clusters with a minimum of five discriminating alleles each (176/233/173/171/231 and 257/266/269/262/255, respectively). Overall, MCOA approach is a valuable tool for the selection of efficient markers and the removal of those which are less informative [37].

Our results also suggested that there was a low level of geographic isolation between the geographic groups. The major part of genetic diversity was observed at the population level within geographical groups for both the B (88.24% of the total variance) and the Q (89.96% of the total variance) biotypes. Only 60% of defined genotypes were assigned to the geographic group where they were collected. Genotypes therefore displayed variability within the same geographic group and could be closely related to genotypes in distinct geographic groups.

DAPC results show that Tunisian populations are split into three well-differentiated genetic groups with low overlap between them. Patterns of population structure inferred from nuclear SSRs are in line with the biotype distributions derived from mitochondrial mt COI analysis. Biotype B occurs within the three delineated geographical groups whereas Biotype Q is excluded from the Sahel and seems to be restricted to the north and the south. Indeed, the first
genetic group obtained by the DAPC was composed of 51.72% of biotype Q and 48.28% of biotype B, the third group was composed of 65% of biotype Q and 35% of biotype B whereas the second group was composed only of populations of biotype B. Each genetic group seems to possess several specific alleles that may indicate long-distance dispersal followed by infrequent long-term isolation [59,60]. The predicted climatic suitability map of *B. tabaci* using the ENFA model showed geographical patterns of distribution similar to those obtained by [7] in the Northern Mediterranean using a temperature-dependent physiologically-based demographic model. Although in this study we developed a correlative SDM that statistically linked environmental data to species distribution records, its predictions fit well with those derived from the mechanistic model described by [7] that integrated the effect of temperature on physiological and demographic processes. By using only spatial observations, the correlative SDM we developed provided a good estimate of the physiological limits of *B. tabaci*. The kernel density estimates of annual mean temperature values within occurrence data show thermal tolerance limits between 10–37°C with an optimum at 26°C (Fig 4A). These limits fit well with the biodemographic parameters previously assessed (Fig 4C and 4D) that mechanistically linked *B. tabaci* populations to their environments [8]. Indeed, the survival and development rates of *B. tabaci* adult stage drop to zero at temperatures values colder than 11° C and warmer than 37° C.

The climatic suitability map predicted by the SDM, after being filtered by vegetation areas (reflecting host and/or prey availability), was used to create a resistance surface (Fig 5A) and then accumulated cost surfaces between each pair of genetically related samples (Fig 5B). This led to the creation of potential migration corridors between localities that take into account habitat suitability, geographic distance, and genetic relatedness. The resulting map of potential migration corridors of *B. tabaci*, when interpreted in conjunction with population structure analysis, suggests that *B. tabaci* migrates along South-North and Sahel-North potential corridors. The biotype distributions also suggested that the Sahel-North migration corridors concern exclusively the B biotype as the Q one is absent in the Sahel. Otherwise, both B and Q biotypes seem to follow South-North migration corridor. Our results also show that the genetic groups characterizing the southern populations are absent in the Sahel region even though no physical barrier hampers the migration of whiteflies from neighboring areas. The low genetic flow between the Southern and the Sahel regions could be explained by the direction of the prevailing winds in these regions. Small insects are generally thought to be passive fliers and are therefore dependent on air currents to carry them to new sites. *B. tabaci* weighs approximately 33 µg and it has been shown that its patterns of migration would be strongly influenced by the wind [61]. This is also in line with our previous work [21]. Our data suggested a possible migration pattern of *B. tabaci* from the South to the North regions of the country which was correlated with the spread of *Begomoviruses* vectored by whiteflies. Other studies [9] suggested a South to North trajectory reflecting the direction of prevailing winds during warmer seasons when *B. tabaci* is more active. Otherwise, the genetic structure identified in this study may also be explained by plant trade which can act as a homogenizing factor, erasing *B. tabaci* population genetic structure between geographic groups [28]. Indeed, the presences of northern populations in the South or the southern population in the North are additional pieces of evidence that human activities may deeply impact the population genetics of *B. tabaci* in Tunisia.

The knowledge gained during this study have several applications for pest management. First, by developing a map of the potential distribution of *B. tabaci*, we have provided the first pest risk analysis of *B. tabaci* to plant health in Tunisia. The predictions of areas suitable for establishment can be used to support pest management tactics and strategies such as incursion monitoring in areas not yet infested but having a high risk of infestation (cf. as is the case for north-western regions). However, it is important to note that predictions made in this study reflect only the potential distribution of *B. tabaci*. Thus, we propose in the future to refine the
method by also mapping the current distribution of *B. tabaci*. Accurately mapping such distribution and projecting the potential distribution of *B. tabaci* are both relevant tasks for managing this plant pests. For instance, control and eradication efforts should focus on the current distribution while containment efforts should focus on the interface between the current and potential distributions [62]. Mapping the current distribution will require a better knowledge of host species and *B. tabaci* interactions and the integration of wind effects on dispersion using for instance Lagrangian dispersion models.

The thorough understanding of whitefly population structure and migration corridors of whiteflies is also particularly important for proper management guidance to plan an effective integrated pest management strategy (IPM) in Tunisia. Today, particular attention is paid to enhancing biological control methods. This would require the use of *B. tabaci* predators in crops to decrease both pest population densities and dispersal. Indeed, the use of host plants infested by predators early in the season, in the proximity of crops which are often infested by whiteflies and taking into account preferential routes of migration would permit predators populations to be established prior to whitefly attack. To fully exploit an IPM over a long time scales, knowledge of the parasitic or predatory natural enemies of *B. tabaci* should be further investigated. Generalists 'predators have been shown to reduce pest pressure and to enhance natural enemy activities [63]. They act through strong competition for resources reducing pest pressure on crops and preventing new infestations by invasive pests [64]. Employing multiple natural enemies could ensure continued suppression of a target pest throughout its lifecycle [65,66]. Indeed, the effectiveness of mixed releasing of *B. tabaci* biological control agents corresponding to a mixture of the polyphagous predatory *Coleoptera* and two whitefly-specific parasitoids have been tested to be implemented in the design of biological pest control plans [67]. Exploring scenarios that include different natural enemy combinations with different migration patterns (South-North and Sahel-North) may greatly help biological control strategies against this damaging and invasive pest species in Tunisia.

**Supporting information**

S1 Fig. The PCR-RFLP (TaqI) of a fragment of the mitochondrial gene (mtCOI). L: Ladder 100 bp, Biomatik. Pattern 1: 634bp, 144bp corresponding to Q biotype. Pattern 2: 440 bp, 220bp, 144 bp corresponding to B biotype. (-): Negative control: mtCOI amplified fragment (778 bp) without TaqI digestion.
(TIF)

**Acknowledgments**

Special thanks to Dr. Anne Murray, Assistant Professor of American and British History and Dr. Benjamin Field, University of Aix-Marseille for their conscientious reading and meticulous corrections of the manuscript.

**Author Contributions**

**Conceptualization:** Ahmed Ben Abdelkrim, Faten Gorsane.

**Formal analysis:** Ahmed Ben Abdelkrim, Tarek Hattab, Hatem Fakhfakh.

**Funding acquisition:** Hatem Fakhfakh.

**Methodology:** Ahmed Ben Abdelkrim, Tarek Hattab, Faten Gorsane.

**Resources:** Mohamed Sadok Belkadhi.
**Software:** Ahmed Ben Abdelkrim, Tarek Hattab.

**Validation:** Hatem Fakhfakh, Faten Gorsane.

**Writing – original draft:** Ahmed Ben Abdelkrim, Tarek Hattab, Hatem Fakhfakh, Faten Gorsane.

**Writing – review & editing:** Faten Gorsane.

## References

1. De Barro PJ, Liu S-S, Boykin LM, Dinsdale AB. *Bemisia tabaci:* A Statement of Species Status. Annu Rev Entomol. 2011; 56: 1–19. https://doi.org/10.1146/annurev-ento-112408-085504 PMID: 20690829

2. Abd-Rabou S, Simmons AM. Survey of reproductive host plants of *Bemisia tabaci* (Hemiptera: Aleyrodiidae) in Egypt, including new host records. Entomol News. 2010; 121: 456–465.

3. Navas-Castillo J, Fiallo-Olivé E, Sánchez-Campos S. Emerging virus diseases transmitted by white-flies. Annu Rev Phytopathol. 2011; 49: 219–248. https://doi.org/10.1146/annurev-phyto-072910-095235 PMID: 21568700

4. Boykin LM. *Bemisia tabaci* nomenclature: lessons learned. *Bemisia tabaci* nomenclature. Pest Manag Sci. 2013; 70: 1454–1459. https://doi.org/10.1002/ps.3709 PMID: 24338873

5. Liu S, Colvin J, De Barro PJ. Species Concepts as Applied to the Whitefly *Bemisia tabaci* Systematics: How Many Species Are There? J Integr Agric. 2012; 11: 176–186. https://doi.org/10.1016/S2095-3119(12)60002-1

6. Brown JK. Phylogenetic biology of the *Bemisia tabaci* sibling species group. *Bemisia*: bionomics and management of a global pest. Springer; 2010. pp. 31–67.

7. Gilioli G, Pasquali S, Parisi S, Winter S. Modelling the potential distribution of *Bemisia tabaci* in Europe in light of the climate change scenario. Pest Manag Sci. 2014; 70: 1611–1623. https://doi.org/10.1002/ps.3734 PMID: 24458692

8. Naranjo SE, Castle SJ, De Barro PJ, Liu S-S. Population dynamics, demography, dispersal and spread of *Bemisia tabaci*. *Bemisia*: Bionomics and management of a global pest. Springer; 2009. pp. 185–226.

9. De Barro PJ. Genetic structure of the whitefly *Bemisia tabaci* in the Asia-Pacific region revealed using microsatellite markers: genetic structure of bemisia tabaci. Mol Ecol. 2005; 14: 3695–3718. https://doi.org/10.1111/j.1365-294X.2005.02700.x PMID: 16202090

10. Xu J, Lin KK, Liu SS. Performance on different host plants of an alien and an indigenous *Bemisia tabaci* from China: Performance of *Bemisia tabaci* on different plants. J Appl Entomol. 2011; 135: 771–779. https://doi.org/10.1111/j.1439-0418.2010.01581.x

11. Fresia P, Silver M, Mastrangelo T, Azeredo-Espin AMLD, Lyra ML. Applying spatial analysis of genetic and environmental data to predict connection corridors to the New World screwworm populations in South America. Acta Trop. 2014; 138: S34–S41. https://doi.org/10.1016/j.actatropica.2014.04.003 PMID: 24742908

12. Bouyer J, Dicko AH, Cecchi G, Ravel S, Guerrini L, Solano P, et al. Mapping landscape friction to locate isolated tsetse populations that are candidates for elimination. Proc Natl Acad Sci. 2015; 112: 14575–14580. https://doi.org/10.1073/pnas.1516778112 PMID: 26553973

13. Holderegger R, Häner R, Csencsics D, Angelone S, Hoebee SE. S-allele diversity suggests no mate limitation in small populations of a self-incompatible plant. Evolution. 2008; 62: 2922–2928. https://doi.org/10.1111/j.1558-5646.2008.00498.x PMID: 18752611

14. Manel S, Holderegger R. Ten years of landscape genetics. Trends Ecol Evol. 2013; 28: 614–621. https://doi.org/10.1016/j.tree.2013.05.012 PMID: 23769416

15. Hall LA, Beissinger SR. A practical toolbox for design and analysis of landscape genetics studies. Landsc Ecol. 2014; 29: 1487–1504. https://doi.org/10.1007/s10980-014-0082-3

16. Spear SF, Balkenhol N, Fortin M-J, Mcrae BH, Scribner K. Use of resistance surfaces for landscape genetic studies: considerations for parameterization and analysis: resistance surfaces in landscape genetics. Mol Ecol. 2010; 19: 3576–3591. https://doi.org/10.1111/j.1365-294X.2010.04657.x PMID: 20723064

17. Russo I-RM, Sole CL, Barbato M, von Bramann U, Bruford MW. Landscape determinants of fine-scale genetic structure of a small rodent in a heterogeneous landscape (Hluhluwe-iMfolozi Park, South Africa). Sci Rep. 2016; 6. https://doi.org/10.1038/srep29168 PMID: 27406468

18. Beheregaray LB, Caccone A. Cryptic biodiversity in a changing world. J Biol. 2007; 6: 9. https://doi.org/10.1186/jbiol60 PMID: 18177504
19. Hadjistylli M, Roderick GK, Brown JK. Global Population Structure of a Worldwide Pest and Virus Vector: Genetic Diversity and Population History of the *Bemisia tabaci* Sibling Species Group. Zhang Y, editor. PLOS ONE. 2016; 11: e0165105. https://doi.org/10.1371/journal.pone.0165105 PMID: 27855173

20. Dinsdale A, Cook L, Riginos C, Buckley YM, De Barro P. Refined Global Analysis of *Bemisia tabaci* (Hemiptera: Sternorrhyncha: Aleyrodoidea: Aleyrodidae) Mitochondrial Cytochrome Oxidase 1 to Identify Species Level Genetic Boundaries. Ann Entomol Soc Am. 2010; 103: 196–208. https://doi.org/10.1603/AN09061

21. Gorsane F, Ben Halima A, Ben Khalifa M, Bel-Kadhi M, Fakhfakh H. Molecular characterization of *Bemisia tabaci* populations in Tunisia: genetic structure and evidence for multiple acquisition of secondary symbionts. Environ Entomol. 2011; 40: 809–817. https://doi.org/10.1603/EN10162 PMID: 22251681

22. Rubinoff D. Essays: Utility of Mitochondrial DNA Barcodes in Species Conservation: DNA Barcodes and Conservation. Conserv Biol. 2006; 20: 1026–1033. https://doi.org/10.1111/j.1523-1739.2006.00372.x PMID: 16922219

23. Brower AVZ. Problems with DNA barcodes for species delimitation: “Ten species” of *Astraptes fulgerator* reassessed (Lepidoptera: Hesperiidae). Syst Biodivers. 2006; 4: 127–132. https://doi.org/10.1017/S147720000500191X

24. Ballard JWO, Whitlock MC. The incomplete natural history of mitochondria. Mol Ecol. 2004; 13: 729–744. https://doi.org/10.1046/j.1365-294X.2003.02063.x PMID: 15012752

25. Tsagkarakou A, Roditakis N. Isolation and characterization of microsatellite loci in *Bemisia tabaci* (Hemiptera: Aleyrodidae). Mol Ecol Notes. 2003; 3: 196–198.

26. Delatte H, David P, Granier M, Lett JM, Goldbach R, Peterschmitt M, et al. Microsatellites reveal extensive geographical, ecological and genetic contacts between invasive and indigenous whitefly biotypes in an insular environment. Genet Res. 2006; 87: 109. https://doi.org/10.1017/S0016672306008135 PMID: 16709274

27. Tsagkarakou A, Tsigenopoulos CS, Gorman K, Lagnel J, Bedford ID. Biotype status and genetic polymorphism of the whitefly *Bemisia tabaci* (Hemiptera: Aleyrodidae) in Greece: mitochondrial DNA and microsatellites. Bull Entomol Res. 2007; 97: 29. https://doi.org/10.1017/S000748530700466X PMID: 17298679

28. Dalmon A, Halkett F, Granier M, Delatte H, Peterschmitt M. Genetic structure of the invasive pest *Bemisia tabaci*: evidence of limited but persistent genetic differentiation in glasshouse populations. Heredity. 2008; 100: 316–325. https://doi.org/10.1038/sj.hdy.6801080 PMID: 18073781

29. Gauthier N, Dalleau-Clouet C, Bouvret M-E. Twelve new polymorphic microsatellite loci and PCR multiplexing in the whitefly, *Bemisia tabaci*. Mol Ecol Resour. 2008; 8: 1004–1007. https://doi.org/10.1111/j.1755-0998.2008.02163.x PMID: 21585955

30. Bickford D, Lohman DJ, Sodhi NS, Ng PKL, Meier R, Winker K, et al. Cryptic species as a window on diversity and conservation. Trends Ecol Evol. 2007; 22: 148–155. https://doi.org/10.1016/j.tree.2006.11.004 PMID: 17129636

31. Elith J, Leathwick JR. Species distribution models: ecological explanation and prediction across space and time. Annu Rev Ecol Evol Syst. 2009; 40: 677–697.

32. Bosco D, Loria A, Sartor C, Cenis JL. PCR-RFLP identification of *Bemisia tabaci* biotypes in the Mediterranean basin. Phytoparasitica. 2006; 34(3):243–51. ZOOREC:ZOOR14209057890.

33. De Barro PJ, Scott KD, Graham GC, Lange CL, Schutze MK. Isolation and characterization of microsatellite loci in *Bemisia tabaci*. Molecular Ecology Notes. 2003; 3(1):40–3. ISI:000181651400014.

34. Frohlich D, Torres-Jerez I, Bedford I, Markham P, Brown J. A phylogeographic analysis of the Bemisia tabaci species complex based on mitochondrial DNA markers. Molecular Ecology. 1999; 8:1593–602.

35. Mckenzie CL, Hodges G, Osborne LS, Byrne FJ, Shatters RG Jr. Distribution of *Bemisia tabaci* (Hemiptera: Aleyrodidae) Biotypes in Florida–Investigating the Q Invasion Journal of Economic Entomology 2009; 102(2): 670–676. https://doi.org/https://doi.org/10.1603/029.102.0227 PMID: 19449648

36. Mckenzie CL, Bethke JA, Byrne FJ, Chamberlin JR, Dennehy TJ, Dickey AM, Gilrein D, Hall PM, Ludwig S, Oetting RD, Osborne LS, Schmale L, Shatters RG Jr. Distribution of *Bemisia tabaci* (Hemiptera: Aleyrodidae) Biotypes in North America After the Q Invasion Journal of Economic Entomology 2012; 105(3):753–766. https://doi.org/https://doi.org/10.1603/EC11337 PMID: 22812110

37. Laloë D, Jombart T, Dufour AB, Moazami-Goudarzi K. Consensus genetic structuring and typological value of markers using multiple co-inertia analysis. Genet Sel Evol. 2007; 39: 545–567. https://doi.org/10.1186/1297-9686-39-5-545 PMID: 17897596

38. Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC genetics 2010; 11(1): 94.

39. Nei M. Genetic distance between populations. Am Nat. 1972; 106: 283–292.
40. Hirzel AH, Hausser J, Chessel D, Perrin N. Ecological-Niche Factor Analysis: How to Compute Habitat-Suitability Maps without Absence Data? Ecology. 2002; 83: 2027. https://doi.org/10.2307/3071784

41. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005; 25: 1965-1978. https://doi.org/10.1002/joc.1276

42. Thuiller W. Patterns and uncertainties of species' range shifts under climate change. Glob Change Biol. 2004; 10: 2020-2027. https://doi.org/10.1111/j.1365-2486.2004.00859.x

43. Hattab T, Albouy C, Lasram FBR, Somot S, Le Loc'h F, Leprieur F. Towards a better understanding of potential impacts of climate change on marine species distribution: a multiscale modelling approach: Threatened coastal region under global change. Glob Ecol Biogeogr. 2014; 23: 1417-1429. https://doi.org/10.1111/geb.12217

44. Boyce MS, Vernier PR, Nielsen SE, Schmiegelow FK. Evaluating resource selection functions. Ecol Model. 2002; 157: 281-300.

45. Hirzel AH, Le Lay G, Helfer V, Randin C, Guisan A. Evaluating the ability of habitat suitability models to predict species presences. Ecol Model. 2006; 199: 142-152. https://doi.org/10.1016/j.ecolmodel.2006.05.017

46. Bayhan E, Ulusoy MR, Brown JK. Host range, distribution, and natural enemies of *Bemisia tabaci* "B biotype" (Hemiptera: Aleyrodidae) in Turkey. J Pest Sci. 2006; 79: 233-240. https://doi.org/10.1007/s10340-006-0139-4

47. Stevenson-Holt CD, Watts K, Bellamy CC, Nevin OT, Ramsey AD. Defining landscape resistance values in least-cost connectivity models for the invasive grey squirrel: a comparison of approaches using expert-opinion and habitat suitability modelling. PloS One. 2014; 9: e112119. https://doi.org/10.1371/journal.pone.0112119 PMID: 25380289

48. Wright S. The Interpretation of Population Structure by F-Statistics with Special Regard to Systems of Mating. Evolution. 1965; 19: 395. https://doi.org/10.2307/2406450

49. Dijkstra EW. A note on two problems in connexion with graphs. Numer Math. 1959; 1: 269-271.

50. Etherington TR, Perry GLW, Cowan PE, Clout MN. Quantifying the Direct Transfer Costs of Common Brushtail Possum Dispersal using Least-Cost Modelling: A Combined Cost-Surface and Accumulated-Cost Dispersal Kernel Approach. Driscoll DA, editor. PLoS ONE. 2014; 9: e88293. https://doi.org/10.1371/journal.pone.0088293 PMID: 24505467

51. R Core Team. R foundation for statistical computing. Vienna Austria. 2013:3.

52. Zidon R, Tsueda H, Morin E, Morin S. Projecting pest population dynamics under global warming: the combined effect of inter-and intra-annual variations. Ecol Appl. 2016; 26: 1198-1210. PMID: 27509758

53. Laarif A, Saleh D, Clouet C, Gauthier N. Regional co-occurrence between distinct *Bemisia tabaci* species in Tunisia with new insights into the role of host plants. Phytoparasitica. 2015; 43: 135-150. https://doi.org/10.1007/s12600-014-0437-y

54. Simón B, Cenis J, De La Rúa P. Distribution patterns of the Q and B biotypes of *Bemisia tabaci* in the Mediterranean Basin based on microsatellite variation. Entomol Exp Appl. 2007; 124: 327-336.

55. Cornuet JM, Luikart G. Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics. 1996; 144: 2001-2014. PMID: 8978083

56. Saleh D, Laarif A, Clouet C, Gauthier N. Spatial and host-plant partitioning between coexisting *Bemisia tabaci* cryptic species in Tunisia. Population Ecology. 2012; 54(2):261-74. https://doi.org/10.1007/s10144-012-0303-z

57. Tsagkarakou A, Mouton L, Kristoffersen JB, Dokianakis E, Grispou M, Bourtzis K. Population genetic structure and secondary endosymbionts of Q *Bemisia tabaci* (Hemiptera: Aleyrodidae) from Greece. Bulletin of Entomological Research. 2012; 102(3):353-65. https://doi.org/10.1017/S0007485311000757 PMID: 22280837.

58. Gauthier N, Clouet C, Perrakis A, Kapantaidaki D, Peterschmitt M, Tsagkarakou A. Genetic structure of *Bemisia tabaci* Med populations from home-range countries, inferred by nuclear and cytoplasmic markers: impact on the distribution of the insecticide resistance genes. Pest Management Science. 2014; 70 (10): 1477-91. https://doi.org/10.1002/ps.3733 PMID: 24458589.

59. Lowe WH, Allendorf FW. What can genetics tell us about population connectivity?: genetic and demographic connectivity. Mol Ecol. 2010; 19: 3038-3051. https://doi.org/10.1111/j.1365-294X.2010.04688.x PMID: 20618697

60. Horne JB, Momigliano P, Welch DJ, Newman SJ, Van HERWERDEN L. Limited ecological population connectivity suggests low demands on self-recruitment in a tropical inshore marine fish (*Eleutheronema tetradactylum*: Polynemidae): *Eleutheronema tetradactylum*. Mol Ecol. 2011; 20: 2291-2306. https://doi.org/10.1111/j.1365-294X.2011.05097.x PMID: 21518062
61. Byrne DN, Rathman RJ, Orum TV, Palumbo JC. Localized migration and dispersal by the sweet potato whitefly, *Bemisia tabaci*. Oecologia. 1996; 105: 320–328. https://doi.org/10.1007/BF00328734 PMID: 28307104

62. Hattab T, Garzon Lopez CX, Ewald M, Skowronek S, Aerts R, Horen H, et al. A unified framework to model the potential and realized distributions of invasive species within the invaded range. Divers Dis-trib. 2017; https://doi.org/10.1111/ddi.12566

63. Jaworski CC, Bompard A, Genies L, Amiens-Desneux E, Desneux N. Preference and Prey Switching in a Generalist Predator Attacking Local and Invasive Alien Pests. Dickens JC, editor. PLoS ONE. 2013; 8: e82231. https://doi.org/10.1371/journal.pone.0082231 PMID: 24312646

64. Ragsdale DW, Landis DA, Brodeur J, Heimpel GE, Desneux N. Ecology and Management of the Soybean Aphid in North America. Annu Rev Entomol. 2011; 56: 375–399. https://doi.org/10.1146/annurev-ento-120709-144755 PMID: 20868277

65. Boivin G, Brodeur J. Intra-and interspecific interactions among parasitoids: mechanisms, outcomes and biological control. Trophic and guild in biological interactions control. Springer; 2006. pp. 123–144.

66. Mirande L, Desneux N, Haramboure M, Schneider MI. Intraguild predation between an exotic and a native coccinellid in Argentina: the role of prey density. J Pest Sci. 2014; 88: 155–162. https://doi.org/10.1007/s10340-014-0597-z

67. Tan X, Hu N, Zhang F, Ramirez-Romero R, Desneux N, Wang S, et al. Mixed release of two parasitoids and a polyphagous ladybird as a potential strategy to control the tobacco whitefly *Bemisia tabaci*. Sci Rep. 2016; 6. https://doi.org/10.1038/srep28245 PMID: 27312174