Introduction

The spread of pests outside of their native range occurs when the new environment and the introduced species are compatible (Facon et al., 2006). Hence, according to theses authors, three scenarios may account for biological invasions: (1) evolutionary change of the invader, (2) change in the migration process, possibly as a result of human activity and (3) environmental change, such as global warming or, more locally, glasshouses. Several tropical and sub-tropical pests are now established in northern countries due to glasshouses (Kiritani, 2006). The patchy distribution and the enclosed environment of glasshouses are expected to induce genetic drift and decrease genetic diversity, as was observed with glasshouse populations of the native aphid Aphis gossypii in France (Fuller et al., 1999). The genetic structure of invasive pests in glasshouses outside of their native range has rarely been investigated.

Here we report such a study with Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae). This pest, originating from warm to hot climates, has recently become established in northern countries due to glasshouses as well as increased human trade activities, particularly of ornamentals (Brown et al., 1995b). Although described for more than 100 years, it was only in the 1980s that B. tabaci became important as a worldwide pest and a vector for several viruses (Jones, 2003). Increases in plant damage were linked to a change in agricultural practices, namely the intensification of crop regimes, including increased monoculture production. In spite of high insecticide use, these practices lead to the surge in whitefly populations (Brown et al., 1995a). In the United States, such surging populations were reported to be composed of invasive individuals (Costa and Brown, 1991), which were morphologically indistinguishable from the indigenous whiteflies. They could, however, be distinguished by some biological features, including efficiency of virus transmission, induction of silvering symptoms and host range. The conclusive test came from a biochemical approach, which revealed different esterase profiles (Costa and Brown, 1991). According to biological differences, the indigenous population was called ‘biotype A’ and the invasive population, ‘biotype B’ (Brown et al., 1995b). Other populations were subsequently distinguished in different geographic regions and were referred to as distinct biotypes, although later distinctions were predominantly based on genetic profile (Perring, 2001). The evidence of limited gene flow reported in various crossing experiments (for review see De Barro et al., 2005) supported the concept of a species complex (Brown et al., 1995b). Global relationships of B. tabaci were recently revealed by comparing cytochrome oxidase 1 (CO1) sequences of 366 specimens collected worldwide using a Bayesian phylogenetic technique (Boykin et al., 2007). Twelve major well-resolved genetic groups were obtained. In particular, there was a New World group containing biotype A, a Mediterranean/Asia Minor/Africa group containing biotype B, and a Mediterranean group containing biotype Q. The groups were associated with particular geographic regions, suggesting allopatric divergence (De Barro et al., 2005). The two latter groups, referred to as group B and Q, were also detected in Europe, America and Asia, as a result of recent invasion events (Boykin et al., 2007).

A simulation based on weather data and developmental temperature requirements showed that the northern limit of the open-field development area of B. tabaci in Europe was across Spain, Italy, Croatia, Albania and Greece (Reynaud, 2000), but not France, which was north of this limit. However, B. tabaci was detected in glasshouses in all European countries where it cohabitates with Trialeurodes vaporariorum. Based on morphological identification, B. tabaci was reported for the first time in Spain in the 1940s in several crops (Gomez-Menor, 1943) and in France in 1989 in glasshouses (Giustina et al., 1989). B. tabaci populations that induced silvering symptoms were later observed in cucurbit and cucumber glasshouses in France (Villevieille and Lecoq, 1992). These symptoms were similar to those reported during the same time period for the biotype B infestation in the United States (Yokomi et al., 1990). The presence of biotype B on the French Mediterranean coast was confirmed in 1995 (Guirao et al., 1997). B. tabaci was only sporadically identified in glasshouses or open-field crops and weeds in southern France until the summer of 2001, when large populations were observed in two regions, PACA and Pyrénées orientales (Figure 1). Interestingly, silvering symptoms were not detected at that time (Chabrière, Terrentroy, Trottin and Schoen, personal communication). In 2003, during an unusually hot summer, populations increased tremendously in the same areas, and several B. tabaci-transmitted viruses were reported (Desbiez et al., 2003; Dalmon et al., 2003, 2005).

Figure 1
figure 1

Sample site locations. The three main vegetable and ornamental producing areas studied in southern France are highlighted with grey colour, and their names are indicated in bold. Sampling sites (see Table 1 for details) are symbolized by black spots.

Using microsatellite and mitochondrial markers, which have consistently proven to be powerful tools for population genetics, our aim was to investigate the genetic diversity and structure of invasive B. tabaci populations on a regional and local scale, beyond the northern boundary of its outdoor development area. In particular, we wanted to test the putative impact of the patchy distribution of glasshouses on the genetic structure of 22 B. tabaci populations collected from vegetable and ornamental glasshouses in southern France over 3 years.

Materials and methods

B. tabaci sampling

B. tabaci individuals were collected from 2003 to 2005 (Table 1) from various locations in southern France (Figure 1). They were sampled from some of the major vegetable and ornamental host species grown in commercial glass or plastic houses (subsequently referred to as glasshouses) in south of France (Table 1). Each glasshouse consisted of a unique crop except in site 3 in Aquitaine, where hibiscus and begonia were grown in the same glasshouse. The population from Nice was sampled from poinsettia in the glasshouse of a botanical garden from which B. tabaci was identified for the first time in France (Bertaux, personal communication). Each B. tabaci population consisted of individuals sampled at the same time, in the same site and on the same host species (Table 1). Each population was from a single glasshouse except for population Aq1 (Aquitaine region), which was collected from adjoining glasshouses. The number of individuals tested per population was about 30, except for 7 populations, which consisted of 10 or 15 individuals. Individuals were collected from distinct plants to avoid related individuals. Whiteflies were conserved in 100% ethanol at −20 °C before DNA extraction.

Table 1 Description of the populations of Bemisia tabaci

To check for genetic structure within glasshouses, the sampling location of each individual was recorded for 10 populations (Table 1) as follows: the glasshouse was virtually divided into 30 equal plots and one female was analysed from each plot.

To study the evolution of the genetic characteristics of B. tabaci populations over time at the same site, some sites were repeatedly sampled at the beginning or end of the crop cycle (Etang de Berre), or over several years in the same village from two distinct glasshouses (Pezilla la rivière) or within the same glasshouse (St Martin de Crau, Chateaurenard and Toulouges).

Microsatellite genotyping

A total of 531 individuals were genotyped using eight microsatellite loci (Table 2), three of which are described for the first time. They were isolated from the library enriched for AC nucleotide repeats described by Delatte et al. (2006). As B. tabaci is a haplo-diploid species, only females were analysed due to their diploid state. DNA was extracted from individual B. tabaci as previously described by Delatte et al. (2005). Two microlitres of the total DNA extract were used for amplification with a QIAGEN multiplex detection kit, according to the manufacturer's instructions, in a final volume of 10 μl. One primer for each microsatellite marker was fluorescently end labelled with Fam, Hex or Ned (Table 2). Loci 11, 145, 177, BT4 and BT159 were amplified in PTC100 thermocyclers (MJ Research, Waltham, MA, USA) as follows: 15 min at 94 °C, followed by 40 cycles of 30 s at 94 °C, 1 min 30 s at 57 °C, 1 min at 72 °C, ended by 30 min at 60 °C. Loci 53, 68 and Bem23 were amplified as above except that the annealing temperature was increased from 57 to 60 °C.

Table 2 Characteristics of the microsatellite markers used

Diluted PCR products were run on a Megabace (Amersham now GE Healthcare, Orsay, France) sequencer. A 100–400 bp sizing standard was run along with samples to determine allele sizes. Microsatellite alleles were scored using the software Genetic profiler (Amersham). Individuals who failed to amplify for more than half of the loci were excluded from the genetic data analysis.

CO1 sequencing

The mitochondrial CO1 gene was partially sequenced for several individuals analysed with the microsatellite markers. For some individuals, an 817 bp fragment was amplified and sequenced following cloning as described by Tahiri et al. (2006). For others, a 755 bp sequence was determined directly from a PCR product amplified with primers C1 (Frohlich et al., 1999) and 801c, 5′-CACAMCTCTTTAAAACTRTGA-3′; amplification conditions were as described by Tahiri et al. (2006) but with 35 cycles at 95 °C for 60 s, 54 °C for 45 s and 72 °C for 30 s. Sequences were aligned with Genbank CO1 sequences that showed the highest nucleotide identity, using the Optimal Alignment method of DNAMAN (version 5.0; Lynnon BioSoft, Quebec, Canada). A phylogenetic tree was set up with a Jukes and Cantor distance matrix using the neighbour-joining method (NJ; Saitou and Nei, 1987) as implemented in the DNAMAN software.

Genetic data analysis

Intra-population genetic diversity

Micro-Checker (Brookfield, 1996; Van Oosterhout et al., 2004) was used to check microsatellite data for scoring errors and null alleles. Putative linkage disequilibrium was assessed between all pairs of loci using Genepop (Raymond and Rousset, 1995; Genepop on the Web available at http://genepop.curtin.edu.au/) following the Markov chain method with the default settings. P-values were combined across loci using Fisher's method.

Mean number of alleles per locus (Ā), observed heterozygosities (HO) and unbiased expected heterozygosities (HE, calculated following Nei, 1973) were calculated using Genetix (Belkhir et al., 2001). To test for Hardy–Weinberg (HW) equilibrium within populations, Weir and Cockerham's (1984) estimators of FIS were calculated using FSTAT 2.9.3 (Goudet et al., 1996) for each locus and overall loci. Significant difference of FIS from 0 (departure from HW equilibrium) was tested using a randomization procedure (10 000 permutations).

Individual-based inferences on population structure

To assess the level of population structure, an individual-based clustering analysis was implemented using the software Structure 2 (Pritchard et al., 2000). Assuming panmixia and the population being at HW equilibrium for each locus, this software assigned individuals to K clusters (to be set), following an iterative approach. We used the admixture model where the individual's genome may be partitioned to match the frequency distribution of alleles among the K clusters. Iteration parameters were set to a burning-in period of 60 000 iterations followed by 600 000 iterations. Five independent simulations were performed to test for the consistency of the results. The likelihood of the assignments (ln p(D)) was evaluated for K varying from 1 to 10. To test for the reliability of the K clusters partition, the ΔK estimator described by Evanno et al. (2005) was computed after the simulation results.

Population differentiation tests

Global and pairwise FST estimates were calculated using FSTAT 2.9.3 according to Weir and Cockerham (1984). The significance of population differentiation was assessed with permutation tests (10 000 permutations). In addition, a hierarchical analysis of molecular variance (AMOVA) was performed, using ARLEQUIN 3.0 (Excoffier et al., 2005), to test for geographical structure effect. This test was performed on all individuals sampled in 2003 (Table 1) except the Nice population (see the ‘Results’ section) using three hierarchical levels: individuals within populations, among populations within regions and between geographical regions. Associated P-values were calculated with 10 000 permutations.

The test for a geographical structure was complemented with an isolation by distance analysis, using the same collection restricted to the 2003 samplings. The relationship between genetic and geographic distance was assessed using a regression of FST/(1–FST) values against the logarithm of geographic distance between populations, performed with the Isolde program implemented in Genepop (Raymond and Rousset, 1995). The significance of the correlation between matrices was assessed using a Mantel test with 1000 permutations.

The isolation by distance analysis was also performed at the smallest spatial scale, within glasshouses, using the positions recorded for each of the 30 individuals (10 glasshouses). For this test, we used the isolation by distance between individuals algorithm of the Isolde program (Rousset 2000).

Mutation-drift equilibrium analysis

To check populations for a recent severe reduction in population size, we used the software Bottleneck (Piry et al., 1999) that calculates the expected equilibrium gene diversity (Heq) from the observed numbers of alleles, under the assumption of a constant population size (mutation-drift equilibrium). We chose the two-phased model recommended for microsatellite loci (parameters for TPM: variance=30.00, probability=70.00%, 1000 simulations). A possible significant heterozygosity excess (signature of bottleneck) was detected using a Wilcoxon signed-rank test, which performs well even for a limited number of loci (Luikart et al., 1998).

Results

Microsatellite variability

The number of alleles over all individuals ranged from 5 to 23 with an average (Ā) of nine alleles per locus (Table 2). Heterozygote deficiency was especially high for locus 68 (Table 2). For this locus, Micro-Checker analysis revealed the presence of null alleles, and a possible stuttering (slight changes in the allele sizes during PCR). This marker was excluded from further analyses. All other markers showed no evidence for null alleles and satisfied HW equilibrium.

Among the 531 individuals analysed, 11 amplified for fewer than four loci and were removed from subsequent analysis. With the remaining seven microsatellite markers, 99.8% of individuals showed a unique genotype, indicating the high discriminatory power of the microsatellite markers used. No significant linkage disequilibrium was detected between loci; these markers were, therefore, considered to be independent.

Intra-population genetic diversity

The HO values ranged from 0.343±0.270 to 0.527±0.218; expected values ranged from 0.354±0.253 to 0.553±0.146 (Table 3). For all populations, we observed HOHE, which lead to positive values of FIS, except for populations Aq2 and Hy03. Nevertheless, FIS values were not significantly different from 0 (P>0.001), except for population Pal03 (Table 3). With an increased threshold (P>0.05), the FIS values of populations EB03sp, SMC03, SMC05, Tou03 and Nic04 were also significant. This slight heterozygote deficiency could not be attributed to particular loci, as HOHE for all loci (Table 2).

Table 3 Genetic characteristics of the populations

The mean number of alleles (Ā) within populations was higher than 3.3 except for CH05, which showed a specific pattern (see below), and for populations with the smallest sampling size (Nic, Aq5, SLV04, N=10).

Global population structure

The Bayesian clustering analysis with STRUCTURE was performed for K varying from 1 to 10 clusters. The likelihood of assignment increased from ln p(D)=−7604 for K=1 to ln p(D)=−6856 for K=5, and then decreased until K=10. The highest likelihood was, thus, obtained at K=5, but when considering the differences in ln p(D) between successive runs, a maximum at K=4 was revealed. Adding a cluster (from K=4 to K=5) did not provide much more information on the genetic structure, and the best modal value according to ΔK index was found for K=4. Assignment results, therefore, are only presented for K=4 (Figure 2). Cluster 1 consisted of all individuals sampled in Nice, which were all assigned to this cluster with very high confidence (assignment coefficient >0.9). This partition was confirmed by a high FST value (0.522) between the Nice population and all others. On the other hand, assignment of the other individuals to clusters 2, 3 and 4 was much less pronounced. In addition, the assignment profiles largely differed between individuals of the same population, except for those of the Chateaurenard location where 69% of the individuals showed an assignment coefficient 0.9 to cluster 4. The distinction of the Chateaurenard location was also supported by a significant FST value between the Ch03, -04 and -05 populations and the others (0.067, P<0.001). Except for the Nice population, which was sampled from a glasshouse in a botanical garden, the Bayesian analysis only tentatively supported that the remaining locations (that is, glasshouses) were valid units to delineate populations. Chateaurenard populations were sampled from tomato host plants. The other six populations in group Q, which were sampled from the same host plant species, did not exhibit similarity with Chateaurenard populations.

Figure 2
figure 2

Bayesian clustering results from STRUCTURE for all samples. Each population is separated by a white vertical line. The names of the populations are indicated at the bottom of the figure (see Table 1 for details); S, SLV04 population; N, Nic04 population.

Phylogenetic analysis

A 744 bp fragment was sequenced within the CO1 gene of 26 individuals representing the different geographic regions (accession numbers AM691050-AM691059, AM691062, AM691063, AM691068-AM691077, AM691079-AM691082). A total of 11 of these 26 individuals were chosen because of their high-confidence assignment (membership probability >0.9) to one of the four clusters derived from the Bayesian analysis (Figure 2): for cluster 1, two individuals of population Nic-04; for cluster 2, one individual from EB-03sp, one from Pez-05 and two from Aq1; for cluster 3, two individuals from Pal-03 and one from Pez-05; for cluster 4, one individual from Pez-03 and one from Ch-05. The Genbank sequences that showed the highest nucleotide identity with the 26 individuals analysed were from the two invasive genetic groups, B and Q. The phylogenetic relationships between all these individuals were figured by the NJ tree performed from pairwise Jukes and Cantor distances (Figure 3). The CO1 sequence of cluster 1 individuals (Nic-04) was 100% identical to a reference individual of the B group from Antibes (31 km from Nice) collected in 1999 (Delatte et al., 2005), indicating that cluster 1 corresponded to the B group. All of the other individuals sampled in 2004 clustered with group Q reference individuals, indicating they belonged to the same genetic group, regardless of the level of admixture they displayed to clusters 2–4. The CO1 sequences of the French individuals belonging to the genetic group Q were at least 98.8% identical within the 744 bp region. They were grouped with western Mediterranean individuals and not with an Eastern group, represented here by an individual from Turkey. The sequence of nine individuals from Gat-04 and a sequence of a Moroccan individual collected in 2005 were 100% identical and formed a distinct group supported with a 97% bootstrap. Except for the individual from EB-03sp, which was slightly different, all other Q individuals, including one individual from Gat-04, could not be differentiated based on either their geographic origin or their assignment to any of the Bayesian clusters. On the contrary, it appeared that individuals from different geographic regions in southern France exhibited 100% nucleotide identity not only among themselves but also with individuals from abroad, collected in the same years. For example, there was 100% identity between individuals from Ch-05, Aq4, SLV-04, Japan (2004) and United States (2005). However, none of the group Q individuals collected previously from France (1999) or from contiguous Spain (1992–1999) were identical to the French individuals collected in this study. Unfortunately, sequences of Q individuals from Italy were not available in Genbank.

Figure 3
figure 3

Rooted NJ tree showing the genetic distance among 744 bp cytochrome oxydase 1 fragments of Bemisia tabaci individuals, either sequenced in this study (in bold) or selected from Genbank. The geographic origin of the individuals is followed by the year of sampling (when available) and the accession number of the corresponding sequence. The star (*) means that several individuals are 100% identical within the sequenced CO1 region. Thus, Gat04-01* is a representative individual of eight more individuals of the Gat-04 population; Pal03-02* is representative of individuals Aq1-02, Aq1-20, Aq5-07, and Pez05-16; Ch-05-27* is representative of Aq4-09 and SLV04-09, and Nic04-02* is representative of Nic04-08. When individuals were assigned to one of the four Bayesian clusters with a high membership probability (>90%, Figure 2), the cluster is indicated between brackets. The Jukes and Cantor distance between sequences is indicated with a scale representing a difference of one nucleotide. Numbers associated with nodes represent the percentage of 1000 bootstrap iterations supporting the nodes; only percentages 90% are indicated.

Geographic structuring

Since the population from Nice belonged to a distinct genetic group (see above), we only considered the non-Nice populations for this analysis. A hierarchical AMOVA showed that the variance within populations explained 96.7% of the total variance, 2.7% was expressed between locations within a region and less than 0.6% between geographical regions. Only the individual and location components were significant. Looking for an isolation of the populations by distance did not show evidence of geographic structuring either between locations or within glasshouses.

Evolution of genetic diversity over time

The Bayesian clustering analysis suggested that the populations of Chateaurenard sampled over 3 years (Ch03, -04, -05) evolved differently than the others. The mean number of alleles (Ā) and the expected heterozygosity (HE) both decreased over time from 3.9 to 3.0 and from 0.48 to 0.41, respectively. Significant FST values were observed between the population sampled in 2005 (Ch05) and the populations sampled in 2003 and 2004 (Ch03 and -04, respectively, 0.111 and 0.100, P<0.001). The so-called ‘evolved’ Ch05 population exhibited relatively high FST values with non-Chateaurenard populations, 0.171 on average (P<0.001). Although all populations seemed to satisfy mutation-drift equilibrium (data not shown), the likeliest cases of heterozygosity excesses (that could reveal a recent decrease in population size) were obtained for the Chateaurenard location (Ch04, P=0.22 and Ch05, P=0.28).

No similar decrease in genetic diversity was detected in the other populations sampled over the years from Toulouges, St Martin de Crau, and Pezilla la rivière.

In the same year, no genetic differences were observed at the Etang de Berre location between a sampling done at the beginning (EB03sp) and at the end (EB03au) of a tomato crop, as the FST value (0.002) was not significant (P>0.05). These results were consistent with the genetic stability of whitefly populations over years in most sampled sites.

Discussion

Unbalanced frequency of the two invasive genetic groups B and Q in southern France

We performed a genetic analysis with two kinds of molecular markers, which revealed the presence of two genetic groups in the French populations of B. tabaci. First, Bayesian clustering analysis performed with STRUCTURE software clearly divided B. tabaci populations into two entities (cluster 1 vs 2, 3 and 4), differentiated by a very high FST value (FST=0.52). Second, individuals from the two entities were assigned to genetic groups B or Q according to the match between their CO1 sequences and Genbank CO1 sequences of individuals from known genetic groups taken as references. The high FST value we found between genetic groups B and Q was consistent with the FST value documented in Moya et al. (2001), FST=0.34) using random amplified polymorphic DNA markers, and between genetic group B and the Indian Ocean group Ms (FST=0.32) (Delatte et al., 2006) using microsatellite markers.

Genetic group B was previously described in France (Guirao et al., 1997), but no exhaustive data were available regarding its distribution. Here, we provide evidence for a very limited distribution of this genetic group in southern France. In spite of a large geographical sampling from several host plant species, group B was only detected in a single population, which was shown to consist of only group B individuals. This peculiar population was from the only botanical garden glasshouse in the sampling. The group B population was further restricted to two Poinsettia plants 10 m apart. Consistent with this small population size, which can result in numerous sibling matings, the group B population showed the highest FIS value (0.255) of all populations studied. Individuals belonging to group B were not detected elsewhere in southern France. The genetic groups B and Q were also detected in other Mediterranean countries but with various ratios and sometimes overlapping distribution areas as shown in Italy (Cavalieri et al., 2004; Parella et al., 2004), southern Spain (Pascual and Callejas, 2004), Canary Islands (Gobbi et al., 2003) and Morocco (Tahiri et al., 2006). The lack of any specific and distinct geographical pattern is consistent with the invasiveness of these groups and stands in contrast with the large geographic structuring detected with Asia-Pacific populations in their native region (De Barro, 2005).

The apparently unpredictable distribution of these invasive groups may partly be related to distinct biological properties. In experimental studies without insecticide pressure, group B was more competitive than group Q (Pascual and Callejas, 2004). However, the higher resistance to insecticides of group Q (Horowitz et al., 2005, Prabhaker et al., 2005) was inferred to have lead to the outdoor displacement of group B in Spain (Pascual, 2006). Similarly, it is possible that the high insecticide pressure within French glasshouses may have selected for group Q populations rather than group B. In addition, the developmental studies done on French group Q populations showed a greater tolerance to both hot and cold conditions than reported for group B (Bonato et al., 2007). Therefore, both extreme temperature tolerance and insecticide resistance might have favoured (alone or in combination) the development of group Q individuals to the detriment of group B individuals and could reasonably account for the predominance of group Q we observed in glasshouses in southern France. In particular, outdoor conditions and non-heated glasshouses may have favoured group Q while putting group B individuals at a disadvantage.

High genetic mixing within group Q

The large predominance of group Q gave us the opportunity to study the genetic structure among populations of the same genetic group, in particular the putative impact of glasshouses on this structure in a region where year-round outdoor development was limited by cold winters. This was only possible because we used microsatellites that proved to give high resolution since nearly all individuals were discriminated. Contrary to our expectations, the patchy distribution of glasshouses did not result in high-population differentiation between glasshouses, even though some population samples were about 600 km apart. FST values between populations were low and the major part of genetic diversity was observed at the individual level, within glasshouses (96.7% of the total variation). No isolation by distance was observed suggesting that gene flow did not erode with geographical distance. These findings are consistent with a large mixing of B. tabaci individuals between French glasshouses, whatever the spatial scale considered. Intense gene flow between glasshouses suggested that they are leaky. Indeed, the migration of individuals between glasshouses could be favoured by the movement of plants between glasshouses and other human activities. This might account for the rapid spread of this pest after its recent introduction. In Crete, on the other hand, Tsagkarakou et al. (2007) showed evidence for higher genetic differentiation within group Q. FST values for six populations from Crete ranged from 0.025 to 0.437, whereas we observed a maximum of 0.171 between Ch05 and other group Q populations in France. This might reflect a longer evolutionary history of B. tabaci in Greece where it was originally described by Gennadius in 1889; it was only 100 years later that it was reported in France (Giustina et al., 1989), probably following its introduction and maintenance in glasshouses.

Our study also provides evidence for intense mixing within glasshouses. First, using a continuous model, no isolation by distance was detected within glasshouses. Second, most populations satisfied HW equilibrium. Human activities in glasshouses (for example pesticide spraying, pruning, harvest and so on) may favour whitefly flight and serve as passive transport, hence increasing within-glasshouse mixing. Conspicuously high levels of FIS were reported in previous B. tabaci population genetic studies (De Barro, 2005; Delatte et al., 2006), which contrast with our results at first glance. De Barro (2005) argued that these high heterozygote deficiencies could at least partially result from high insecticide pressures. Although intense mixing in large and stable populations should maintain the populations at HW equilibrium, our population genetic analysis may support De Barro's hypothesis since the populations showing the highest, significant (P<0.005 level) values of FIS were from two chemically controlled 2003 crops in Pyrénées Orientales (Pal03, Tou03). This region experienced large B. tabaci populations during the unusually hot summer of 2003 with a generally high B. tabaci-transmitted virus incidence in tomato plants (Desbiez et al., 2003; Dalmon et al., 2005). Interestingly the third population sampled from this region in 2003 exhibited the lowest, insignificant (P>0.05 level) value of FIS and was from a crop in which B. tabaci was biologically controlled (Pez03). Apart from some individuals of the Gat population, the CO1 sequence data did not differentiate the group Q individuals in France, which were all clustered within the same group. The 100% identity of the CO1 sequence between some of these French individuals and the cognate sequence of American, Japanese and other Mediterranean individuals was indicative of worldwide gene flow. This is consistent with the recent reports of invasions of the Mediterranean group Q into Asia and America (Zhang et al., 2005; Dennehy et al., 2006). Although some inter-glasshouse spread can possibly be explained by natural dispersal of B. tabaci (Byrne, 1999), the international and intercontinental spread likely results from anthropogenic movement. Interestingly, although the population of Gat was not differentiated by microsatellite markers, we found that 9 of the 10 individuals whose CO1 genes were sequenced formed a distinct cluster within the French group Q individuals. This discrepancy might result from the different tempo and evolutionary mode of microsatellite markers compared to the mitochondrial CO1 marker. Thus, the differences within the CO1 sequences might reveal an ancient evolutionary history that has been scrambled by microsatellite allele reshuffling fostered by intense current gene flow. It is noteworthy that 1 of the 10 sequenced individuals belonged to the main French CO1 group, suggesting that these two groups were already mixed. The previously reported differentiation between western and eastern Mediterranean populations of group Q (Tsagkarakou et al., 2007) is shown in Figure 3 by the differentiation of the Turkish CO1 sequence from the other group Q sequences from the western regions. Here, we tentatively provide evidence for additional CO1 differentiation within the western Mediterranean group Q population, with significant differentiation of a group represented by individuals from the Gat population and from Morocco. The evolutionary mechanisms that putatively lead to differentiation among group Q sequences have to be elucidated, with particular regard to whether this differentiation would have resulted from former geographic isolation among western group Q populations.

Although the differentiation between glasshouses was generally low, microsatellite analysis revealed some differentiation within group Q in Chateaurenard. The populations collected there from 2003 to 2005 not only differentiated from other group Q populations (FST=0.171 on average between the most differentiated 2005 population and all other non-Chateaurenard populations) but also over time (significant FST values in that location between the population sampled in 2005 and the populations sampled in 2003 and 2004 of 0.111 and 0.100, respectively). Although the Bottleneck test was not significant, a decrease in Ā and HE was detected. Consistently, reductions in population sizes were observed at this site from 2003 to 2004, and particularly from 2004 to 2005 (Chabrière and Trottin, personal communication), probably as a consequence of crop-free periods between the yearly crops. The reduction in population sizes could have gone undetected by the test either because of a lack of statistical power or because of some gene flow with populations outside the glasshouse (Luikart et al., 1998; Leblois et al., 2006). This limited differentiation process may be due to local genetic drift and/or selection, for example, for resistance to insecticide treatments.

Concluding remarks

The invading potential of B. tabaci is very high, as shown by its worldwide distribution (Frohlich et al., 1999; De Barro, 2006; Boykin et al., 2007). The genetic structure we observed beyond the northern limit of outdoor development of this species in this study suggested that its dispersal could not result solely from natural migration in a favourable environment, as we did not detect isolation by distance patterns. Instead, we detected large gene flow between sites, which is inconsistent with the local genetic differentiation we would expect if genetic drift was exerted on isolated glasshouses. This limited genetic differentiation, therefore, is likely the result of stochastic long-distance movement related to intense human trade activity. This is further supported by 100% nucleotide identity of CO1 sequences between B. tabaci individuals on different continents. Ornamental trades are known to have contributed to group B spreading (Brown et al., 1995b). Likewise, plantlet trade can act as a homogenizing factor, erasing B. tabaci population genetic structure among group Q individuals. The only genetic differentiation we observed in a unique location was likely due to efficient cultural practices aimed at controlling the pest (for example a significant crop-free period, integrated pest management strategy). The predominance of group Q and the putative displacement of group B, which is more sensitive to pesticide spraying, are additional pieces of evidence that human activities deeply impact the population genetics of B. tabaci in southern France. Furthermore, the putative higher tolerance of group Q populations to extreme temperature conditions (Bonato et al., 2007) can foster its further spread northward in Europe, due to better outdoor survival. These data are of primary importance for prophylactic policies against the spread of this invasive pest throughout the world.