Evaluating Gene Flow Using Selected Markers: A Case Study

Corresponding author: Thomas Lenormand, Laboratoire Génétique et Environnement, Institut des Sciences de l’Evolution (UMR 5554), Université Montpellier II, CC065 Place E. Bataillon, F-34095 Montpellier Cedex 5, France, lenorman{at}isem.univ-montp2.fr (E-mail).

Communicating editor: M. SLATKIN

	ABSTRACT

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

The extent to which an organism is locally adapted in an environmental pocket depends on the selection intensities inside and outside the pocket, on migration, and on the size of the pocket. When two or more loci are involved in this local adaptation, measuring their frequency gradients and their linkage disequilbria allows one to disentangle the forces—migration and selection—acting on the system. We apply this method to the case of a local adaptation to organophosphate insecticides in the mosquito Culex pipiens pipiens in southern France. The study of two different resistance loci allowed us to estimate with support limits gene flow as well as selection pressure on insecticide resistance and the fitness costs associated with each locus. These estimates permit us to pinpoint the conditions for the maintenance of this pocket of adaptation as well as the effect of the interaction between the two resistance loci.

ALTHOUGH evolutionary theory attempts mainly to explain past changes, its predictions can be tested by examining actual evolutionary processes in natural populations. To do so we must quantify the deterministic processes causing genetic evolution, namely, selection and gene flow, and take into account the unpredictable changes due to stochastic processes such as random drift and mutation. Among these factors only gene flow (and stabilizing selection) will oppose genetic differentiation between populations. Its evaluation is therefore required for the understanding of the evolution of populations in their "adaptive landscape" (for review, see SLATKIN 1987 ). However, methods are lacking to evaluate gene flow independently of selection, mutation, or drift.

In most cases, it is possible to determine the relative magnitude of gene flow vs. drift, and thus estimate the degree of isolation of populations. This determination enables the evaluation of the effects of different kinds of selection, of the geographic scale of a local adaptation (NAGYLAKI 1975 ), or whether such populations might be able to cross an "adaptive valley" (LANDE 1985 ). These methods, in principle, are valid for neutral and independent genes at equilibrium and for a given rate and mode of mutation. When averaged over many loci, these methods may be robust to slight departures from the assumptions (SLATKIN and BARTON 1989 ). However, they provide an estimate of the number of "effective" migrants but not of migration variance (²). Estimating this finite variance requires knowledge of population sizes and of the patterns of isolation by distance (ROUSSET 1997 ).

Another approach is to analyze directly the relative magnitude of gene flow vs. selection through clinal patterns that have been extensively studied theoretically for various selection models (FELSENSTEIN 1976 ; ENDLER 1977 ). In an infinite environment and at equilibrium, the cline slope is a robust estimate of the relative magnitude of selection vs. migration (BARTON and GALE 1993 ). Different methods can be used to infer the absolute value of each term, and in all cases, they require extra information about the system such as (1) a direct measure of dispersal, thus giving an indirect estimate of selection (ENDLER 1977 ; BARTON and HEWITT 1985 ); (2) a genotypic parameter such as heterozygote deficiency at one locus or linkage disequilibria when several loci are involved, the latter being more reliable (MALLET and BARTON 1989 ); (3) the variation of the cline through time, for instance, the speed of a wave of advance of an advantageous gene (FISHER 1937 ) or the rate of modification of the cline shape when selection or migration is not constant. Despite their potential for the understanding of the evolution of populations over their "adaptive landscape," these methods have mainly been used for the analysis of tension zones (BARTON 1982 ; SZYMURA and BARTON 1986 ; MALLET et al. 1990 ; SITES et al. 1995 ). The aim of this article is to show that these methods can also be useful for understanding the dynamics of local adaptation.

We have investigated the case of local adaptation of the mosquito Culex pipiens pipiens to organophosphate insecticides in the Montpellier area in France. This adaptation is conferred by resistance alleles at two major loci. Insecticide selection varies geographically, creating a pocket of adaptation. We have analyzed clinal patterns at these two loci to estimate selection intensities and gene flow. These estimates were used to evaluate the role of interaction between the two loci for the maintainance of the pocket of adaptation. Finally, we compared our estimates to direct or indirect estimates of selection intensity and gene flow in other studies.

	MATERIALS AND METHODS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Culex pipiens and its environment:
Larval development of the mosquito C. p. pipiens takes place mainly in anthropic pools where insecticide control occurs. Females are presumably fertilized at emergence (WEIDHASS et al. 1973) and then search for their blood meal and a site to lay ~150–200 eggs. Insecticides are applied during the breeding season, that is, approximately from April to October near Montpellier in southern France (see CHEVILLON et al. 1995 for details) and are restricted to a 20-km coastal belt. Between 1968 and 1990–91, organophosphate (OP) insecticides were exclusively used for mosquito control and have been replaced since then by the Bacillus sphaericus toxin. However, in the Montpellier area, OP insecticides are still used at large doses to control other Culicids and even to control C. p. pipiens in some situations (ANONYMOUS 1990 –1995). Additionally, residual doses of OP may be of the order of lethal concentrations for susceptible mosquitoes in many places of the treated area, possibly due to other pest controls (R. ERITJA, personal communication).

Genetics of resistance:
Two main loci are responsible for OP resistance in C. p. pipiens. The first locus, Ace.1, codes for an acetylcholinesterase (AChE1), the OP target (BOURGUET et al. 1996A ; MALCOLM et al. 1998 ) and has three alleles: Ace.1^R that codes for an insensitive AChE1; Ace.1^S that codes for a sensitive AChE1; and Ace.1^RS that corresponds to a duplication of Ace.1 and codes for both enzymes (RAYMOND et al. 1986 ; BOURGUET et al. 1996B , BOURGUET et al. 1996C ). The second locus corresponds to a "super locus," that is, two closely linked loci, Est-3 and Est-2, coding for esterases A and B, respectively (DE STORDEUR 1976 ; PASTEUR et al. 1981A , PASTEUR et al. 1981B ). Only 2 to 6 kb of DNA separate Est-3 from Est-2 (ROOKER et al. 1996 ; GUILLEMAUD et al. 1997 ). Resistance alleles at Est-3 and Est-2 induce an overproduction of esterase, resulting in gene amplification or gene regulation (ROOKER et al. 1996 ). Due to their close proximity, esterase genes are often coamplified as a single unit, which explains the complete association of resistance alleles at both loci (ROOKER et al. 1996 ; GUILLEMAUD et al. 1997 ). This coamplification justifies considering them as a single "super locus," which will hereafter be designated as Ester. In southern France, three resistance alleles have been identified at this locus. Ester¹ corresponds to an increased expression of the esterase A1, whereas Ester² and Ester⁴ correspond to coamplification of esterases A and B genes (A2-B2 and A4-B4, respectively). The recombination rate between Ester and Ace.1 has been estimated to be 14.5% (our unpublished results). The nomenclature used in this article is indicated in Table 1.

Fitness of resistant mosquitoes:
The different resistance alleles contribute unequally to OP resistance: in southern France, insensitive acetylcholinesterase alleles confer in general a resistance higher than overproduced esterases (RAYMOND et al. 1986 ; POIRIE et al. 1992 ; SEVERINI et al. 1993 ; RAYMOND and MARQUINE 1994 ; RIVET et al. 1994 ). When overproduced esterases and insensitive acetylcholinesterase are present together in the same mosquito, the insecticide resistance combines additively (RAYMOND et al. 1989 ). At both loci, resistance alleles can be assumed to be codominant for resistance (RAYMOND et al. 1987 ; POIRIE 1991 ), although the dominance of the resistance conferred by the Ace.1^R allele is environment dependent (BOURGUET et al. 1996D ). In natural populations, resistance alleles at both loci are associated with fitness costs, in the absence of insecticides, through decreases in fecundity and adult survival and an increase in larval developmental time. These fitness costs tend to be higher for insensitive acetylcholinesterase than for overproduced esterases (CHEVILLON et al. 1997 ).

Data collection:
Pupae were sampled on July 5, 1995 in 10 breeding sites along a 50-km north-south transect (Figure 1) across the treated and untreated areas studied by GUILLEMAUD et al. 1998 . They were reared until emergence, and adults were stored at -80°.

View larger version (24K): In this window In a new window Download PPT slide   Figure 1. Samples location. The dashed line represents the limit between the treated and untreated area.

For each mosquito, resistance alleles at the Ester and Ace.1 loci were determined as follows. The thorax and the abdomen were used to detect overproduced esterases using starch-gel electrophoresis (Tris-Maleate-EDTA 7.4 buffer; PASTEUR et al. 1988 ). The head was used to characterize AChE1 using the Témoin-Propoxur-Propoxur (TPP) test described by BOURGUET et al. 1996E . Overproduced esterases are dominant markers under our electrophoretic conditions, and the TPP test determines individuals displaying sensitive, resistant, or both types of acetylcholinesterase. Thus, these methods do not allow complete genotype identification. Table 1 indicates the correspondence between each genotype and its simplified code (parentheses) as well as the corresponding identified phenotype [brackets]. In addition, a phenotypic class for individuals that carry at least one resistance allele has been defined at each locus {braces}. To identify individuals at both loci, Ace.1 is indicated first, followed by Ester separated by a comma.

Theoretical expectations:
Let us consider first the case of one locus. Let us note 1-si and 1-c, the probability that a susceptible and a resistant homozygote survive exposure to insecticide, and 1 and 1-c, the probability that they survive in the absence of insecticide. Further, si represents the fitness decrease due to insecticide exposure, and c the fitness cost of resistance. In the Montpellier area, insecticide treatments are restricted to the coastal belt (between 0 and L kilometers from the sea). For one locus with two alleles, the fitness of each genotype can be written as follows:

For codominance (d = 0) and = 1, NAGYLAKI 1975 showed that a cline may be maintained given that k > /4, with k² = where is the standard deviation of parent-offspring distance measured along one dimension. Such clines cannot be characterized only by their slope, in contrast to numerous other cases (BARTON and GALE 1993 ), because they can be very asymmetric. However, NAGYLAKI 1975 (Equations 32–33) showed that they can be described by the maximum gene frequency and the gene frequency at the transition between the two environments and that the relative magnitude of selection vs. migration can be deduced from these characteristics.

The full analytical treatment in the case of two loci in a semi-infinite environment has not been performed, although SLATKIN 1975 worked out numerically the case of an infinite environment. When two loci with two alleles each are considered, fitness interactions between genes as well as their linkage must be considered. Additionally, linkage disequilibria are generated by migration that steepen each of the clines.

Descriptive analysis:
In order to test for the presence of frequency gradients at each resistance allele along the transect, data were fitted to descriptive cline models. Allelic distributions were fitted according to a scaled negative exponential. For instance, the frequencies of the four esterase alleles were modeled as follows:

• Ester¹: f1(x) = h1.e^-a1.x2

• Ester²: f2(x) = h2.(1 - h1).e^-a2.x2

• Ester⁴: f4(x) = h4.(1 - h1 - h2.(1 - h1)).e^-a4.x2

• Ester⁰: f0(x) = 1 - f1(x) -f2(x) - f4(x),

where a1, a2, a4, h1, h2, and h4 are estimated parameters. The phenotypic distributions were computed by using these allelic distributions and by assuming each locus at Hardy-Weinberg equilibrium. The phenotype was considered to be a three-state or seven-state random variable for the Ace.1 and Ester locus, respectively (see Table 1). The likelihood of a sample was computed from the phenotypic multinomial distribution.

Departure from Hardy-Weinberg proportions was tested in each population at the Ester locus by a likelihood ratio test. For an overall test, P values of each test were combined across populations using Fisher's method (MANLY 1985 ). At the Ace.1 locus, departure from Hardy-Weinberg cannot be evaluated because only three phenotypes are identified for three alleles. The presence of the Ace.1^RS allele creates an apparent excess of [RS] when only the two alleles Ace.1^R and Ace.1^S are considered. If Hardy-Weinberg proportions are assumed at this locus, Ace.1^RS frequencies can be computed from this apparent excess of [RS], and an additional cline of allele frequency can be fitted.

A linkage disequilibrium measure D = freq{S,O} - freq{S} x freq{O} was computed for each population and tested by an exact test on the contingency table ({S},{R}) x ({O},{E}) using the Genepop software (ver. 3.1a; RAYMOND and ROUSSET 1995 ). P values of each test were combined across populations using Fisher's method.

Simulations:
In order to estimate migration and selection, we used deterministic simulations to infer the allelic distribution at equilibrium because the analytical solution is intractable and requires the assumption of weak selection. One-dimension clines were simulated by a series of demes connected by migration as described in MALLET and BARTON 1989 . The migration distribution was reflected at one edge of the stepping stone to simulate a semi-infinite environment. The probability P of an individual in deme a migrating into deme a - t + i was calculated using a symmetric binomial distribution B(2t, 1/2) corrected by the reflecting condition when a < t.

The migration variance was measured by ²t/2, which is the variance of this distribution when a > t and where is the distance between demes. Selection coefficients were combined additively. The order of the processes was assumed to be reproduction-migration-selection, as should be the case for C. pipiens.

Migration and selection estimations:
The method of estimation is based on the principle that all resistance allele frequencies should be clinal, decreasing from south to north (Figure 1). As a consequence, the mixing of genotypes by migration from populations along these clines should create heterozygote deficiencies at each locus (Wahlund effect) and positive linkage disequilibrium between loci. This disequilibrium is predicted to be maximal at a medium distance from the coast (x L) where the most dissimilar genotypes are mixed. Migration and selection parameters were estimated conjointly such that the expected frequencies, computed using the simulation described above, and observed frequencies were as close as possible.

We focused our study on the differences between susceptible and resistance alleles within and between loci rather than on the transient polymorphism or allele replacements at each locus. For such a purpose, we pooled individuals carrying at least one resistance allele at each locus. The phenotype was considered therefore to be a four-state ({S,O}, {S,E}, {R,O} and {R,E}; see Table 1) random variable, and the likelihood of a sample was computed from its multinomial distribution. Eleven parameters are needed to describe the system. Among them, three can be estimated from external data: the recombination rate (r = 14.5%), the size of the treated area (L = 20 km), and the epistasis for resistance (zero). Furthermore, we assumed that epistasis for fitness costs was negligible. These estimations and assumptions allowed us to investigate the selection intensities (s_a, _a for the Ace.1 locus and s_e, _e, for the Ester locus) and the migration variance (²). To evaluate the influence of the dominance level on the estimation of the migration variance, three cases of dominance for both loci were considered: recessivity (d = -1), codominance (d = 0), and dominance (d = 1). The influence of the recombination rate was also investigated for codominance at both loci.

Model comparisons and tests:
Maximum likelihood estimates (MLE) of parameters were computed conjointly using the Metropolis algorithm adapted from N. H. BARTON (SZYMURA and BARTON 1986 ). G-tests were computed between related models and scaled to the dispersion of residual deviance (CRAWLEY 1993 ). The support limits of a particular parameter were defined as the range of values within two units of log-likelihood from the maximum (EDWARDS 1972 ).

	RESULTS

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Resistance allele frequencies and linkage disequilibria:
The frequencies of the different phenotypes combined at both loci are given in Table 2. The Hardy-Weinberg expectation was not rejected at the Ester locus (global test over populations P = 0.87). The linkage disequilibrium estimates between Ester and Ace.1 (D) and their corresponding P value are indicated in Table 2. A positive D is observed (Table 2D > 0, combined test across populations P = 5.10^-5) that peaks (4–6%) near the ecotone transition, which is consistent with a linkage disequilibrium created by migration.

Descriptive models:
At the Ace.1 locus, a clinal pattern is detected for both Ace.1^R and Ace.1^RS alleles (Table 3). The presence at high frequencies of the duplication (0.33 on the coast, Table 4) is thus strongly supported, and its cline explains well the pattern of apparent excess of heterozygotes in the transect. These two similar clines explain 88% of the total deviance at Ace.1 locus. However, the frequency of the duplicated allele Ace.1^RS is underestimated by assuming Hardy-Weinberg proportions because a heterozygote deficiency due to migration is expected. At the Ester locus, the model explains 77% of the total deviance. Significant and similar clinal patterns were found for all esterase resistance alleles, even for Ester², which is rare (Table 3 and Table 4). These results are consistent with the hypothesis that, for each locus, the selection pressures acting on resistance alleles are similar, that is, that the main differences in selection pressure are between susceptible and resistance alleles.

Migration-selection models:
The migration-selection models explain 92–93% of the total deviance (Table 5 and Figure 2). The maximum likelihood estimate of the parent-offspring standard deviation measured on one dimension is = 6.6 km.gen^- (support limits 4.8–8.7 km.gen^-1/2) when codominance is assumed at both loci. It should be underlined that the expected linkage disequilibrium does not peak at the ecotone transition, as in the case of an infinite environment (see SLATKIN 1975 and Figure 2A), but is shifted 8 km into the untreated area.

View larger version (18K): In this window In a new window Download PPT slide   Figure 2. Fitted and observed clines and linkage disequilibrium. (a) Linkage disequilibrium D (see text); (b) {S,O} frequency; (c) {R,O} frequency; (d) {S,E} frequency; (e) {R,E} frequency. Circles represent observed values and lines fitted values for codominance at both loci (see text).

Effect of the linkage: As in the case of a single locus, the ratio of selection intensities in treated and untreated areas (²) and the selection-migration ratio (k) depend only on the relative magnitude of selection vs. gene flow for each locus (data not shown). However, the linkage between the two loci (14.5%) has a noticeable effect on the Ester locus: in order to maintain the Ester cline at the same frequency in the absence of selection on Ace.1, selection (or k²) would have to be 26% higher. In contrast, to maintain the Ace.1 cline at the same level, k² need only be increased by 7% in the absence of selection on the Ester locus.

Dominance effect: The different hypotheses of dominance do not have an important effect neither on the estimation of migration variance ( range 6.6–7.1) nor on its support limits (Table 5). In fact, for given frequency gradients, the estimation depends mainly on the linkage disequilibrium pattern, as previously pointed out by BARTON 1982 .

Recombination effect: The estimate of the migration variance strongly depends on the recombination rate between Ace.1 and Ester: for the same migration variance, the closer the loci the higher the linkage disequilibrium. The recombination rate (r) of 14.5% was estimated between Ester and Ace.1 based on 503 individuals (T. LENORMAND, T. GUILLEMAUD, D. BOURQUET and M. RAYMOND, unpublished results). Figure 3 shows the joint support area for r and , assuming codominance at both loci. Support limits of are not affected by error measurements on the recombination rate (see Figure 3).

View larger version (12K): In this window In a new window Download PPT slide   Figure 3. Maximum likelihood estimates (bold line) of the standard deviation of parent-offspring distance ( in km·gen-1/2) as a function of the recombination rate (r) between the Ace.1 and Ester loci. Outer lines correspond to the support limits of . The small circle is the joint maximum likelihood of r and , and the dashed ellipse the joint support area for both r and .

Selection intensities: Estimations of selection intensities may not be as robust as the estimation of , because they depend on the assumptions of dominance. However, the different cases of dominance that were investigated are not equally likely (Table 5). In particular, models considering dominance at Ace locus (models A–C in Table 5) are ~20 times less likely than those that consider recessivity (models G–I). In contrast, for a given dominance on the Ace.1 locus, there are no noticeable differences between models considering different dominance levels on the Ester locus. When considering only the most likely models (E–I), selection intensity is likely to be ~0.12 on Ace.1 (s_a) and ~0.055 on Ester (s_e). For codominance at both loci (model E), these selection intensities give an estimate of the insecticide selection pressure (si 0.30 for Ace.1 and 0.16 for Ester) and of the intensities of the fitness costs (c 0.11 and 0.06 for Ace.1 and Ester, respectively).

	DISCUSSION

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

Validity of the assumptions:
The analysis of clinal patterns allowed us to infer in a single step the different parameters that are relevant to describe the dynamics of local adaptation, i.e., gene flow and selection coefficients in the different part of the environment for each locus. Additionally, the model developed permits us to explain 92% of the total deviance of the data in a quite economical manner. However, many simplifications were assumed, and external estimations were used for some parameters. We will discuss these points in turn.

Models of selection: We assumed that the different resistance alleles at each locus were subjected to the same selection pressure. This is probably not true since allele replacements were observed over the last 20 years (GUILLEMAUD et al. 1998 ). However, the rate of these replacements is quite low, at least for the Ester locus, and can be explained by fitness differences between resistance alleles that are much lower (1–2%) than those between resistance and susceptible alleles (GUILLEMAUD et al. 1998 ). Additionally, the fitted selection intensities may mainly represent those associated with the most common alleles. This is especially true for the Ester locus: among individuals that carry at least one resistance allele, only 15% lack the Ester⁴ allele. The situation may not be as clear for the Ace.1 locus, where both the Ace.1^R and the duplicate Ace.1^RS alleles are present in non-negligible frequencies and where the rate of allele replacement over time is not documented.

External estimation of parameters: We supposed that the summer clines were observed at migration-selection equilibrium. This is of course not exactly true, because selection intensities vary during the year. However, the high selection pressure and migration variance estimated are consistent with very rapid adjustments of frequencies. Moreover, frequencies, as well as selection intensities, are autocorrelated in time: adjustments to selection intensities require only limited changes in frequency.

We assumed a symmetric binomial migration distribution with a reflecting condition on the sea coast. Departures from this assumption could exist due to several factors. First, the migration distribution may be more leptokurtic. This may not strongly affect the peak of linkage disequilibrium at the ecotone transition (see MALLET et al. 1990 ), but further work is required to settle this issue. Second, it is possible that density variations may cause asymmetric flux of migrants. These variations in density may be caused by the control of mosquito populations in the treated area, which affects population sizes. However, the number of favorable larval breeding sites varies as well and would tend to compensate for this effect: the density of C. p. pipiens larval sites are associated with human activity, which is more important in the treated (urban and peri-urban areas) than in the untreated areas (countryside). Third, because the density of hosts (for blood feeding) and of larval sites is likely to vary along the transect, the hypothesis of a constant migration variance may be violated if most of the migration variance is due to foraging behavior (search for a blood meal and a site to lay eggs; REISEN et al. 1991 ). Fourth, the pattern of cytoplasmic incompatibilities caused by Wolbachia endosymbionts may cause local variation of the effective migration variance (MAGNIN et al. 1987 ). For these reasons, the effective migration estimated in the Montpellier area may not be considered as a general feature of the species.

Comparison with other estimates:
Direct measures of dispersal: Many studies have investigated the active dispersal of Culex species by mark-recapture experiments. Although none consider C. p. pipiens, plenty of data is available for C. p. quinquefasciatus, the tropical subspecies of C. pipiens (MATTINGLY et al. 1951 ), and other Culex species. These dispersal estimations are often biased toward low values either as a consequence of the small areas investigated (4 km, MORRIS et al. 1991 ; REISEN et al. 1991 ; 1.5 km, SCHREIBER et al. 1988 ) or of the absence of correction for dilution of sampling effort with distance (REISEN et al. 1991 , REISEN et al. 1992 ). In all these studies, some individuals were trapped close to the limit of the trapping grid. MORRIS et al. 1991 report a mean distance traveled (mdt) per day for three Culex species (0.73, 0.76, and 0.84 km for C. erraticus, C. nigripalpus, and C. salinarius, respectively), and SCHREIBER et al. 1988 report a mdt of 1.27 km after 36 hr for C. p. quinquefasciatus. Over a period of 12 days, REISEN et al. report an mdt between 0.6 and 1 km (1991) or 2 km (1992) for C. p. quinquefasciatus (which are both strongly biased; see above), and O'DONNELL et al. 1992 estimated an mdt of 6.8 km for C. annulirostris. Given an average of 8–10 days from adult emergence to the first oviposition (LOWE et al. 1973 ; SMITTLE et al. 1973 ; WEIDHAAS et al. 1973 ), these estimates are in agreement with ours (6.6 km·gen^-1/2 corresponds to an mdt of 0.66–0.82 km/day), although it appears that mark-recapture experiments should be performed on a larger scale to give more reliable estimates of dispersal for Culex species.

Migration-drift equilibrium: The relative magnitude of gene flow vs. drift has been evaluated in southern France by CHEVILLON et al. 1995 , using allozymic markers. Their study supports the hypothesis that migration outweighs drift in C. p. pipiens populations, which is in good agreement with our high estimate of migration variance. However, our estimate of allows us to disentangle drift and migration and to estimate average effective population densities. We reanalyzed the allozymic data from CHEVILLON et al. 1995 , using the method described in ROUSSET 1997 . Samples were collected in 31 localities distributed along the Mediterranean coast (southern France and northern Spain) and analyzed at five loci. A regression of F_ST/(1-F_ST) estimates computed for pairs of subpopulations on the logarithm of geographical distance was performed using the Genepop software (ver. 3.1a; RAYMOND and ROUSSET 1995 ). A significant isolation by distance was detected (Mantel test, P = 0.0123). The slope of this regression was used to estimate 1/(4 De.²), where De is the density of mosquitoes (ROUSSET 1997 ). The estimate of De.² is 16.3 individuals. Using our estimate of = 6.6 km·gen^- , the estimate of the density De is 0.37 individuals per km². This result is surprising when compared to the mosquito densities observed in the field during the breeding season (10⁴–10⁷ individuals/km²; REISEN et al. 1991 , REISEN et al. 1992 ; LINDQUIST et al. 1967 ). This strongly suggests that mosquito populations are heterogeneous in space, that they vary seasonally, and that they endure severe bottlenecks. A simple explanation of these results would be the low density of the founders in early spring. High mortality rates during overwintering (up to 80–90%; MINAR and RYBA 1971 ; SULAIMAN and SERVICE 1983 ) and between pupation and oviposition (90%; LOWE et al. 1973 ; WEIDHAAS et al. 1973 ) have been reported, and weather-related mortality is likely to occur in the spring due to a fall of temperature after warming periods. Finally, it is also possible that the De estimate is not totally accurate due to the different scales of the two approaches.

Selection intensities: We found that insecticide selection on Ace.1 locus (si 0.30) was higher than on Ester locus (si 0.16). This is consistent with the resistance ratio associated with these loci: the insensitive AChE1 confers a higher level of resistance. However, even if selection has been clearly associated with OP insecticides, it has never been measured in natural populations. The evaluation of fitness costs is even less straightforward because all fitness components can be influenced during both larval and adult stages. For example, larval development time, fecundity, susceptibility to parasites or predators, ability to blood feed, etc., can be modified by the presence of resistance genes (e.g., WOOD and BISHOP 1981 ; ROUSH and MCKENZIE 1987 ). For C. p. pipiens, some fitness costs on larval development time and female fecundity have been experimentally found in natural populations (BOURGUET 1996 ). The presence of the Ace.1^R allele increased the generation time by 4.3%. In an exponential growth phase of the population (during the spring), this difference gives an estimation of fitness cost (c) between 0.07 and 0.13 for an effective fecundity between five and 25 offspring. This estimation may be conservative, because only larval development time is taken into account and indicates that our estimates are not overestimated (c 0.11 and 0.06 for Ace.1 and Ester, respectively).

Conditions for the maintenance of Ace.1 and Ester clines:
Both clines maintain each other. Their concomitant presence makes the conditions for their maintenance less strict than if they were alone. However, this effect concerns mainly the least selected locus, Ester. We computed that the frequencies of Ester⁰ and Ace.1^S in coastal populations would be 0.074 and 0.013 higher, respectively, if each locus was considered independently from the other. Using the estimates of migration and selection provided by model E (codominance at both loci), the Ester cline would disappear if the width of the treated area (L) was reduced to 11 km. Similarly, the Ace.1 cline is not maintained when L < 7 km. In an infinite and uniform environment, the minimum size of a potential adaptive pocket would therefore be ~15 km for codominance at both loci.

Estimating gene flow from selected loci:
When estimating gene flow from selected loci, the selection pressure is taken explicitly into account. This situation presents different advantages. First, there is no need to formulate ad hoc hypotheses concerning neutrality of markers; second, few markers are needed; and third, predictable frequency patterns are expected and can be tested. However, this method requires that some genes be identified that are subjected to clear selection pressures and that some conclusion be made a priori concerning these selection pressures. Additionally, this method permits working at a restricted scale in time and space where assumptions of constant population sizes and homogeneity of space are the most reliable. In particular, it is possible to take explicitly into account specific features of the environment (e.g., presence of geographic barriers), if needed. The drawback is that such estimates can hardly be representative of other environmental conditions because they are not averaged over a long period of time and over large geographic areas. However, they provide an "instantaneous" measure of dispersal that is the most pertinent for the area, the period, and the scale considered, especially in the case of recent local adaptation. In this respect, these measures may be comparable to mark-recapture estimates. However, direct measures of dispersal do not provide estimates of effective gene flow and may miss long-distance migrants because individuals are often trapped at the limit of the trapping grid.

We have estimated gene flow using two selected genetic markers. This estimation is an essential step in understanding the dynamics of selected genes when selection pressures vary in space and time. We have focused on the selection-migration equilibrium on a local scale, considering gene flow as a "constraining force" reducing the potential for local adaptation. However, on a much wider scale, gene flow is also responsible for the spread of resistance alleles across the species range mainly by passive migration (QIAO and RAYMOND 1995 ; GUILLEMAUD et al. 1996 and reference therein). At these two scales, gene flow may not be comparable, playing either a conservative or a creative role.

	FOOTNOTES

¹ Present address: Station de Recherche de lutte biologique, INRA La Minière, 78285 Guyancourt Cedex, France.

	ACKNOWLEDGMENTS

We are very grateful to C. CHEVILLON, N. PASTEUR, F. ROUSSET and J. BRITTON-DAVIDIAN for their helpful comments and discussion about the manuscript. This work was financed by Groupement de Recherche 1105 du programme Environnement, Vie et Sociétés du Centre National de la Recherche Scientifique, the Région Languedoc-Roussillon (no. 963223) and ACC SV3 (no. 9503037). T.L. was supported by an ASC from Institut National de la Recherche Agronomique. Contribution 98.062 of the Institut des Sciences de l'Evolution de Montpellier (UMR CNRS 5554).

Manuscript received September 16, 1997; Accepted for publication March 9, 1998.

	LITERATURE CITED

TOP ABSTRACT MATERIALS AND METHODS RESULTS DISCUSSION LITERATURE CITED

ANONYMOUS, 1990–1995 Rapport d'activité technique et scientifique. Entente Interdépartementale pour la Démoustication du Littoral Méditerranéen, Montpellier, France.

BARTON, N. H., 1982  The structure of the hybrid zone in Uroderma bilobatum (Chiroptera: Phyllostomatidae). Evolution 36:863-866.

BARTON, N. H., and K. S. GALE, 1993 Genetic analysis of hybrid zones, pp. 13–45 in Hybrid Zones and the Evolutionnary Process, edited by R. G. HARRISON. Oxford University Press, Oxford, UK.

BARTON, N. H. and G. M. HEWITT, 1985  Analysis of hybrid zones. Ann. Rev. Ecol. Syst. 16:113-148.

BOURGUET, D., 1996 Mécanismes biochimiques, plasticité et évolution de la dominance. Ph.D. Thesis. Université Montpellier II, Sciences et Techniques du Languedoc, Montpellier, France.

BOURGUET, D., M. RAYMOND, D. FOURNIER, D. MALCOLM, and P. TOUTANT et al., 1996a  Existence of two acetylcholinesterases in the mosquito Culex pipiens complex (Diptera: Culicidae). J. Neurochem. 67:2115-2123[Medline].

BOURGUET, D., R. CAPELA, and M. RAYMOND, 1996b  An insensitive acetylcholinesterase in Culex pipiens L. mosquitoes from Portugal. J. Econom. Entomol. 89:1060-1066[Medline].

BOURGUET, D., M. RAYMOND, J. BISSET, N. PASTEUR, and M. ARPAGAUS, 1996c  Duplication of the Ace.1 locus in Culex pipiens from the Caribbean. Biochem. Genet. 34:351-362[Medline].

BOURGUET, D., M. PROUT, and M. RAYMOND, 1996d  Dominance of insecticide resistance presents a plastic response. Genetics 143:407-416[Abstract].

BOURGUET, D., N. PASTEUR, J. BISSET, and M. RAYMOND, 1996e  Determination of Ace.1 genotypes in single mosquitoes: toward an ecumenical biochemical test. Pest. Biochem. Physiol. 55:122-128[Medline].

CHEVILLON, C., N. PASTEUR, M. MARQUINE, D. HEYSE, and M. RAYMOND, 1995  Population structure and dynamics of selected genes in the mosquito Culex pipiens.. Evolution 49:997-1007.

CHEVILLON, C., D. BOURGUET, F. ROUSSET, N. PASTEUR, and M. RAYMOND, 1997  Pleiotropy of adaptive changes in populations: comparisons among insecticide resistance genes in Culex pipiens.. Genet. Res. 70:195-204[Medline].

CRAWLEY, M. J., 1993 GLIM for Ecologists. Blackwell Scientific Publications, Oxford, UK.

DE STORDEUR, E., 1976  Esterases in the mosquito Culex pipiens pipiens L.: formal genetics and polymorphism of adult esterases. Biochem. Genet. 14:481-493[Medline].

EDWARDS, A. W. F., 1972 Likelihood. Cambridge University Press, Cambridge.

ENDLER, J. A., 1977 Geographic variation, speciation and clines. Princeton University Press, Princeton, NJ.

FELSENSTEIN, J., 1976  The theoretical population genetics of variable selection and migration. Annu. Rev. Genet. 10:253-280[Medline].

FISHER, R. A., 1937  The wave of advance of advantageous genes. Ann. Eugen. 7:353-369.

GUILLEMAUD, T., S. ROOKER, N. PASTEUR, and M. RAYMOND, 1996  Testing the unique amplification event and the worldwide migration hypothesis of insecticide resistance genes with sequence data. Heredity 77:535-543.

GUILLEMAUD, T., N. MAKATE, M. RAYMOND, B. HIRST, and A. CALLAGHAN, 1997  Esterase gene amplification in Culex pipiens.. Insect Mol. Biol. 6:319-327[Medline].

GUILLEMAUD, T., T. LENORMAND, D. BOURGUET, C. CHEVILLON, and N. PASTEUR et al., 1998  Evolution of resistance in Culex pipiens: allele replacement and changing environment. Evolution. in press.

LANDE, R., 1985  The fixation of chromosomal rearrangements in a subdivided population with local extinction and colonization. Heredity 54:323-332.

LINDQUIST, A. W., T. IKESHOJI, B. GRAB, B. DE MEILLON and Z. H. KHAN, 1967 Dispersion studies of Culex pipiens fatigans tagged with ³²P in the Kemmendine area of Rangoon, Burma. Bull. Org. mond. Santé 36: 21–37.

LOWE, R. E., H. R. FORD, B. J. SMITTLE, and D. E. WEIDHAAS, 1973  Reproductive behavior of Culex pipiens quinquefasciatus released in natural populations. Mosq. News 33:221-227.

MAGNIN, M., N. PASTEUR, and M. RAYMOND, 1987  Multiple incompatibilities within populations of Culex pipiens L. in southern France. Genetica 74:125-130[Medline].

MALCOLM, C., D. BOURGUET, A. ASCOLILLO, S. J. ROOKER, and C. F. GARVEY et al., 1998  A sex-linked Ace gene, not linked to insensitive acetylcholinesterase-mediated insecticide resistance in Culex pipiens.. Insect. Mol. Biol. 7:107-120[Medline].

MALLET, J. and N. H. BARTON, 1989  Inference from clines stabilized by frequency-dependent selection. Genetics 122:967-976[Abstract/Free Full Text].

MALLET, J., N. H. BARTON, G. M. LAMAS, J. C. SANTISTEBAN, and M. M. MUEDAS et al., 1990  Estimates of selection and gene flow from measure of clines width and linkage disequilibrium in Heliconius hybrid zones. Genetics 124:921-936[Abstract].

MANLY, B. F. J., 1985 The Statistics of Natural Selection. Chapman and Hall Ltd., London.

MATTINGLY, P. F., L. E. ROZEMBOOM, K. L. KNIGHT, H. LAVEN, and S. R. DRUMMOND et al., 1951  The Culex pipiens complex. Trans. R. Ent. Soc. Lond. 102:331-342.

MINAR, J. and J. RYBA, 1971  Experimental studies on overwintering conditions of mosquitoes. Folia parasitologica (Praha) 18:255-259.

MORRIS, C. D., V. L. LARSON, and L. P. LOUNIBOS, 1991  Measuring mosquito dispersal for control programs. J. Am. Mosq. Cont. Assoc. 7:608-615[Medline].

NAGYLAKI, T., 1975  Conditions for the existence of clines. Genetics 80:595-615[Abstract].

O'DONNELL, M. S., G. BERRY, T. CARVAN, and J. H. BRYAN, 1992  Dispersal of adult females of Culex annulirostris in Griffith, New South Wales, Australia. J. Am. Mosq. Cont. Assoc. 8:159-165[Medline].

PASTEUR, N., A. ISEKI, and P. GEORGHIOU, 1981a  Genetic and biochemical studies of the highly active esterases A and B associated with organophosphate resistance in mosquitoes of the Culex pipiens complex. Biochem. Genet. 19:909-919[Medline].

PASTEUR, N., G. SINEGRE, and A. GABINAUD, 1981b  Est-2 and Est-3 polymorphisms in Culex pipiens L. from southern France in relation to organophosphate resistance. Biochem. Genet. 19:499-508[Medline].

PASTEUR, N., G. PASTEUR, F. BONHOMME and J. BRITTON-DAVIDIAN, 1988 Practical Isozyme Genetics. Ellis Horwood Ltd., Chichester, U.K.

POIRIÉ, M., 1991 Mécanisme de détoxification impliqués dans la résistance aux insecticides: glutathion transférases et estérases non spécifiques. Ph.D. Thesis. Université Montpellier II, Sciences et Techniques du Languedoc, Montpellier, France.

POIRIÉ, M., M. RAYMOND, and M. PASTEUR, 1992  Identification of two distinct amplifications of the esterase B locus in Culex pipiens (L.) mosquitoes from Mediterranean countries. Biochem. Genet. 30:13-26[Medline].

QIAO, C.-L. and M. RAYMOND, 1995  The same esterase B1 haplotype is amplified in insecticide resistant mosquitoes of the Culex pipiens complex from the Americas and China. Heredity 74:339-345.

RAYMOND, M. and M. MARQUINE, 1994  Evolution of insecticide resistance in Culex pipiens populations: the Corsican paradox. J. Evol. Biol. 7:315-337.

RAYMOND, M. and F. ROUSSET, 1995  Genepop version 1.2: population genetics software for exact tests and ecumenicism. J. Hered. 86:248-249[Free Full Text].

RAYMOND, M., D. FOURNIER, J.-M. BRIDE, A. CUANY, and J. BERGÉ et al., 1986  Identification of resistance mechanisms in Culex pipiens (Diptera: Culicidae) from southern France: insensitive acetylcholinesterase and detoxifying oxidases. J. Econom. Entomol. 79:1452-1458[Medline].

RAYMOND, M., N. PASTEUR, and G. P. GEORGHIOU, 1987  Inheritance of chlorpyrifos resistance in Culex pipiens L. (Diptera: Culicidae) and estimation of the number of genes involved. Heredity 58:351-356.

RAYMOND, M., D. HECKEL, and J. G. SCOTT, 1989  Interaction between pesticide genes: model and experiment. Genetics 123:543-551[Abstract/Free Full Text].

REISEN, W. K., M. M. MILBY, R. P. MEYER, A. R. PFUNTNER, and J. SPOEHEL et al., 1991  Mark-release-recapture studies with Culex mosquitoes (Diptera: Culicidae) in southern California. J. Med. Entomol. 28:357-371[Medline].

REISEN, W. K., M. M. MILBY, and R. P. MEYER, 1992  Population dynamics of adult Culex mosquitoes (Diptera: Culicidae) along the Kern river, Kern county, California, in 1990. J. Med. Entomol. 29:531-542[Medline].

RIVET, Y., M. RAYMOND, J.-A. RIOUX, A. DELABRE, and N. PASTEUR, 1994  Resistance monitoring in Culex pipiens L. from central eastern France. J. Med. Entomol. 31:231-239[Medline].

ROOKER, S., T. GUILLEMAUD, J. BERGÉ, N. PASTEUR, and M. RAYMOND, 1996  Coamplification of esterase A and B genes as a single unit in the mosquito Culex pipiens.. Heredity 77:555-561.

ROUSH, R. T. and J. A. MCKENZIE, 1987  Ecological studies of insecticide and acaricide resistance. Ann. Rev. Entomol. 32:361-380[Medline].

ROUSSET, F., 1997  Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics 145:1219-1228[Abstract].

SCHREIBER, E. T., M. S. MULLA, J. D. CHANEY, and M. S. DHILLON, 1988  Dispersal of Culex quinquefasciatus from a dairy in southern California. J. Am. Mosq. Cont. Assoc. 4:300-304[Medline].

SEVERINI, C., R. ROMI, M. MARINUCCI, and M. RAYMOND, 1993  Mechanisms of insecticide resistance in field populations of Culex pipiens from Italy. J. Am. Mosq. Cont. Assoc. 9:164-168[Medline].

SITES, J. W., N. H. BARTON, and K. M. REED, 1995  The genetic structure of a hybrid zone between two chromosome races of the Sceloporus grammicus complex (Sauria, Phrynosomatidae) in central Mexico. Evolution 49:9-36.

SLATKIN, M., 1975  Gene flow and selection in a two locus system. Genetics 81:787-802[Abstract/Free Full Text].

SLATKIN, M., 1987  Gene flow and the geographic structure of natural populations. Science 236:787-792[Abstract/Free Full Text].

SLATKIN, M. and N. H. BARTON, 1989  A comparison of three indirect methods for estimating average levels of gene flow. Evolution 43:1349-1368.

SMITTLE, B. J., R. E. LOWE, H. R. FORD, and D. E. WEIDHAAS, 1973  Techniques for ³²P labeling and essay of egg rafts from field-collected Culex pipiens quinquefasciatus Say. Mosq. News 33:215-220.

SULAIMAN, S. and M. W. SERVICE, 1983  Studies on hibernating populations of the mosquito Culex pipiens L. in southern and northern England. J. Nat. Hist. 17:849-857.

SZYMURA, J. M. and N. H. BARTON, 1986  Genetic analysis of a hybrid zone between the fire-bellied toads, Bombina bombina and Bombina variegata, near Cracow in southern Poland. Evolution 40:1141-1159.

WEIDHAAS, D. E., B. J. SMITTLE, R. S. PATTERSON, R. E. LOWE, and C. S. LOFGREN, 1973  Survival, reproductive capacity and migration of adult Culex pipiens quinquefasciatus Say. Mosq. News 33:83-87.

WOOD, R. J., and J. A. BISHOP, 1981 Insecticide resistance: populations and evolution, pp. 97–127 in Genetic Consequences of Man-made Changes, edited by J. A. BISHOP, and L. M. COOK. Academic Press, London.

This article has been cited by other articles:

				C. F. BOSIO, L. C. HARRINGTON, J. W. JONES, R. SITHIPRASASNA, D. E. NORRIS, and T. W. SCOTT GENETIC STRUCTURE OF AEDES AEGYPTI POPULATIONS IN THAILAND USING MITOCHONDRIAL DNA Am J Trop Med Hyg, April 1, 2005; 72(4): 434 - 442. [Abstract] [Full Text] [PDF]

				M. H. Kohn, H.-J. Pelz, and R. K. Wayne Locus-Specific Genetic Differentiation at Rw Among Warfarin-Resistant Rat (Rattus norvegicus) Populations Genetics, July 1, 2003; 164(3): 1055 - 1070. [Abstract] [Full Text] [PDF]

				C. Berticat, M. Marquine, M. Raymond, and C. Chevillon Recombination Between Two Amplified Esterase Alleles in Culex pipiens J. Hered., July 1, 2001; 92(4): 349 - 351. [Abstract] [Full Text] [PDF]