Two Aedes aegypti laboratory strains were compared for the purpose of this study: the Bora Bora reference strain, known to be susceptible to most insecticides, and a strain artificially selected for resistance to a decomposed tree leaf litter showing a high toxicity for mosquito larvae. This leaf litter had been collected in a mosquito pond in Eastern France three months after treatment for mosquito control, and has been proved to contain Bacillus thuringiensis subsp. israelensis (Bti) spores from commercial origin . The use of this Bti-contaminated leaf litter in the selection experiments allowed mimicking the evolution of resistance to Bti in a situation close to field conditions. The susceptible and resistant strains were separated by 18 generations of selection (and by 20 generations in total), with a strong bottleneck at generation 10 (see Additional File 1 for more details on the demographic history of the resistant strain). The selection experiment and the bioassays implemented to monitor the evolution of resistance are described in . At each generation, the lethal dose for 50% of the individuals after a 24 h-exposure (24 h-LD50) was determined for each strain using the Probit software . After 18 generations of selection, the resistance ratios RR of the resistant strain (i.e., the ratio between the 24 h-LD50 values for the resistant and the susceptible strains, respectively) were 3.4-fold, 30.2-fold, 13.7-fold, 6.3-fold and 3-fold for the toxic leaf litter, Cry4A, Cry4B, Cry11A and Cyt1A toxins, respectively.
DNA extraction and genome scan
The genomic DNA used for all subsequent molecular work was extracted from fresh fourth-instar mosquito larvae using the Qiagen DNeasy Tissue Kit and protocol (Qiagen). Prior to extraction, the larvae midgut was removed carefully to avoid bacterial contamination.
The classical protocol of the Diversity Array Technology (DArT)  was slightly modified so as to provide hundreds of good-quality markers scattered in the genome of Ae. aegypti and possibly associated to gene-rich regions . Briefly, in a first step, genomic DNA was digested with restriction enzyme Bsp1286I and a specific adaptor was ligated to compatible ends. Restriction fragments including a particular Ae. aegypti MITE called Pony were PCR-amplified using a primer annealing to the adaptor sequence and a primer complementary to a conserved motif of the Pony element. PCR products obtained for all individuals of the two strains were pooled together and cloned to construct a DArT library containing a total of 6144 MITE-based clones. In a second step, a labelled target produced for each Ae. aegypti individual as described in the first step was hybridized to the library fragments spotted on a glass slide in order to reveal the polymorphic ones. Details of the protocol, in particular the adaptor and primer sequences used and the reproducibility rates, can be found in .
Identification and sequencing of outlier loci potentially under selection
For each DArT marker obtained, allelic frequencies were estimated with the Bayesian method with non-uniform prior distribution  implemented in AFLP-SURV 1.0 . Among those markers, we tracked those for which alternative phenotypes (fragment presence/absence) were fixed or nearly fixed in the two strains (for example, fragment present/absent for all individuals or all individuals except one). Our assumption was that such a pattern of extreme inter-strain genetic differentiation could be explained by the spread, in the resistant strain, of an advantageous allele initially present at low or intermediate frequency in the susceptible strain (standing variation), from which it is eventually purged by genetic drift and/or because it is slightly deleterious.
For 70 such markers, bacterial cultures were sent to Genome Express®
http://www.genome-express.com for insert amplification and sequencing with M13 forward and M13 reverse primers. Raw sequence files were edited with BioEdit 7.0.9  and purged from Pony and primer sequences. The obtained sequenced were blasted against the full genomic sequence of Ae. aegypti (consisting of 4758 supercontigs and available at http://aaegypti.vectorbase.org/GetData/Downloads?type=Genome).
Identification of candidate genes
Although mechanisms of resistance to Bti are still unknown in dipterans, resistances to several Cry toxins have been intensively studied, especially in lepidopteran pests resistant to transgenic crops expressing Bacillus thuringiensis Cry toxins genes [20, 30]. Because these toxins share similar three-dimensional structures, similar modes of action and resistance mechanisms can be expected between lepidopteran and dipteran insects [16, 61]. We therefore considered as candidate genes for Bti resistance those belonging to families previously proved to be involved in Cry resistance [20, 30, 31]. To this list, we added genes potentially implicated in activation of Bti toxins (aminopeptidases, e.g. ; and trypsins and chymotrypsins, e.g. ), in toxin binding (alkaline phosphatases, e.g. [49, 63]; aminopeptidases, e.g. ; cadherins, e.g. ; galactosidases and glycosyltransferases, e.g. [20, 65]); or immune defense (mitogen-activated protein kinases, e.g. ). A keyword search was conducted in the VectorBase database http://aaegypti.vectorbase.org/index.php and a total of 160 candidates located on 98 different supercontigs were identified out of the 15,419 putative genes (16789 transcripts in total) referenced in the Aedes aegypti genome.
Cadherin (CAD) and leucine aminopeptidase (LAP) gene sequencing
The complete genomic sequences of the CAD and LAP genes (VectorBase Gene IDs AAEL001196 and AAEL001649, respectively) were downloaded from the VectorBase website to help design sequencing primers (Additional file 4) with the software package Lasergene 7.2 (DNASTAR Inc.). The sequencing strategy for the CAD gene targeted exon 5 and more specifically the membrane-proximal subdomains (subdomains 4 to 6) of the protein which are the preferential binding sites of Bti Cry toxins in lepidopterans (Figure 1). For the LAP gene, three different primer pairs were selected to amplify the two main exons.
PCR amplifications were conducted for each gene in a 25-μl total volume with 2 mM MgCl2, 0.1 mM of each dNTP (Roche), 0.2 μM of each primer, 5 μg of BSA, 0.6 U of AmpliTaq Gold DNA polymerase (Applied Biosystems) and 10-30 ng of DNA. The PCR program included an initial 10-min denaturation step at 95°C; 40 cycles of denaturation at 95°C for 45s, annealing at the optimal temperature indicated in Additional file 4 for 45s and elongation at 72°C for 60s; followed by a final extension step at 72°C for 5 min. PCR products were purified with the QIAquick PCR purification kit (Qiagen) and sequencing reactions were performed in both directions using the amplification primers and the BigDye Terminator Cycle Sequencing Kit 3.1 (Applied Biosystems), following the manufacturer's indications. Fluorescently labelled sequencing products were run on an ABI PRISM 3100 capillary DNA sequencer (Applied Biosystems) and sequences were analyzed with SeqMan Pro 7.1.0 (DNASTAR Inc.). Overall, we obtained sequences for 11 and 12 individuals of the susceptible and resistant strains, respectively.
Sequence analysis and neutrality tests
The software DnaSP 5.0  was used to infer haplotype phase and to assess a variety of genetic diversity and differentiation parameters (e.g., nucleotide diversity p, haplotype diversity Hd, number of segregating sites S, Fst, etc.) for each gene. Several statistics were also calculated based on the observed polymorphism data to test for deviation from neutral evolution in the resistant strain, including Tajima's D , Fu and Li's D* and F* , and Fay and Wu's H . To assess whether these statistics significantly departed from a neutral scenario of evolution given the known demographic history of the resistant strain, we performed coalescent simulations using the program ms . This program generates random independent samples according to a Wright-Fisher neutral model allowing population size changes in the past. For each gene, the mutation rate μ was estimated from the per-locus mutation parameter ? observed for the susceptible strain (θ = 4Neμ and Ne = 6000) and used as the starting value for the simulations (that is, as the value at present). Then 1000 neutral samples consisting of 24 haplotypes were simulated based on the known demographic history of the resistant strain (Additional file 1).
RNA extraction and gene expression analyses
For the two candidate genes, real-time reverse-transcription PCR (RT-PCR) analyses were performed on three biological replicates for each strain, with each replicate consisting of 30 larvae reared in standard insectary conditions up to the fourth-instar stage (5 days). Total RNAs were extracted using TRIzol (Invitrogen) following the manufacturer's instructions and their quality was assessed with a 2100 Bioanalyzer (Agilent) after DNase I (Invitrogen) treatment. Four micrograms of total RNA were digested with DNase I (Invitrogen) and then used for first-strand cDNA synthesis with SuperScript III (Invitrogen) reverse transcriptase and oligo-dT20 primers for 60 min at 50°C, according to the manufacturer's instructions. Real-time RT-PCR reactions were performed on an iQ5 system (Bio-Rad) in a 25-μL total reaction volume with 0.3 μM of each primer, 12.5 μL of iQ SYBR Green supermix (Bio-Rad) and 5 μl of cDNA diluted 25 times. The real-time RT-PCR program included an initial 3-min denaturation step at 95°C and 40 cycles of denaturation at 95°C for 15s and annealing for 30s at the optimal temperature indicated in Additional file 4. For each gene, real-time RT-PCR efficiency was estimated from a serial dilution of cDNA (5 to 500 times) and taken into account in the data analysis performed with the ΔΔCT method . Two housekeeping genes encoding ribosomal protein L8 (RPL8, GenBank accession number: DQ440262) and S7 (RPS7, GenBank accession number: AY380336) were used for normalization. Results were represented as mean expression ratios between Bti-resistant and susceptible larvae (± SE).