http://orcid.org/0000-0002-6701-6184
Figures
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.
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, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: 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 consisting 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 compatibility, 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 distribution 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 heterogeneity shape the distribution of pest populations and gene flow between landscape patches is crucial 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 resistance or its inverse landscape permeability can be modeled to (i) identify landscape and environmental features that constrain genetic connectivity, (ii) and to perform a spatially explicit exploration of genetic diversity and population structure, mainly to inform resource management 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 environmental 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 cytochrom 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 including 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 questions 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.
(A) across the geographical extent used for model calibration and (B) across Tunisia. Blue dots represent recorded occurrences.
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)31 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), 1Uof 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. 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 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 phytophagous 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 FST paire-wise genetic distance matrix, we extracted likely genetically identical samples: pairs with Wright's FST 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 micosatellite 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 29were 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).
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 (Cos2) for the two components of the first axis.
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).
(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.
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 area. This indicates that a species-specific habitat selection process takes place. The occurrences records used to calibrate the ENFA model indicate that B. tabaci is limited by an annual 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 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).
(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].
Identification of dispersal paths
A maximum FST value of 0.27 was set as thresholds to select genetically related samples used for corridors’ construction, thus pairs of localities with pair-wise FST> 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).
(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).
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 cytochrom oxidase I (mtCOI). It was assumed that the largest part of variability could be attributed to differences between 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 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 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.
https://doi.org/10.1371/journal.pone.0185724.s001
(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.
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. pmid:20690829
- 2. Abd-Rabou S, Simmons AM. Survey of reproductive host plants of Bemisia tabaci (Hemiptera: Aleyrodidae) 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 whiteflies. Annu Rev Phytopathol. 2011;49: 219–248. pmid:21568700
- 4. Boykin LM. Bemisia tabaci nomenclature: lessons learned: Bemisia tabaci nomenclature. Pest Manag Sci. 2013;70: 1454–1459. 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.
- 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. 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. 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.
- 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. 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. 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. pmid:18752611
- 14. Manel S, Holderegger R. Ten years of landscape genetics. Trends Ecol Evol. 2013;28: 614–621. pmid:23769416
- 15. Hall LA, Beissinger SR. A practical toolbox for design and analysis of landscape genetics studies. Landsc Ecol. 2014;29: 1487–1504.
- 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. 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. pmid:27406468
- 18. Beheregaray LB, Caccone A. Cryptic biodiversity in a changing world. J Biol. 2007;6: 9. 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. 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.
- 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. pmid:22251681
- 22. Rubinoff D. Essays: Utility of Mitochondrial DNA Barcodes in Species Conservation: DNA Barcodes and Conservation. Conserv Biol. 2006;20: 1026–1033. 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.
- 24. Ballard JWO, Whitlock MC. The incomplete natural history of mitochondria. Mol Ecol. 2004;13: 729–744. 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. 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. 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. 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. 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. 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. 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 inNorth America After the Q Invasion Journal of Economic Entomology 2012;105(3):753–766. 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. 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.
- 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.
- 42. Thuiller W. Patterns and uncertainties of species’ range shifts under climate change. Glob Change Biol. 2004;10: 2020–2027.
- 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.
- 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.
- 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.
- 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. pmid:25380289
- 48. Wright S. The Interpretation of Population Structure by F-Statistics with Special Regard to Systems of Mating. Evolution. 1965;19: 395.
- 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. 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.
- 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.
- 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. 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. 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. 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. 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. 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 Distrib. 2017;
- 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. 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. 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.
- 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. pmid:27312174