An A. funestus BAC library from The Institute for Genomic Research, Notre Dame University, was screened by PCR using primers from nine P450 genes located within the boundaries of the rp2 QTL using the synteny projection with the A. gambiae chromosomal map: these genes are CYP6Z3, CYP6Z1, CYP6Y1, CYP6M8, CYP6S1, CYP6M1, CYP4H18, CYP6N2 and CYP6M7. The primers used are listed in Supplementary Table S1. DNA of whole 384-well plates was pooled and a PCR carried out for each plate. The positive plates were then subdivided into six column pools and 4 row pools and the PCR screen repeated. Finally, individual colonies from the set of 16 identified from the pooled column and row screen were used as template to identify the individual clone containing the markers of interest. The BAC clone was grown at 37°C overnight and harvested in a glycerol solution and stored at −80°C. The size of the BAC clone was estimated after a restriction digestion using the BamH1 restriction enzyme to separate the insert from the vector. The BAC clone was then fragmented by sonication into shorter, random sequences of around 2–5 kb. These small fragments were cloned into the plasmid vector (pC31) and sequenced. The reads were trimmed using Seqman (http://www.dnastar.com/web/r13.php) and assembled using Phrap (P. Green unpublished; http://www.genome.washington.edu/UWGC/analysistools/Phrap.cfm).

The programme Genemark.hmm version 2.2 was used to locate genes in the assembled sequence of the BAC clone. Putative genes were annotated by using BlastX (http://www.ncbi.nlm.nih.gov/blast) and the predicted transcripts were also compared with A. gambiae transcripts in Vectorbase (www.vectorbase.org). Further detailed annotation of the P450 genes was aided by the P450 site (http://p450.sophia.inra.fr/). Sequence alignments of A. funestus and A. gambiae genes were carried out using ClustalW (^{Thompson et al., 1994}). Exact boundaries of the exon/intron were manually inspected and compared with that of A. gambiae using the DNAstar sequence analysis package. MEGA 4.0 (^{Tamura et al., 2007}) was used to construct a Neighbour-Joining tree of the A. funestus genes in comparison with A. gambiae.

Two AIL families at F₆ generations were used for this study. These families (Family 1 and Family 10) were generated from reciprocal crosses between the susceptible and the resistant strains as previously described (^{Wondji et al., 2009}). Briefly, a two exposure times approach of 30 min and 2 h was used to select the most susceptible and the most resistant mosquitoes, respectively, as described previously in order to minimise the level of phenotype misclassification and increases the power of QTL detection (^{Lander and Botstein, 1989}).

In order to carry out a fine-scale mapping of the rp2 QTL, SNPs were identified in all the genes detected in the BAC clone and in other genes spanning the rp2 QTL boundaries in the 2L chromosome such as the glutathione-s-transferase GSTe2, the cuticular protein genes CPLC5, CPLC9 and CPLC8, the P450s CYP9M1, CYP303A1 and CYP4H18, and AGAP007980 and AGAP009073. The cuticular protein genes were amplified using primers listed in Supplementary Table S2. Genomic DNA was extracted using the LIVAK method (^{Livak, 1984}) and amplified using the primers listed in Supplementary Table S3. For each family, the F₀ parental female was amplified with 3F₁ in order to detect the informative SNPs to be used for genetic mapping. The PCR were carried out using 10 pmol of each primers and 30 ng of genomic DNA as template in 25 ml reactions. PCR products were purified using the QIaquick PCR purification kit (Qiagen, Hilden, Germany) and directly sequenced on both strands. Sequences for each gene were analysed to detect the polymorphic sites manually using BioEdit and as sequence differences in multiple alignments using ClustalW (^{Thompson et al., 1994}).

SNPs were identified as transitions or transversions in coding and non-coding regions. Genetic diversity analysis was performed using DnaSP 5.1 (^{Rozas et al., 2003}). The nucleotide diversity π, was calculated for each gene as well as the haplotype diversity. The average number of synonymous substitutions per synonymous site (Ks) and non-synonymous substitutions per non-synonymous site (Ka) was computed.

The expression pattern of 28 genes located in the rp2 BAC clone or in 2l chromosome was compared between the resistant strain FUMOZ-R and the susceptible strain FANG using the GeXP multiplex gene expression profiling method from Beckman Coulter as previously described (^{Wondji et al., 2009}). Total RNA was extracted from three batches of ten (1–3-day old) female and male adult mosquitoes using a PicoPure RNA isolation kit (Arcturus, Mountain View, CA, USA) according to manufacturer's instructions. Total RNA quantity and quality were assessed using Nanodrop spectrophotometer (Nanodrop Technologies, Oxfordshire, UK) and Bioanalyzer (Agilent, Santa Clara, CA, USA), respectively. The GenomeLab eXpress Profiler programme was used for automated primer design, calculation of relative gene expression values, data checking and data analysis. The quantitative PCR reaction was carried out using the GenomeLab GeXP Start Kit (Beckman and Coulter, Brea, CA, USA) according to the protocol provided and using the primers listed in Supplementary Table S5. The amplified qPCR products were diluted 10 times and added in a 96-well microplate with the DNA size standard-400 and ran on the GenomeLab GeXP genetic analysis system. The expression level of the RSP7 ribosomal gene was used to normalise for variation in total cDNA concentration. A two-sample t-test was used to compare the results between the two strains.

Genes with a significant overexpression in the FUMOZ-R strain from the GeXP multiplex expression profiling was further analysis using individual quantitative reverse transcription PCR (qRT–PCR). To increase the specificity of the amplification, reverse primers were designed in the 3′UTR region for all the genes. Primers are listed in Supplementary Table S6. One microgram of total RNA from each of the three biological replicates for FUMOZ-R and FANG was used as template for cDNA synthesis using the Superscript III (Invitrogen, Carlsbad, CA, USA) with oligo-dT20 and RNase H, according to the manufacturer's instructions. A serial dilution of cDNA was used to establish standard curves for each gene in order to assess PCR efficiency and quantitative differences between samples. The qRT–PCR amplification was carried out in a MX 3005 real-time PCR system (Agilent) using Brilliant III Ultra-Fast SYBR Green QPCR Master Mix (Agilent). Ten nanogram of cDNA from each sample was used as template in a three-steps programme involving a denaturation at 95°C for 3 min followed by 40 cycles of 10 s at 95°C and 10 s at 60°C and a last step of 1 min at 95°C, 30 s at 55°C and 95°C at 30 s. The relative expression and fold change of each target gene in FUMOZ-R relative to FANG was calculated according to the 2^−ΔΔCT method incorporating PCR efficiency after normalisation with the housekeeping RSP7 ribosomal protein S7 (AGAP010592) and the actin 5C (AGAP000651) genes. A two-sample t-test was used to compare the results between the two strains.

In order to detect any potential polymorphism associated with resistance, the full length of the candidate genes detected after qRT–PCR was amplified using the cDNA synthesised from total RNA as template between the resistant and the susceptible strains. Amplification was performed with the Phusion polymerase with the following conditions: 1 cycle at 95 °C for 5 min; 35 cycles of 94 °C for 20 s, 57 °C for 30 s and elongation at 72 °C for 60 s; followed by 1 cycle at 72 °C for 5 min. The PCR products were purified using the QIaquick PCR purification kit (Qiagen) and cloned into the pJET1.2/blunt cloning vector using the CloneJET PCR cloning kit (Fermentas, Burlington, ON, Canada). Positive clones were purified and sequenced on both strands. Sequences were analysed as described above.

Nine positive 384-well plates containing one or more of the loci spanning rp2 were identified from the A. funestus BAC library. A single individual BAC clone presented a positive PCR result for eight of the nine genes tested (CYP4H18 PCR was negative). The estimation of the size of this clone indicated that it was above 100 kb and likely to span the rp2 QTL.

A complete genomic organisation of this P450 cluster is presented in Figure 1. The genomic arrangement of the fifteen P450 genes is exactly the same as in A. gambiae. The size of the intergenic spaces is also similar between the two species with the lowest intergenic space in this cluster being 78 bp between AfCYP6S1 and AfCYP6R1 (This is also the lowest in A. gambiae with 71bp), while the longest intergenic space is between AfCYP6M4 and the salivary protein gSG7-2 (24 739 bp), which is similar in A. gambiae.

All the 15 P450 genes possess two exons separated by short introns from 58 bp for AfCYP6S1 to 87 bp for AfCYP6M7. These intron sizes are also similar to that of A. gambiae (Supplementary Table S8). The amino acid sequences of the fifteen P450 genes varied from 492 amino acids for AfCYP6Z1 and AfCYP6Z3 to 514 for AfCYP6Y2, similar to A. gambiae with the only difference that AfCYP6Y2 is 520 amino acids long in this species. The highly conserved motifs C-helix WxxxR, ExLR and PERF are all present as indicated in the multiple alignment. Similarly, the haem-binding loop (the P450 signature sequence) PFxxGxRxCxG/A (^{Feyereisen, 2011}) and the membrane targeting/anchor hydrophobic N-terminal region are present in all fifteen P450s (Supplementary Figure S1). The GxE/DTT/S is also present in all genes apart from AfCYP6R1 as observed in A. gambiae. The orthology of the 15 P450 genes to the respective gambiae genes is confirmed in the Neighbour-joining tree in Supplementary Figure S2 showing that the A. funestus P450 genes are primarily closer to their orthologues in A. gambiae. This is also the case for the three AfCYP6M1 genes all closer to AgCYP6M1 in A. gambiae with similarity rates of 79%, 82% and 83%, respectively, for AfCYP6M1b, AfCYP6M1a and AfCYP6M1c.

The triplication of AfCYP6M1 was confirmed by PCR amplification and sequencing of the 5′flanking region of each copy of the three genes using genomic DNA from laboratory and field samples. The deduced amino acid sequences of the three genes are all of the same length with 498 amino acids each. There are 59 amino acid substitutions between AfCYP6M1a and AfCYP6M1b among which 13 replacements between amino acid of opposite physicochemical properties resulting in an 88% similarity (Supplementary Figure S3). AfCYP6M1a and AfCYP6M1c have 45 variant amino acids (90% similarity) among which 8 replacements between amino acid of opposite physicochemical properties (Supplementary Figure S3). AfCYP6M1b and AfCYP6M1c have 57 variant amino acids (88% similarity) among which 12 replacements between amino acid of opposite physicochemical properties (Supplementary Figure S3).

Four CPLCG genes were successfully amplified and sequenced in A. funestus samples. As some of these genes have been previously associated with insecticide resistance in A. gambiae (^{Vontas et al., 2005}; ^{Djouaka et al., 2008}), the aim was to assess the correlation between these genes and rp2 QTL. Detailed information about these genes is presented in Supplementary Table S9. Overall, high similarity is observed between the four CPLCG genes and their orthologues in A. gambiae (Supplementary Table S9) with percentage of protein sequence similarity between 67 and 86%. The CPLCG3 protein in A. funestus, previously associated with possible reduced penetration resistance mechanism in A. gambiae (^{Vontas et al., 2005}; ^{Djouaka et al., 2008}), had seven amino acid differences with the A. gambiae orthologue.

In total, 565 SNPs from 28 genes screened were identified after analysing 20 431 bp of sequencing traces indicating that one SNP is present every 36 bp. The patterns of polymorphism observed in this study is similar to the patterns previously described in this species from 50 genes (^{Wondji et al., 2007a}). This is the case for the ratio of transition/tranversion substitutions with transition SNPs significantly predominant in the total (60.9%) over transversion substitutions (39.1%) compared with a ratio of 62% vs 38% in the previous study (^{Wondji et al., 2007a}). This bias towards transition substitutions was also more pronounced in the coding region with 62.7% (271 SNPs out of 432) of transition and 37.3% (161 SNPs) of transversion slightly lower (but not significantly) than the 66.3% vs 33.7% seen previously (^{Wondji et al., 2007a}). More transversion SNPs were significantly observed in non-coding regions (47.5%) than in coding region (37.3%) as also observed in other species such as A. gambiae (^{Morlais et al., 2004}) and Aedes aegypti (^{Morlais and Severson, 2003}), confirming that SNPs occur more frequently as transition in coding regions than in non-coding regions. This difference between coding and non-coding regions is also reflected in the frequency of SNPs, which is one SNP every 39 bp in coding regions vs one SNP every 24 bp in non-coding regions (5′UTR, 3′UTR and intron combined).

A higher frequency of SNPs was observed at the third codon position (71.9%) than at the first or second position (Table 1). More than 2/3 (308) of the 432 coding SNPs were synonymous, whereas <1/3 (124) were non-synonymous or replacement SNPs. Indels were also detected (1–6 bp long) or inferred owing to overlapping peaks in the sequencing traces. Twenty-five triallelic and six tetrallelic SNPs were detected confirming that multiallelic SNPs are present in this mosquito species although at low frequency.

The average nucleotide diversity per gene (π=0.0105) was higher than the estimate of 0.0072 in the previous study (^{Wondji et al., 2007a}). The average nucleotide diversity in non-coding DNA (0.013) was lower than in synonymous sites of the coding regions (0.028). This is an indication that non-coding regions are under greater purifying selection than synonymous sites within coding regions. The Ka/Ks ratio was 0.11, which is lower than the estimate of 0.18 seen in the previous study.

In order to detect the candidate genes associated with the rp2 QTL, the expression profile of the 15 P450 genes from the rp2 BAC clone was assessed in comparison with the reference ribosomal protein S7 gene (RSP7) between the resistant strain FUMOZ-R and the susceptible FANG strain. This multiplex qPCR experiment detected no expression signal for AfCYP6M1b, AfCYP6S1 and AfCYP6S2 either in resistant or in susceptible samples. The gene with the highest expression level was AfCYP6M8 orthologue of AgCYP6M2 gene regularly associated with pyrethroid resistance in A. gambiae (^{Djouaka et al., 2008}; ^{Muller et al., 2008}; ^{Stevenson et al., 2011}). However, this high expression of AfCYP6M8 is observed in both resistant and susceptible samples indicating that contrary to its A. gambiae AgCYP6M2 orthologue, AfCYP6M8 is not associated with pyrethroid resistance in this FUMOZ-R strain. A significant differential expression (P<0.01; t-test) was observed at two genes, AfCYP6M7 and AfCYP6Z1, between the resistant FUMOZ-R and the susceptible FANG, with, respectively, a 2.4- and 2.82-fold-change overexpression in the resistant strain (Figure 4a). No significant difference was observed for the remaining genes, although a 1.7-fold change was also observed for AfCYP6Z3 but at P=0.065 (t-test).

The analysis of six other genes of interest on the 2L chromosome away from the rp2 BAC clone indicated a significant downexpression of the cuticular protein CPLCG3 (CPLC8) gene in the resistant strain contrary to previous observation made in A. gambiae associating this gene with pyrethroid resistance in ^{Djouaka et al., 2008}. A significant overexpression of the GSTe2 gene (P=0.017) was also observed in the resistant FUMOZ-R strain. No informative SNP was identified in this gene and therefore was not included in the QTL mapping to assess its correlation with resistance phenotype. No differential expression was observed for the cytochrome c oxidase subunit I gene (orthologue of the A. gambiae AGAP009537 gene), which exhibited a significant genotype/phenotype association.

The expression profile of AfCYP6Z1, AfCYP6M7 and GSTe2 upregulated in FUMOZ-R through the multiplex GeXP expression profiling was further assessed individually using the qRT–PCR. This was also done for AfCYP6Z3, which showed a 1.7-fold change although not significant. All genes tested had PCR efficiency between 90 and 102%. A significant overexpression was observed for all the four genes when their (2^−ΔΔCt) relative expression was compared with that of the two housekeeping genes RSP7 and Actin either individually or when combined (Figure 4b). The highest fold change is observed for CYP6Z1 with an 11.16-fold change upregulation in FUMOZ-R compared with the susceptible FANG. AfCYP6Z3 is 6.3-fold overexpressed in FUMOZ-R, while AfCYP6M7 is also upregulated in this resistant strain at 5.5-fold. GSTe2 also showed a significant upregulation in FUMOZ-R but at a lower fold change of 2.82.

The polymorphism of the AfCYP6Z1 and AfCYP6Z3 was analysed in order to identify potential mutations associated with pyrethroid resistance. The full length of the coding region of AfCYP6Z1 was sequenced for eight clones (four FANG and four FUMOZ-R), whereas 11 clones (five FANG and six FUMOZ-R) were sequenced for AfCYP6Z3. The two genes exhibited a high polymorphism with a total of 49 and 59 polymorphic sites observed, respectively, for AfCYP6Z1 and AfCYP6Z3, while 19 amino acid changes were observed for AfCYP6Z1 (Supplementary Figure S5) and 9 for AfCYP6Z3 (Supplementary Figure S6). This level of polymorphism is similar between the FANG and the FUMOZ-R strain with, nevertheless, a slightly higher number of mutations in the susceptible strain for the two genes. Other parameters of this polymorphism are summarised in Table 3. A Neighbour-joining tree of the haplotypes of both genes indicated a cluster of haplotypes specific to each strain for both genes (Figure 5) but there are also haplotypes similar between the two strains.