The sample includes nine Zaprionus species for which we had 10 or more individuals per geographical location (Table 1): Z. africanus, Z. indianus, Z. aff. proximus, Z. aff. vittiger, Z. davidi, Z. taronus, Z. tuberculatus, Z. sepsoides, Z. vittiger. They were collected in Senegal, Guinea, Cameroon and Gabon (Fig. 1) between 2005 and 2010 using banana traps or by sweeping insect nets over different substrates and over the vegetation. The taxonomic determination of three species remains to be confirmed. They were identified by barcoding (Hebert et al., 2003) and using the Internal Transcribed Spacer (ITS) nuclear sequence (Eickbush and Eickbush, 2007). They could be new species. Wild-caught samples were kept in 70% ethanol until DNA was extracted from single males, using a Qiagen DNeasy Blood and Tissue Kit (Qiagen GmbH, Hilden, Germany), following the manufacturer's instructions. Males were used in order to avoid genetic contamination by mates in fertilized females.

We PCR-amplified a 546 bp fragment of the mitochondrial DNA cytochrome c oxidase 1 gene (COI) using primers Lep F1 (5′-WTTCWACCAATCATAAAGATATTGG-3′) and Lep R1 (5′-TAAACTTCTGGRTGTCCAAAAAATCA-3′), redesigned from those used by Hebert et al. (2004). Each PCR mix (30 μL) contained 3 μL of Fermentas Dream Taq Buffer, 200 μM of each dNTPs, 0.5 μM of each primer, 0.9 u of Fermentas Dream Taq DNA polymerase and 1 μL of DNA. Amplification included an initial denaturation step at 94 °C for 5 min, 10 cycles (94 °C for 30 s, 45 °C for 1 min and 30 s, 72 °C for 1 min), followed by 35 cycles (94 °C for 30 s, 51 °C for 1 min and 30 s, 72 °C for 1 min) and a final elongation step at 72 °C for 5 min. DNA sequencing was performed by Sanger Sequencing using a 96-capillary ABI 3730×l DNA sequencer (Applied Biosystems).

Sequence quality was checked using CodonCode Aligner 2.0.6 (CodonCode Corporation, Centerville, MA, USA) and data were analyzed using MEGA 5.05 (Tamura et al., 2011). Intra-population molecular variation was characterized using DnaSP 5.10 (Librado and Rozas, 2009) by the number of segregating sites (S), the number of haplotypes (h), haplotype diversity (H_d) and nucleotide diversity (using Tajima's π and Watterson's θ_s). Two species, Z. vittiger and Z. sepsoides, were monomorphic and were removed from further analyses. Departure from neutrality was assessed using Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) and their significance was computed from 1000 coalescent simulations using Arlequin 3.5 (Excoffier and Lischer, 2010). Tajima's D is based on the difference between two coalescent estimators of the effective population size, which under neutrality is distributed around zero. Departure from neutrality skews this difference towards positive values in case of a genetic bottleneck, and towards negative values in case of a population expansion (Tajima, 1989). Fu's Fs is based on the expected number of alleles under a certain value of the effective population size, and it assumes values close to zero under stationarity (Fu, 1997). Population expansions skew this statistic towards negative values; bottlenecks skew it towards positive values.

In species represented by more than one population, we computed the pairwise Φ_st and its significance (from 1000 permutations) using Arlequin 3.5 (Excoffier and Lischer, 2010).

Parameters of the instantaneous demographic expansion model (Rogers and Harpending, 1992) were estimated by means of DnaSP 5.10 using the distribution of pairwise differences between chromosomes (the mismatch distribution). Comparing this distribution to predictions of the coalescent theory can be used to detect changes in the effective population size (Schneider and Excoffier, 1999). Expansion times were scaled in calendar units assuming a mutation rate of 5.7 × 10⁻⁸ per site per year (Tamura, 1992) and 10 generations per year. We discuss these assumptions in the discussion section below.

The variation of the effective population size over time was reconstructed by means of the Extended Bayesian Skyline Plot (EBSP, Heled and Drummond, 2008) as implemented in Beast 1.7 (Drummond et al., 2012). Beast 1.7 is a software using Monte Carlo Markov Chain (MCMC), a simulated random walk through parameter space for the purpose of sampling from the posterior probability distribution of interest (Gilks et al., 1996). Each analysis was run until the posterior parameters converged (between 10,000,000 and 50,000,000 iterations) with a 10% burn-in and a thinning of 1000. We assumed the HKY + γ model of sequence evolution.

Out of 21 Zaprionus species identified from their morphology and their molecular barcode in the studied region, nine were represented by population samples, including ten individuals or more and were retained for the molecular variation study. Results are shown in Table 1. A striking contrast appears in the level of variation across species, with three different patterns. Firstly, two species, Z. vittiger (n = 81 individuals in three populations) and Z. sepsoides (n = 82 in four populations) were nearly monomorphic despite a large sample size. The second species is known to be infected by Wolbachia, a commensal intracellular bacterium (Cordaux et al., 2008). The lack of variation in the COI mitochondrial gene probably results from hitchhiking by this bacterium, which is transmitted matrilinearly in oocytes, and is known to spread in host populations through a transmission bias conferring an advantage to its host matrilines. We do not know whether or not a similar argument applies to Z. vittiger. Second, two species, Z. africanus (n = 21 in one population) and Z. indianus (n = 99 in 6 populations) present a very high level of variation, with a value of π = 1.3% in the first species, and π = 1.2–1.8% among populations of the other species. They will be referred to below as the high-variation species. Third, the remaining five species (Z. aff. proximus, Z. aff. vittiger, Z. davidi, Z. taronus, and Z. tuberculatus) represented by a total of ten populations, all show a low level of variation, ranging from π = 0.2% to π = 0.4%, depending on the population. They will be referred to below as the low-variation species. Tajima's D was significant in all low- and high-variation species (but in only seven populations out of ten in low-variation species, and not in a single population of the high-variation species), whereas Fu's F_s was significant in all species and populations for both low- and high-variation species. In all cases, Tajima's D and Fu's F_s were negative, including non-significant cases.

Significant pairwise Φ_st values were found in two species, Z. indianus and Z. taronus. Data are shown in Table 2 for these two species. Zaprionus taronus showed high and significant values for all three comparisons. As expected, the two populations from the Cameroon volcanic line, Mt Kupe and Mt Fako, were closer to each other than either was to the Guinean population from Mt Nimba. The fact that all populations were significantly different is evidence of substantial divergence between populations from this species, which lives almost exclusively in montane forests. For instance, the collection sites at Mt Kupe and Mt Fako are only 100 km distant, but they are separated by a broad lowland area. In Z. indianus, populations were not significantly different from each other, despite an extended sampling covering distant areas. A surprising fact is that the Mt Nimba population (Gouan), differed from the Mt Kupe and Libreville populations, but did not differ from the populations that are geographically closer to them, such as Mbiyeh (for Mt Kupe) and Lope (for Libreville). This can be explained by the fact that the Mt Kupe and Libreville populations are much less variable than the other populations whereas Gouan is the most variable of all. This increases significance in the random genetic drift of the smaller populations increases genetic distance in these two pairwise comparisons. We have no explanation for the lower variation of these populations. Libreville is a suburban population, whereas Mt Kupe is a forested montane environment. Both locations are not typical of this species, which is generally found in open areas.

Since negative values of Tajima's D and Fu's F_s support a demographic expansion, we checked the plausibility of this interpretation by visual inspection of the mismatch distribution (Fig. 2). In agreement with this, it appeared that the distribution was unimodal for almost all species and also for each population of Z. taronus (data not shown). The modes of the distributions fell into two classes corresponding respectively to the low-variation and the high-variation species. The average mismatch value in the low-variation group was k = 1.57 (range 1.038–1.845). In the high-variation group, the average mismatch was k = 7.75 (6.933–8.567). Thus, the average mismatch in the latter group was about five times higher than in the former.

For each of these species, the output of the EBSP is shown in Fig. 3, while the parameters inferred by the instantaneous population expansion model are given in Table 3. Using Tamura's molecular clock, the time span since the beginning of population growth was estimated to be 138 kya for Z. indianus and 111 kya for Z. africanus. Even though the uncertainty around these estimates is very large, these values are very close to each other. This estimate corresponds roughly to the onset of the last interglacial. For the low-variation group, the estimates were also very close to each other (average: 25.2 kya, range: 17–30 kya). This period approximately matches the time elapsed since the last glacial maximum (22–26 kya).

Some limits inherent to the methods used need to be mentioned. Firstly, the generation time of drosophilids in the wild is unknown. While the time of development as a function of temperature is well established for several Drosophila species in the laboratory, generation time in the wild is difficult to estimate, since it depends on several factors, including the temperature of development which can be different in lowlands and in highlands, and on the survival of adults in natural conditions. Second, the neutral mutation rate in Drosophila mitochondrial DNA is difficult to evaluate. Several estimates of the mutation rate are available in Drosophila melanogaster. Tamura (1992) estimate of the mitochondrial DNA molecular clock (5.7 × 10⁻⁸ per site per year) is based on the assumption from geological evidence that the radiation of the Drosophila genus occurred 40 million years ago. This date could be an underestimate due to a lack of geological record, thus, the neutral mutation rate could be overestimated and lead to underestimated time inferences. Haag-Liautard et al. (2008)’s estimate of the mitochondrial DNA mutation rate (6.2 × 10⁻⁸ per site per generation) is based on laboratory mutation accumulation experiments. It includes some proportion of heteroplasmic non-synonymous changes, which could represent mildly deleterious mutations subsisting in experimental lines. Thus this is probably an overestimate of the neutral mutation rate. Obbard et al. (2012)’s estimates of the nuclear molecular clock in Drosophila are based, either on divergence data between species endemic to Hawaiian Islands of known age, or on experimentally derived mutation rates. The authors note that the rate calibrated by Hawaiian Islands is probably an overestimate since it assumes that the divergence of new species started soon after the emergence of a new island. The second estimate, calibrated by the mutation rate, is compatible with geological data on Drosophila radiation. Below, we use Tamura's estimate of the neutral clock, but it should be remembered that the most recent account on the nuclear mutation rate in Drosophila concludes that estimates of divergence dates in Drosophila should still be treated with some caution.

Attention should also be paid to cautiously interpreting mitochondrial DNA variation data. Since mitochondrial DNA is a non-recombining genetic system, any bias in its genetic transmission affects it as a whole. This includes infection by Wolbachia, an intracellular endosymbiotic bacterium. Differential genetic transmission between matrilines affects the transmission of mitochondrial DNA and results in hitchhiking. Hitchicking (Maynard-Smith and Haigh, 1974) results from the simultaneous increase in frequency of advantageous alleles and of neighboring neutral variants, since they do not segregate as independent Mendelian factors. This can result in the “selective sweep” of local variation in whole chromosomal fragments (Kaplan et al., 1989). Selective weeps can be detected by a number of selective neutrality tests (Depaulis et al., 2004), which are also used to detect demographic events (Depaulis et al., 2003; Mousset et al., 2004). Mitochondria and Wolbachia co-occur in matrilines, thus, the endosymbiotic bacterium can induce a selective sweep of variation in the mitochondrion. Interestingly, a Wolbachia strain has been isolated in Z. sepsoides (Cordaux et al., 2008), and its presence can explain the low mitochondrial DNA variation observed for this species in our data. The very low level of mitochondrial DNA variation in Z. vittiger could be explained in the same way, but this species has never been tested for Wolbachia infection. For the other species, the selective sweep hypothesis is less plausible. In this study, the evolutionary history drawn is that of mitochondrial DNA rather than that of the species. Mitochondrial DNA is a very useful genetic system for putting forward new hypotheses on population history, but genetically independent data from nuclear genes should be used to confirm these conclusions.

The most striking result observed in our study is the long-lasting population expansion of Z. indianus and Z. africanus, the expansion of the former species being associated with a large increase in effective population size, whereas that of the latter species is moderate. Zaprionus indianus is a worldwide distributed species (Yassin et al., 2008b) whereas Z. africanus was recently described as a new species (Yassin and David, 2010) on the basis of mitochondrial analysis and experimental evidence of reproductive isolation of two lines from Gabon and Uganda with the former species. The striking increase in the population size of Z. indianus over a very long period provides new evidence on the remarkable evolutionary success of this species. This is a widely studied invasive species which succeeded in extending its geographical range from Africa to India, the Middle-East, southern Europe, and the Americas over the fifty last years (Chassagnard and Kraaijeveld, 1991; David et al., 2006; Yassin et al., 2008b), revealing its ability to spread into new environments. In tropical Africa, it is more abundant in open environments than in rain forests (Lachaise, 1974; Yassin et al., 2008b), which could explain its ability to invade seasonally dry areas, such as Brazilian savannas (Loh and Bitner-Mathé, 2005) and arid areas, such as oases in the Egyptian desert (Yassin et al., 2009). Our results suggest that the expansion time of both species (about 130 kya according to Tamura (1992)’s molecular clock) is consistent with an expansion beginning in the last interglacial period. Since there is no indication of a bottleneck coinciding with the period of the last glacial maximum (22–26 kya), our findings suggest either that these two species adapted to an environment which remained sufficiently widespread to support large populations, or that this environment was fragmented into numerous patches with large carrying capacities. In line with this, Z. indianus seems to be adapted to a wide range of environments. Our knowledge of the ecology of Z. africanus is limited. The samples used in this study, for which a species identification is confirmed by a molecular barcode, originated from five populations occurring at variable elevations between lowlands and highlands, suggesting that Z. africanus is adapted to a wide spectrum of environments.

The samples of the remaining five species, Z. aff. proximus, Z. aff. vittiger, Z. davidi, Z. taronus and Z. tuberculatus, were all caught in forested environments. The fact that all of these species show evidence of low-variation and of a population expansion beginning with the last glacial maximum remarkably complements the conclusions drawn for the previous two species, Z. africanus and Z. indianus. Thus, the same rationale applies to both cases, with the appropriate changes, while species typical of open environments were expanding during the last glacial period, when savanna was the dominant ecosystem; forest species began expanding after the Last Glacial Maximum. A thorough characterization of the ecological distribution of all five species is underway.

Observations similar to those made on Zaprionus were previously made in the melanogaster subgroup of Drosophila. From nuclear DNA variation, Li and Stephan (2006) estimated the East African population of Drosophila melanogaster to have undergone a population growth having begun 60,000 years ago. Baudry et al. (2006) found evidence of population expansion in several D. simulans populations from eastern Africa, but not from central or southern Africa. Since no reduction of population size during the last glaciation was observed in D. melanogaster, Stephan and Li (2007) suggested that the wild-to-domestic habitat shift of D. melanogaster postulated by Lachaise et al. (1988) occurred before the Last Glacial Maximum thus sheltering non-forest-dwelling populations from a decrease in population size. The fact that no reduction in population size is observed in Z. indianus and Z. africanus rather suggests that the expansion event observed in D. melanogaster was general among species living in open environments.

In the forest species D. teissieri, Cobb et al. (2000) showed that African populations are genetically strongly differentiated, in relation with the history of African forests. Thus, among these three sibling species from the melanogaster subgroup, D. melanogaster and D. simulans are preferentially adapted to open environments, and have become invading species throughout the world, whereas D. teissieri is limited to forest ecosystems. It is striking that, like Zaprionus indianus, the eastern populations of the first two species probably underwent an old population expansion, whereas the last one, which is ecologically closer to Z. taronus, is likewise geographically structured between African forests. These cases could represent two main classes of recurrent conditions, reflecting two different kinds of responses of African drosophilids to recent climate change, depending on their adaptation to the main biomes, forests and savannas.