Mosquito strains
Three mosquito laboratory strains were used in this study: SLAB [31], SR [32] and SRQ (this study). SLAB is fixed for a single-copy susceptible allele (SSLAB, isolated in California, C. quinquefasciatus). SR is fixed for R1 [24], a resistance allele isolated from Southern France and found in C. pipiens s.s. all over Europe and around the Mediterranean Sea [23, 33]. SRQ is fixed for R2 [24], a resistance allele found worldwide in C. quinquefasciatus [23] and isolated from a population from Martinique Island. The two resistance alleles were introgressed into the genetic background of the SLAB strain through at least 15 rounds of back-crossing. All strains thus share the same genetic background (> 99%) and differ from each other almost only in their ace-1 locus (although recombination around the ace-1 gene is not complete, most of the background effects would be eliminated).
All strains were regularly checked for contamination: DNA was extracted from pools of first-instar larvae (~ 200 individuals per pool) and molecular tests, specific for each ace-1 allele (detailed below), were used to check the homogeneity of each strain.
Genotyping
The various strains can be easily distinguished using a single PCR and different restriction fragment length polymorphisms (RFLPs). After DNA extraction, following the protocol in [34], a ~ 600-bp fragment of the ace-1 gene, including intron 2 and most of exon 3 (with the resistance G119S mutation), was amplified using two generalist primers, Intron2dir1 and CpEx3rev, according to [35].
Susceptible vs. resistant
The G119S mutation creates an AluI restriction site [33] so that three genotypes can be distinguished (AluI RFLP test): susceptible homozygote (SS; one fragment, 597 bp), resistant homozygote (RR; two fragments, 496 and 101 bp) and heterozygote (RS; three fragments, 597, 496 and 101 bp); 5 µl of the PCR product was incubated for 2 h at 37 °C.
R1 vs. R2
The two different resistance alleles can be further distinguished by taking advantage of another single-nucleotide polymorphism (SNP) between R1 and R2 creating a Bfa1 restriction site in R2 (Additional File 1). This second RFLP test (Bfa1 RFLP test) distinguishes three genotypes, the homozygotes R1R1 (one fragment, 597 bp) and R2R2 (three fragments, 73, 132 and 392 bp) and the heterozygotes R1R2 (four fragments, 597, 73, 132 and 392 bp); 5 µl of the PCR product was incubated for 2 h at 37 °C.
Gene copy number quantification
ace-1 gene copy number was estimated for ten individuals of each resistant strain using quantitative real-time PCR (qRT-PCR). Two individuals from the SLAB-susceptible strain were also used as controls. After DNA extraction, we dispensed 250 ng of genomic DNA and 1.5 μl of reaction mixture containing specific primers, each at a concentration of 0.8 μM and 0.75 μl of Master Mix (LightCycler 480 SYBR Green I Master, Roche), into the wells of a 384-well plate with a Labcyte Echo525 dispenser. We performed qPCR as follows: activation at 95 °C for 8 min followed by 45 cycles of 95 °C for 4 s, 67 °C for 13 s and 72 °C for 19 s. Melting curves were generated by a post-amplification melting step between 70 °C and 95 °C for Tm analysis. All quantifications were replicated three times for each DNA template. Two loci were amplified for each individual: ace-1 (primers: ‘Culexace1univdir3’ AGA AGG TGG ACG CAT GGA TG; ‘Culexace1univrev3’ ATC TGG ACG CAG GAG TTG G) and ace-2, a locus known to be in a single copy in these species (primers: ‘acequantidir’ GCA GCA CCA GTC CAA GG; ‘acequantirev’ CTT CAC GGC CGT TCA AGT AG) [36]. ace-1 over ace-2 copy-number ratios were determined by the advanced quantification method (LightCycler 480 software v.1.5.0). Standard reference curves were constructed with tenfold dilutions of a PCR product previously amplified with specific primers for each locus from SLAB DNA.
Phenotyping
Protein activity
We measured acetylcholinesterase (AChE1) activity for 48 individuals of each resistance strain, using spectrophotometry [37]. Adult mosquitoes were decapitated, and each head was individually homogenized in 400 μl of a phosphate buffer (0.25 M, pH7) supplemented with 1% Triton X-100. Homogenates were centrifuged (9.3 g for 3 min), and 100 μl of the supernatant was dispensed into each of two wells of a 96-well microtitration plate. We added 10 μl of propoxur, a carbamate insecticide, at 10−3 M and 10−1 M (diluted in ethanol) into the first and second well, respectively. The plate was incubated for 15 min at room temperature. We then added 100 μl of substrate solution (25 mM sodium phosphate, pH 7.0, 0.2 mM DTNB, 0.35 mM sodium bicarbonate, 2.5 mM acetylthiocholine) to each well. AChE1 activity was estimated by measuring the change in optical density following the cleavage of acetylthiocholine, as described by [38]. Optical density at 412 nm was recorded every minute for 15 min with an EL 800 microplate reader (Bio-Tek Instruments, Inc.). The mean slope of each reaction was calculated with KCjunior v1.41.4 analysis software (Bio-Tek Instruments, Inc.) and was used as a measurement of AChE1 activity in each well. Individual AChE1 activity was computed as the average activity between the two wells. To avoid any block or sex confounding effects, individuals from both sexes and the two strains were evenly distributed in the plates.
Resistance level and bioassays
We used bioassays to assess the three strains’ resistance to an OP insecticide, temephos (PESTANAL®,96% purity). We incubated 20 late third-instar larvae for 24 h at 27 °C ± 2 °C in plastic cups containing 99 ml of distilled water to which we added 1 ml of insecticide solution at the required concentration (1 ml of ethanol in controls). Four replicates were performed for each concentration (from 0 to 0.07 g\(\upmu\).ml−1 see Additional File 2 for the complete dataset). Larval mortality was recorded after 24 h of exposure. We used the BioRssay R package (v.1.0.0 [39], https://CRAN.R-project.org/package=BioRssay) to analyze the dose-mortality responses of the different ace-1 alleles and calculate the LD50 of the different strains, i.e. the lethal dose for 50% of the sample.
Experimental evolution in population cages
Population cages were used to set up a competition experiment between the two resistance alleles in absence of insecticides. R1R1 and R2R2 individuals were crossed, and the resulting F1 (100% R1R2 individuals) was reared until adulthood under standard conditions (25 °C, > 60% humidity, 12:12 h light:dark). Adults were released into a new cage to mate freely and reproduce. Their offspring were raised and released in new cages to ensure discrete generations. The process was repeated 11 times (i.e. 11 generations) with three independent cages (i.e. replicates). Almost each generation, and for each cage, about a hundred second-instar larvae were genotyped using the Bfa1 RFLP test (see above) to measure the frequency of each genotype (R1R1, R1R2, R2R2). Allelic frequencies were then computed from genotypic frequencies.
We estimated the relative fitness of the various genotypes (R1R1, R1R2 and R2R2) using a deterministic genetic model (reproduction and selection, 11 cycles, no drift). The model was adjusted to the data and optimized using a maximum-likelihood approach as in Milesi et al. [5, 23]. For the reproduction step, the frequency of each genotype in the larvae of generation i was computed from the allelic frequencies (p) in the gametes of the previous generation, assuming panmixia (Eq. 1):
$$\begin{array}{l} f ( {R_{1} R_{1} }) = p_{1}^{ 2} \\ f ( {R_{1} R_{2} } ) = 2 \times p_{1} \times p_{2} \\ f ( {R_{2} R_{2} } ) = p_{2}^{ 2}\end{array}$$
(1)
For each genotype g selection was then computed between larval and adult stages of generation i using the following genotype fitness: wR1R1 = 1, wR1R2 = 1 + h.s and wR2R2 = 1 + s, with h the dominance coefficient and s the selection coefficient, both varying between − 1 and 1 (Eq. 2):
$$f^{\prime}_{gi} = \frac{{f_{gi} \times W_{g} }}{{\sum (f_{gi} \times W_{g} )}}$$
(2)
The genotypic frequencies after selection were used to calculate the allelic frequencies in the gametes produced by the surviving adults (Eq. 3).
$$\begin{array}{l} {p}_{1}^{^{\prime}}=f\left({R}_{1}{R}_{1}\right)+\frac{f\left({R}_{1}{R}_{2}\right)}{2} \\ {p}_{2}^{^{\prime}}=(1-{p}_{1}^{^{\prime}})\end{array}$$
(3)
The first run of 100,000 simulations was used to explore the parameter space and provide the likelihood profile associated with different random pairs of h and s values (Eq. 4):
$$L = \mathop \sum \limits_{g} \mathop \sum \limits_{i} \left( {n_{gi} \times \ln \left( {f_{gi} } \right)} \right)$$
(4)
where n is the number of individuals of genotype g observed in the cages at generation i, and f is the frequency of the genotype g at generation i calculated for a given pair of h and s values using our deterministic genetic model. One million additional simulations were run with parameter ranges more limited around the maximal likelihood h and s pair to precisely estimate the coefficients and their support limits (rough equivalents to 95% confidence intervals), defined as h and s maximal and minimal values, resulting in a likelihood equal to the maximum likelihood minus 1.96, as in [5].
Statistical analyses
All the statistical analyses were conducted using the R software (R Core Team, https://www.r-project.org/):
We used the following linear model to compare ace-1 copy number between the various strains:
$${\gamma }_{ij}= \mu +{\alpha }_{j}+{\varepsilon }_{ij }\quad (\mathrm{mod}. 1)$$
with \({\gamma }_{ij}\) the number of copies of the ace-1 gene in replicate i from strain j, \(\mu\) the population mean, \(\alpha\) is the fixed effect of strain j (SLAB, SR or SRQ) and \({\varepsilon }_{ij}\) the error term following a normal distribution \(\mathcal{N}(0, 1)\).
We used the following linear model to test the significance of the difference in AChE1 activity between the SR and SRQ strains:
$${\gamma }_{ijkl}= \mu +{\alpha }_{j}+{\beta }_{k}+{\delta }_{l}+{\varepsilon }_{ijkl} \quad (\mathrm{mod}. 2)$$
with \({\gamma }_{ijkl}\) the AChE1 activity for individual i of strain j and sex k measured in plate l, \(\mu\) the population mean and \(\alpha\) the fixed effect of strain j (SR or SRQ). \(\beta\) and \(\delta\) are control for the fixed effects of sex k and plate l, respectively, and \({\varepsilon }_{ijkl}\) is the error term following a normal distribution \(\mathcal{N}(0, 1)\).
For both models, the significance of the various terms was tested using likelihood ratio tests (LRTs) comparing the full model with a model without the tested effect (‘anova’ function, R [40]). For both models, we also confirmed the absence of significant heteroskedasticity (‘bptest’ function, ‘lmtest’ R package [41]) and that the models’ residuals followed a normal distribution (‘shapiro.test’ function, ‘stats’ R package).
Finally, we used binomial proportion tests (‘prop.test’ function, ‘stats’ R package) to assess whether the allele frequencies at the end of the experimental evolution in cages (i.e. after 11 generations) differed from initial frequencies of 0.5 (100% R1R2 individuals).