Evidence for population-specific positive selection on immune genes of Anopheles gambiae

Abstract

Host-pathogen interactions can be powerful drivers of adaptive evolution, shaping the patterns of molecular variation at the genes involved. In this study, we sequenced alleles from 28 immune-related loci in wild samples of multiple genetic subpopulations of the African malaria mosquito Anopheles gambiae, obtaining unprecedented sample sizes and providing the first opportunity to contrast patterns of molecular evolution at immune-related loci in the recently discovered GOUNDRY population to those of the indoor-resting M and S molecular forms. In contrast to previous studies that focused on immune genes identified in laboratory studies, we centered our analysis on genes that fall within a quantitative trait locus associated with resistance to Plasmodium falciparum in natural populations of A. gambiae. Analyses of haplotypic and genetic diversity at these 28 loci revealed striking differences among populations in levels of genetic diversity and allele frequencies in coding sequence. Putative signals of positive selection were identified at 11 loci, but only one was shared by two subgroups of A. gambiae. Striking patterns of linkage disequilibrium were observed at several loci. We discuss these results with respect to ecological differences among these strata as well as potential implications for disease transmission.

Keywords: Anopheles gambiae; Plasmodium falciparum; host−pathogen interaction; immune gene; positive selection.

Figure 1

Barplot distributions of genetic differentiation among groups of A. gambiae at immune genes. Genetic differentiation was estimated using Weir and Cockerham’s unbiased estimator of F_ST between groups of A. gambiae at each gene separately. Loci are arranged according to their genomic coordinates with the 2La inversion presented in the inverted arrangement. The vertical bar colors specify genomic region according to the legend. The schematic indicates the physical distribution of the genes on 2L, with C and T representing the centromere and telomere, respectively. The PRI region and 2La are also indicated on the chromosome schematic. Asterisks indicate F_ST values significantly greater than zero as determined by permutation tests (see Methods). Both the M and S molecular forms were compared to GOUNDRY separately and exhibited qualitatively similar levels and patterns of differentiation, so only the S-form comparison is presented here.

Figure 2

Distribution of Tajima’s D. Tajima’s D was calculated for all genes using synonymous sites and both boxplots and data points are presented. The dotted line indicates the expected value of D under neutral equilibrium population models. GNDRY a/a and GNDRY +/+ refer to GOUNDRY 2La^a/2La^a and 2La⁺/2La⁺, respectively.

Figure 3

The distribution of Fay and Wu’s H and the Ewens-Watterson statistic for all loci in A. gambiae populations. Normalized Fay and Wu’s H and Ewens-Watterson statistics were calculated on synonymous variation at each gene and evaluated by the use of coalescent simulations (see Materials and Methods). Red dots indicate that the HEW statistic is significant at a 5% threshold level after correcting for multiple tests. Each panel presents the results for a different group and the gene names are presented next to the corresponding data point.

Figure 4

LD plot and neighbor-joining tree of FBN32 in GOUNDRY 2La^a/2La^a. (A) Linkage disequilibrium (r²) plotted among variant sites in the sequenced fragment of FBN32. Each pixel represents an r² value according to the shade of gray, as indicated in the scale. Greater r² values indicate increased linkage among those sites. The exon structure of FBN32 is placed on the diagonal of the plot to indicate the physical location of each variant site. The frequency of the derived allele frequency (DAF) relative to A. merus is plotted in the barplots above and to the side of the LD plot. (B) Neighbor-joining tree of all FBN32 sequences from GOUNDRY 2La^a/2La^a as well as three additional taxa: A. merus, A. arabiensis, and A. quadriannulatus. The scale bar indicates genetic distance. Large clades of genetically similar taxa were collapsed for presentation.

Figure 5

LD plot and neighbor-joining tree of Toll9 in GOUNDRY 2La⁺/2La⁺. Linkage disequilibrium (r²) plotted among variant sites in the sequenced fragment of Toll9. Each pixel represents an r² value according to the shade of gray, as indicated in the scale. Greater r² values indicate increased linkage among those sites. The exon structure of Toll9 is placed above and beside the plot to indicate the structural location of each variant site. The frequency of the derived allele relative to A. merus is plotted in the barplots above and to the side of the LD plot. The red bars indicate the three nonsynonymous sites and the triangle delineates the block of linked sites. (B) Neighbor-joining tree of all Toll9 sequences from GOUNDRY 2La⁺/2La⁺ as well as three additional taxa: Anopheles merus, Anopheles arabiensis, and Anopheles quadriannulatus. The scale bar indicates genetic distance. Large clades of genetically similar taxa were collapsed for presentation.

Figure 6

Clade specific patterns of divergence at Toll9 intron. Sliding window analysis of Jukes-Cantor corrected divergence (K_JC, range 0 to 1) at all sites relating the two haplotype groups (A and B) to three sister species. Divergence was calculated for 50-bp physical windows shifting 10 bp for every consecutive window. The top horizontal dotted line indicates the average maximum divergence per window for Toll9 in the M and S molecular forms as well as GOUNDRY 2La^a/2La^a. The lower horizontal dashed line indicates the average maximum divergence per window for all other sequenced loci that included an intron in the GOUNDRY 2La⁺/2La⁺ population. The legend indicates the color and line style for each clade/sister taxa comparison. The schematic under the plot depicts the exon structure in the sequenced region of Toll9.