The kdr-bearing haplotype and susceptibility to Plasmodium falciparum in Anopheles gambiae: genetic correlation and functional testing

Mosquito sampling and wild Plasmodium falciparum infection

Mosquitoes were sampled as larvae using the standard dipping method as previously described [28]. In order to collect a random sample of essentially unrelated individuals, fewer than ten larvae were collected from any one larval site and a single collection comprised larvae from >50 larval sites. Collections were done in the Sudan Savanna region of Burkina Faso in the village of Goundry (12°30′N, 1°20′W), 30 km N of the capital city, Ouagadougou, across months of the malaria transmission season of 2007 and 2008 [29]. Larvae were raised in Centre National de Recherche et de Formation sur le Paludisme (CNRFP) insectaries in Ouagadougou to adulthood under standard insectary conditions [29]. Following emergence, 3-day old adults were challenged with wild P. falciparum by experimental infection. Feeding was done on an artificial membrane in a water-jacketed feeding device as described previously using gametocytaemic blood obtained from study participants [29]. Unfed mosquitoes were excluded from analysis and infection levels for fed mosquitoes were determined by counting midgut oocysts 7–8 days post infection. Genomic DNA was extracted from carcasses for SNP genotyping. Population sub-group composition was determined by published molecular diagnostic assays [28, 30]. Mosquitoes used for genetic association analysis came from five successful experimental infections on wild gametocytes using criteria defined previously [6, 28, 29, 31]: oocyst infection prevalence ≥30 % and oocyst intensity in at least one mosquito in the group of ≥10 oocysts. These criteria were applied to the mosquitoes raised together as a cohort and fed on the same blood, and identified the experiments with broad phenotype distributions. This infection-quality filter used for infection with wild gametocytes, which are inherently variable and noisy, assures that all analysed individuals came from an infection session with the power to distinguish between levels of susceptibility, free from confounding noise due to uninfective blood, technical failure or related factors. Low prevalence infections in the mosquitoes fed on gametocytaemic blood can potentially result from multiple causes, but regardless the cause, such groups are not informative for genetic analysis because of poor phenotype distribution. Mosquito genetics can be excluded as a cause, because the fed cohort is comprised of random unrelated individuals.

SNP genotyping and genetic association

DNA samples from dissected mosquito carcasses were subjected to whole genome amplification (WGA) (Genomiphi, GE Health Sciences) using supplied protocols. Following WGA, DNA was ethanol precipitated, concentrations determined by the picogreen method [32], and 500 ng was used for SNP genotyping by the Illumina GoldenGate method [33] by the Boston Children’s Hospital Molecular Genetics Core Facility (IDDRC). The SNPs and their primer sequences are available in Additional file 1: Table S1. Data were analysed using BeadStudio software [33]. Genotypes of mosquitoes for the Pfin6/para haplotype were independently verified using the Sequenom MassArray method, with a set of seven diagnostic SNPs at the following genomic positions in the PEST reference genome: 2L.1272741 [T/C], 2L.1834476 [A/G], 2L.1970368 [C/T], 2L.2081228 [T/C], 2L.2430786 [C/T], 2L.2489023 [G/A], 2L.2489212 [T/C]. As these seven SNPs are in perfect linkage disequilibrium, they are all equally informative and thus comprise a sevenfold redundant assay for a single genotype feature, that is, the haplotype they are carried upon. Because the seven SNPs are all completely cross-informative for the haplotype, additional markers in this region would not increase resolution of the feature. For this, aliquots of the same DNA samples used for Illumina genotyping were genotyped using Sequenom at the University of Minnesota Genomic Center (UMGC). Sequenom genotypes were analysed for genetic association with malaria infection phenotype using the Haploview program [34].

Gene silencing and functional assays

A colony of An. coluzzii mosquitoes (Fd03) was initiated with the eggs of ten mated, gravid females captured in village houses in Goundry, Burkina Faso. The colony is wild-type for para and is maintained using standard laboratory conditions at the Institut Pasteur. Double-stranded RNAs (dsRNA) were synthesized from PCR amplicons using the T7 Megascript Kit (Ambion) as described previously [35]. The sequences of the primers used for amplification of ClipC9 dsRNA were:

T7_CLIPC9F:TAATACGACTCACTATAGGACCGTGCTGCAGAATGACTGC T7_CLIPC9R: TAATACGACTCACTATAGGTATCCCTAGTAGCACTAACCG.

The primers used for amplification of the para/VSC dsRNA were:

T7-4707F:TAATACGACTCACTATAGGAGGGCTATCCGGGAAATTGTGG.

T7-4707 R:TAATACGACTCACTATAGGTGAAGCGTCTGTTCCGCCTCC.

The underlined portion of the primer sequence is the T7 tag. For each targeted gene, 200 ng of dsRNA (but not more than 200 nl, depending on the concentration) were injected into the thorax of cold-anesthetized 1-day old An. coluzzii females using a microinjector (Nanoject II, Drummond). Validation of gene silencing efficiency was done by RT-PCR. In all cases, treatment with dsRNA led to decreases in target transcript abundance compared to dsGFP-injected controls (Additional file 2: Figure S1). The same numbers of mosquitoes were treated with control and experimental dsRNA in each replicate.

For gene silencing assays, mosquitoes were challenged with P. falciparum strain NF54 gametocytes, cultured using an automated tipper-table system implemented in the CEPIA mosquito infection facility of the Institut Pasteur, as previously described [35]. For experimental infection, 10 ml of medium containing mature gametocytes was centrifuged at 2000 rpm, and the cell pellet was resuspended in an equal volume of normal type AB human serum. The infected erythrocytes were added to fresh erythrocytes in AB human serum and transferred to a membrane feeder warmed to 37 °C. Mosquitoes were allowed to feed for 15 min. Unfed mosquitoes were removed, so that only fully engorged females were used for further analysis. Mosquitoes were maintained on 10 % sucrose solution supplemented with 0.05 % para-amino benzoic acid. At 7–8 days after the malaria infection challenge, all surviving mosquitoes were dissected. Thus, the dissected number, given as sample size (n) values below infection histograms, indicates survivorship, and control versus treatment values are compared as an estimate of gene-specific or overall mortality effects. At least two independent replicate infections were performed per condition. Replicates were analysed independently using the tests described below. If at least one replicate met the significance criterion of p ≤ 0.05, a third replicate was done.

Analysis of infection phenotypes

For infections on either wild or cultured gametocytes, infection phenotype was determined at 7–8 days post blood meal. Mosquito midguts were dissected, stained with dibromofluorescein (Sigma Chemical) and the number of midgut oocysts determined by light microscopy. For mosquitoes fed on wild gametocytes, infection prevalence was measured (the fraction of mosquitoes with at least one oocyst), while for functional genomic tests using mosquitoes fed on cultured gametocytes, two infection phenotypes were measured, infection prevalence, and oocyst intensity, counts of oocysts in only those individuals with ≥1 oocyst. Oocyst intensity comparisons were not used for wild gametocyte infections because the sample size of mosquitoes with ≥1 oocyst did not provide sufficient statistical power for intensity.

For statistical analysis, comparisons of infection prevalence were made using the Chi Square test, and comparisons of oocyst intensity using the non-parametric Wilcoxon Mann–Whitney test. Replicates were analysed independently using the tests described above and p values from independent tests of significance were combined using the meta-analytical approach of Fisher [36], and this combined p value is reported here. The threshold for significance was defined as p = 0.01.

Ethical considerations

For collection of blood from P. falciparum gametocyte carriers for experimental membrane feeder infection of mosquitoes, the study protocol was reviewed and approved by the institutional and national health ethical review board (Commission Nationale d’Ethique en Santé) of Burkina Faso (code No 2006-032). The study procedures, benefits and risks were explained to parents or legal guardians of children and their informed consent was obtained. Children of parents or guardians who had given consent were brought to CNRFP the day of the experiment for gametocyte carrier screening. All children were followed and symptomatic subjects were treated with the combination of artemether–lumefantrine (Coartem^®) according to relevant regulations of the Burkina Faso Ministry of Health.

Structure of the para gene locus

Wild mosquitoes were collected as larvae in Burkina Faso, were reared to adults in the insectary, and were challenged with wild P. falciparum parasites by membrane feeding on blood from gametocyte carriers, yielding 363 mosquitoes belonging to the An. gambiae species complex that were used in the study.

The mosquito DNAs were genotyped with 105 SNPs spanning a broad ~15 Mb region containing the para gene across the centromere of chromosome 2, between coordinates chromosome 2R:55 Mb to chromosome 2L:10 Mb (Additional file 1: Table S1; Fig. 1, bottom track of 105 SNPs shown as black vertical lines). Analysis of the individual mosquito genotypes reveals a region of elevated local homozygosity that defines the boundaries of two allelic haplotypes, termed HapA and HapB (Fig. 1). When mosquitoes were categorized by their HapA/HapB genotype, homozygotes are evident as regions of high linkage, with additional unlinked variation seen as dips in homozygosity (Fig. 1, homozygosity plot for each genotype). The core haplotype across both alleles is ~3 Mb in size (Fig. 1, blue horizontal lines superimposed on red homozygosity plot). This corresponds to the location and estimated 3.3 Mb size of the haplotype bearing the kdr allele of the para gene that was reported as an introgression from An. gambiae to An. coluzzii [24, 25].

The current sample set was comprised of three population groups of the An. gambiae species complex (An. gambiae, n = 46, An. coluzzii, n = 53, and Goundry form, n = 264). The Goundry form is a population sub-group within An. gambiae sensu lato with unknown taxonomic status or epidemiology for malaria transmission [28, 37, 38]. The haplotype structure of the para gene region in Goundry was not previously known. Of the 105 genotyped markers in this region, a sub-set of 26 SNPs span the haplotype (Fig. 1, middle track of SNPs shown as purple vertical lines) and were used to examine segregation of haplotype alleles across the three taxa. Consistent with observations that the kdr mutation introgressed from An. gambiae to An. coluzzii upon a haplotype [20, 21, 24, 25], the 26 SNPs in this genomic region did not distinguish between An. gambiae, An. coluzzii and the Goundry form (Fig. 2a). This result indicates that the para-bearing haplotype is identical by descent (IBD) across the three taxa, with strong ancestral sharing as well as the acquisition of additional mutations over time.

Genetic association analysis for Plasmodium falciparum susceptibility

Tests for genetic association are done within panmictic population groups to avoid stratification artifacts [39]. Of the three population groups in the current study, An. gambiae have no power to test for association because the kdr allele is essentially fixed in An. gambiae, consistent with other reports [20, 21], and despite segregation of the kdr allele within An. coluzzii, the sample size was insufficient for an adequately powered statistical test. Thus, in the current study, only the Goundry group sample was appropriate for an association test of the para-bearing haplotype. As shown (Fig. 2A), a haplotype with a common evolutionary origin segregates in all three population groups. Thus, although the genetic background may differ among taxa, the same para-bearing haplotype is present in all of these taxa, and tests for association of the haplotype with P. falciparum infection levels were performed within the Goundry group.

A sub-set of seven completely correlated SNPs were chosen as haplotype-tagging markers for genetic association (Fig. 1, top track of SNPs shown as blue vertical lines). The seven SNPs mark ancestral differentiation, are in perfect linkage disequilibrium (LD) (thus are redundant tests for the same feature), and efficiently tag HapA and HapB alleles. Association of the para-bearing haplotype with susceptibility to P. falciparum was tested in the Goundry sub-group mosquitoes using the haplotype-tagging SNPs combined as a single redundant assay for the haplotype. Mosquitoes homozygous for haplotype allele HapA were more susceptible to infection with P. falciparum than heterozygotes or HapB homozygotes (p = 0.0293 for oocyst infection prevalence; Fig. 2B). The locus was named Plasmodium falciparum infection locus 6 (Pfin6) following naming convention [31].

Functional dissection of the Pfin6 locus

The haplotype that carries the mapped Pfin6 locus as well as the para gene is ~3.3 Mb in length [24, 25] and contains 65 predicted genes (Fig. 2C; Additional file 3: Table S2). Because the haplotype contains long blocks of SNPs in perfect linkage, resolution by recombinational fine-mapping is not expected to be productive. Fine-mapping would be further complicated by the centromeric location of the mapped locus, as centromeres typically have lower than average recombination rates. Instead, in the present study a candidate gene approach is employed, dissecting the locus by functional testing of candidates ascertained by predicted gene function and other evidence. It is important to note that candidate gene assays can provide suggestive evidence, but do not constitute proof of genetic causation.

Of the genes in the locus (AGAP004677-AGAP004742), two are potentially involved in host-defense processes: AGAP004719 (2L:2.71 Mb), with strong prediction of immune function, identified as ClipC9, a clip-domain serine protease homologous with the Drosophila extracellular serine protease Persephone (CG6367), which mediates activation of the Toll signalling pathway in response to fungal infection [40]; and AGAP004699 (2L:1.97 Mb), weakly predicted as potentially immune, which codes for a serine/threonine-protein kinase with a Drosophila orthologue involved in haemocyte differentiation (Raf homologue serine/threonine-protein kinase phl, CG2845, [41]. The locus also contains AGAP004707 (2L:2.35 Mb), the para gene in which the kdr mutant is associated with response to pyrethroid insecticides.

The following lines of evidence prioritized the functional testing of candidates ClipC9/AGAP004719 and para/AGAP004707 for their effect on the outcome of P. falciparum infection, using RNAi-mediated gene silencing assays. ClipC9 is predicted to be catalytically active as a serine protease, bearing all cysteine residues involved in disulfide bonds, as well as the conserved residues of the catalytic triad. Expression of ClipC9 was upregulated following infection with Gram-negative bacteria and down-regulated following infection with Gram-positive bacteria [42, 43], but was unchanged following infection with Plasmodium [43]. The para gene would not be selected by the same criteria as a candidate for functional testing, because it does not display any known or predicted immune function, but it was chosen on an ad hoc basis because of reports of association of kdr with infection [26, 27]. The para gene influence upon any phenotype has not been functionally analysed by gene silencing, according to published reports.

Mosquitoes from an An. coluzzii colony recently initiated from the Burkina Faso study site were injected with double-stranded RNA (dsRNA) targeting ClipC9, para, or control GFP. Four days after dsRNA treatment, decreased transcript levels were verified (Additional file 2: Figure S1), mosquitoes were fed on cultured P. falciparum gametocytes, and 7–8 days later midgut oocysts were counted. Infection intensity (numbers of oocysts in mosquitoes with ≥1 oocyst) was consistently higher in the dsClipC9-treated individuals as compared to controls (Fig. 3; p = 6.23 × 10⁻⁴; p values from three replicate experiments combined using the Fisher Method; individual p values 0.005, 0.089, 0.017). There was no effect of dsClipC9 treatment upon infection prevalence (the proportion of mosquitoes with ≥1 oocyst, p = 0.78; individual p values 0.658, 0.678, 0.450).

For the para voltage-gated sodium channel (VSC) gene, there was no statistically significant difference for either infection prevalence or intensity between dsVSC-treated as compared to control mosquitoes (infection prevalence p = 1.0, infection intensity p = 0.10, Fig. 4). Thus, under these conditions, abolishing para gene transcriptional activity does not appear to influence parasite development, although the effect of transcript silencing on protein levels is not known.

Testing gene function by silencing, as shown here, is distinct from testing allele function by genetic association, as shown above. Gene silencing mediated by a ~500 bp double-stranded RNA targeting the gene is immune to sequence differences greater than the kdr SNP variant between para alleles [44]. The colony used for the functional test was fixed for wild-type para, and transcript silencing is efficient (Additional file 2: Figure S1). Thus, para gene transcripts are abolished by silencing, regardless of the genetic variants they carry.

The effect of gene silencing on mosquito mortality was monitored by comparison of control dsGFP and experimental dsRNA treatments. The same numbers of control and experimental mosquitoes were injected with dsRNA per replicate, and all surviving mosquitoes were dissected, thus the sample sizes (n values shown on Figs. 3, 4) indicate overall survival per treatment. There was no evidence of a consistent difference in survival between mosquitoes treated with dsClipC9 as compared to dsGFP controls (Fig. 3), or between dsVSC-treated as compared to dsGFP controls (Fig. 4). Thus, it appears that neither ClipC9 nor para gene activity is biologically essential for mosquito viability, although minor effects on survival cannot be ruled out.

The Pfin6 malaria infection locus and kdr are carried on the same haplotype

Results presented here show that a ~3 Mb haplotype is shared by introgression as an identical by descent (IBD) feature among three otherwise reproductively isolated groups of the An. gambiae species complex. Two of the three groups, An. gambiae and An. coluzzii, display overall equivalent susceptibility to wild P. falciparum [29, 45–48], while the Goundry form is more susceptible to experimental infection with wild parasites [28]. However, because the haplotype is a shared feature with common evolutionary origin among the three population groups, it is unlikely to explain this phenotypic difference between the groups.

The unit of introgression of the kdr allele was not simply the kdr variant nucleotide, but rather an extended haplotype bearing para and >60 other neighbouring genes. The ~3 Mb size of the introgressed haplotype reported here is consistent with other reports [24, 25]. It has been suggested that exposure to insecticides used in agriculture or on bed nets could be driving the spread of kdr, based on the temporal correlation of the two events [23, 49, 50]. The finding that the pathogen-resistance locus, Pfin6, maps to the same haplotype suggests a potential complementary explanation for the adaptive introgression of this haplotype. Components of the immune systems of An. coluzzii and An. gambiae mosquitoes are under distinct selective pressure that has strongly altered the alleles and diversity at particular immune loci unrelated to Pfin6 [37, 51, 52]. In those cases, it is hypothesized that exposure to distinct pathogen profiles in the different ecological niches of An. coluzzii and An. gambiae mosquitoes could be the source of this strong recent selection. The presence of a functional immunity locus upon this haplotype might explain at least part of the introgression pressure based on differential immunity to environmental microbes.

Functional genomic analysis of the haplotype carrying Pfin6 and para

Based on prediction of potential immune function one high-priority candidate gene was identified within the Pfin6 locus, the clip domain serine protease, ClipC9. This does not rule out that a different or multiple genes could actually underlie the Plasmodium susceptibility phenotype of Pfin6 detected in nature. ClipC9 was functionally confirmed, by gene silencing, as a gene that influences mosquito infection levels with P. falciparum. Functional assays of candidate genes can demonstrate that a gene influences a similar phenotype to a genetically mapped locus, and thus could be a plausible genetic candidate. However, functional assays are not proof of genetic causation. Subsequent work will be required to establish whether naturally-occurring genetic variation in ClipC9 underlies the genetic effect of the Pfin6 locus in wild mosquitoes.

Pfin6 linkage with kdr and malaria transmission

The observed linkage reported here between P. falciparum susceptibility and insecticide resistance, and observations in other recent studies using different methodological approaches show the same direction of association [26, 27]. The functional analysis of para suggests that para itself is not responsible for the malaria susceptibility phenotype, although an effect in nature cannot be ruled out. The influence of transcript loss on protein abundance could not be measured because of the lack of existing antibodies against this highly hydrophobic protein. Protein abundance is often safely ignored in functional genomic experiments that produce a detectable phenotype, although in the case of para the absence of a transcript depletion phenotype for either infection or mortality means that the unknown protein level remains a caveat.

In principle, wild kdr-bearing mosquitoes could contribute disproportionately to malaria transmission if they live longer than wild-type mosquitoes, although the comparative lifespan of kdr mosquitoes in field populations is not known. Here, mosquitoes were not exposed to insecticide treatment, and all individuals were offered a single blood meal and dissected 7–8 days afterwards, so any kdr-based differences in longevity, if they exist, do not influence the results.

Loss of para transcript by gene silencing did not result in higher mosquito mortality. This viability result is consistent with observations from the para orthologue in Drosophila (FlyBase ID FBgn0264255), which harbours many described mutations, some lethal but many viable [53], often displaying temperature-sensitive behavioural defects such as excitability that have been used as models for neurological disorders [54, 55]. Based on the finding that mosquitoes depleted for para transcript are viable, it would be interesting to assay para-silenced mosquitoes for influence on other behavioral phenotypes that could be related to insecticide response.