Effective population size of Anopheles funestus chromosomal forms in Burkina Faso.

BACKGROUND: As Anopheles funestus is one of the principal Afro-tropical malaria vectors, a more complete understanding of its population structure is desirable. In West and Central Africa, An. funestus population structure is complicated by the coexistence of two assortatively mating chromosomal forms. Effective population size (Ne) is a key parameter in understanding patterns and levels of intraspecific variation, as it reflects the role of genetic drift. Here, Ne was estimated from both chromosomal forms, Kiribina and Folonzo, in Burkina Faso.
METHODS: Short-term Ne was estimated by evaluating variation at 16 microsatellite loci across temporal samples collected annually from 2000-2002. Estimates were based on standardized variance in allele frequencies or a maximum likelihood method. Long-term Ne was estimated from genetic diversity estimates using mtDNA sequences and microsatellites.
RESULTS: For both forms, short-term and long-term Ne estimates were on the order of 10(3) and 10(5), respectively. Long-term Ne estimates were larger when based on loci from chromosome 3R (both inside and outside of inversions) than loci outside of this arm.
CONCLUSION: Ne values indicate that An. funestus is not subject to seasonal bottlenecks. Though not statistically different because of large and overlapping confidence intervals, short-term Ne estimates were consistently smaller for Kiribina than Folonzo, possibly due to exploitation of different breeding sites: permanent for Folonzo and intermittent for Kiribina. The higher long-term Ne estimates on 3R, the arm carrying the two inversions mainly responsible for defining the chromosomal forms, give natural selection broader scope and merit further study.

Background

The efficient application of malaria control methods that target the mosquito vector depends upon knowledge of its population genetic structure. This information can improve current insecticide-based strategies and aid in the management of insecticide resistance, but it is also essential to future genetic control or modification strategies that aim to reduce, eliminate or replace vector populations with non-vectors. Unfortunately, present understanding of the population structure of any malaria vector is insufficient to underpin a genetic control programme, and nowhere is this shortfall more critical than in sub-Saharan Africa where three widespread species (Anopheles gambiae, Anopheles arabiensis and Anopheles funestus) are responsible for transmitting most of the 1–3 million fatal cases each year [1].

An. funestus, one of the most anthropophilic vectors known, exploits permanent or semi-permanent breeding sites such as marshes or rice fields. Its population density peaks in the dry season, extending malaria transmission by relay after An. gambiae and An. arabiensis populations have declined [2]. A highly polymorphic species, its population structure appears quite shallow across continental Africa. Evidence from microsatellite and mtDNA markers suggests a division between populations on either side of the Great Rift Valley complex, but little differentiation among populations within these regions, even between locations spanning several thousand kilometers [3]. However, in Burkina Faso, West Africa, analysis of polymorphic chromosomal inversions has revealed cryptic complexities in population structure. A temporally and spatially stable pattern of inversion polymorphism in sympatric and synchronous samples of An. funestus is inconsistent with random mating, suggesting the presence of two assortatively mating chromosomal forms designated Folonzo and Kiribina [4,5]. The relative abundance of Kiribina near rice cultivation and Folonzo near marshy areas suggests some partitioning of larval habitats, and the latter form is more likely to rest indoors, feed on humans, and be infected with malaria parasites [4]. At the molecular level, differentiation between forms is slight and is accounted for mainly by markers mapping to the chromosome arm (3R) bearing the principal inversions whose frequencies define the forms [6,7]. The working hypothesis is that Folonzo and Kiribina forms of An. funestus are incipient species whose distinctions are linked, at least in part, to chromosomal inversions.

Effective population size (Ne) is a central parameter in the description of population structure, though notoriously difficult to estimate with precision [8]. Ne is inversely related to genetic drift, the rate of fluctuation in allele frequencies caused by random sampling [9], and can be defined as "the size of an ideal population that exhibits the same rate of drift as the actual population it characterizes" [10]. In an ideal population, one which is randomly mating, constant in size and uniform in reproductive potential, Ne is equal to N, the census population size. In actual populations, Ne is usually less than N, because of large variance in reproductive success among individuals. This is especially true where population size fluctuates seasonally, as is the case for the primary malaria vectors of Africa, because Ne approximates the harmonic mean of single-generation sizes and, therefore, more closely reflects the lowest values. Ne has been estimated for An. gambiae and An. arabiensis by a number of authors in different parts of Africa using indirect genetic methods and direct mark-release-recapture [11-16]; values were typically at least 103, inconsistent with seasonal bottlenecks. Because no estimates were available for An. funestus, Ne was estimated for Folonzo and Kiribina forms in Burkina Faso using temporal variation at 16 microsatellite loci across 3 consecutive years, and sequence variation in an 834 bp region of the mitochondrial DNA ND5 gene.

Materials and methods

Mosquitoes

Indoor-resting insecticide spray-sheet collections were performed in December of three consecutive years (2000, 2001, and 2002) in Koubri-Kuiti, Burkina Faso (12°11'N; 1°23'W), as previously described [7]. Karyotyping and chromosomal form identification followed methods of Guelbeogo et al .[5].

Microsatellites and mtDNA

Genomic DNA was extracted from single mosquito carcasses [7]. Prior to analysis by PCR, genomic DNA was diluted 1:10 in H2O (~5 ng/ul). Morphological identification of An. funestus was verified on each specimen by a modified rDNA-based PCR assay [7].

For 2001, microsatellite data of Michel et al. [7] were taken from 50 randomly chosen specimens of both chromosomal forms. For 2000 and 2002, the same 16 physically mapped microsatellites were PCR amplified from ~50 specimens of each form (exact sample size given in [Additional file 1]). Products were diluted, pool-plexed (two groups of eight loci each), genotyped on a Beckman-Coulter CEQ8000, and sized with software provided (see [7] for detailed methods). Fstat2.9.3.2 [17] was used to estimate allelic richness (Rs), deviations from Hardy-Weinberg equilibrium (inbreeding coefficient, FIS) and linkage disequilibrium. The impact of suspected null alleles was explored by using MICRO-CHECKER [18] to adjust allele and genotype frequencies based on the Brookfield 2 estimate [19], followed by repeated analyses on the null allele-adjusted data set, as described previously [7]. No significant changes in outcome were observed. Microsatellite differentiation (FST) between Folonzo and Kiribina samples from each year was computed with Microsatellite Analyzer [20].

Mitochondrial sequences (834 bp of the ND5 gene) were taken from Michel et al. [6], based on 90 Folonzo and 96 Kiribina samples collected from La and Pehele, Burkina Faso in 2002 [GenBank, Kiribina: DQ126772–DQ126867; Folonzo: DQ127048–DQ127137]. There was no significant mtDNA or microsatellite differentiation between these locales for samples of the same chromosomal form [6].

Estimates of short-term Ne

Short-term Ne was estimated using temporal changes in microsatellite allele frequencies. Under the assumption of no mutation, selection or migration, changes in allele frequency are the result of genetic drift, whose strength is inversely related to population size. Both a moment estimator (based on the standardized variance in allele frequency change, F [21]) and a probability method (the maximum likelihood method of [22]) were used to calculateshort-term Ne from temporal samples. Although the probabilistic approach has been shown to have higher accuracy and precision [8], the moment method of Waples [21] (hereafter, F-statistic method) was the approach adopted for estimates in other vectors; including it facilitated comparison of Ne between vector species. The F-statistic method was implemented using the software programme NeEstimator [23]. The maximum likelihood method (hereafter, ML method) was implemented using the programme MLNE 2.0 [24]. To reduce computational burden, maximum population size was set initially at 15,000. Subsequently, analyses were repeated using a maximum size of 25,000. In both cases, the upper bound of the 95% confidence interval (CI) always reached the maximum limit for Folonzo (and did so for one sampling period for Kiribina). Thus, the upper bound CI was not determined in these cases. Following other authors [12-15], Ne was evaluated under the conservative assumption of 12 generations per year, but the effect of fewer (10) or more (20) generations per year was also explored.

Estimates of long-term Ne

Long-term Ne was estimated from current genetic variation (θ = 4Neμ for autosomal loci, where μ is mutation rate), under the assumption of mutation-drift equilibrium (MDE). For microsatellites, current genetic variation was represented by Xu and Fu's [25] estimator, θF, based on sample homozygosity under the single-step stepwise mutation model. The average mutation rate for dinucleotide microsatellite repeats is unknown for An. funestus, as for An. gambiae. In the latter species, an upper-bound estimate of 10-4 mutations per locus per generation based on the data of Zheng et al. [26,27] has been used [12], but this value is likely an overestimate. Accordingly, a slower mutation rate of 9.3 × 10-6 mutations per microsatellite locus per generation was adopted, based on the estimates derived from dinucleotide microsatellites in Drosophila melanogaster [28]. Mitochondrial values of θ were estimated from mean pairwise sequence differences per site (π [29]) and from the number of segregating sites among sequences (S [29]) as calculated in DNASP v 4.10.1 [30]. Assuming a mutation rate of 5.7 × 10-8 per base per gamete [31], long-term Ne was estimated using the equation θ = Neμ.

Results and discussion

Sample sizes and summary polymorphism statistics are presented in [Additional file 1]. Allelic richness and heterozygosity were relatively high across years and chromosomal forms (mean Rs per locus ranged from 4.4 to 21.0; mean Ho per locus ranged from 0.41 to 0.81). All loci were found to be in linkage equilibrium. Significant Hardy-Weinberg deficits occurred in 12 of 96 possible tests, but were clustered at 3 loci (AFND12, FunD, and AFUB12) suspected of harboring null alleles [6,7].

For all three sampling years there was slight but significant differentiation between the chromosomal forms, with frequency differences between loci on chromosome 3R (both inside and outside of inversions) accounting for much but not all of the differentiation (Table 1).

Table 1

Microsatellite FST values between An. funestus chromosomal forms at each of three sampling years

	Year 2000	Year 2001	Year 2002
All loci (16)	0.007***	0.006***	0.011***
Non-3R loci (11)	-0.013	0.005*	0.007**
3R Loci (5)	0.011***	0.007***	0.018***

*P < 0.05; **P < 0.01; ***P < 0.001

Short-term Ne

Table 2 presents estimates of short-term Ne across sampling intervals of one year (2000–2001, 2001–2002) and two years (2000–2002) for both chromosomal forms. The number of generations per year is uncertain, but 12 are plausible and the actual number should fall within the extreme values employed of 10 to 20. Between these values, all Ne estimates were on the order of 103, indicating the absence of seasonal bottlenecks. This is consistent with results from An. gambiae and An. arabiensis, even under severe climatic conditions, on islands, and following insecticide spray campaigns (e.g., [13,14]).

Table 2

Short-term Ne based on temporal changes in microsatellite allele frequency

	10 generations per year		12 generations per year		20 generations per year
Sampling Interval	FolonzoNe-Fa (95% CI)Ne-MLb (95% CI)	KiribinaNe-F (95% CI)Ne-ML (95% CI)	FolonzoNe-Fa (95% CI)Ne-MLb (95% CI)	KiribinaNe-F (95% CI)Ne-ML (95% CI)	FolonzoNe-Fa (95% CI)Ne-MLb (95% CI)	KiribinaNe-F (95% CI)Ne-ML (95% CI)
2000–2001	2,083 (665-∞)	811 (397-3,323)	2,500 (798-∞)	973 (476-3,987)	4,167 (1,330-∞)	1,622 (1,330-∞)
	2,736 (807-NDc)	1,081 (507-11,836)	3,248 (967-NDc)	1,288 (605-13,747)	5,300 (1,616-NDc)	2,120 (987-NDc)
2001–2002	2,166 (649-∞)	744 (357-3,324)	2,599 (779-∞)	893 (428-3,989)	4,332 (1,299-∞)	1,488 (714-6,648)
	1,876 (685-NDc)	588 (324-1,338)	2,227 (813-NDc)	597 (347-1,391)	3,637 (1,309-NDc)	982 (571-2,282)
2000–2002	7,363 (1,635-∞)	941 (543–1,990)	8,835 (1,963-∞)	1,129 (651–2,389)	14,726 (3,271-∞)	1,882 (1,086–3,981)
	3,404 (1,430-NDc)	865 (571–1,512)	4,064 (1707-NDc)	1,034 (682–1,807)	3,404 (1,430-NDc)	1,710 (1,127–2,988)

aNe estimated based on the standardized variance in allele frequency change (F) [21]

bNe estimated using the maximum likelihood (ML) approach implemented in MLNE 2.0 [24]

cND, not determined. See methods.

Regardless of the estimation method, Ne values for the Kiribina form at each time point were at least 2–3 times smaller than those for the Folonzo form, though the 95% CIs overlapped. The upper bound of the 95% CI for Folonzo was invariably infinity under the F-statistic method and was not determined under the ML method because it exceeded the maximum value imposed (25,000) during implementation of MLNE 2.0. With one exception (under the unrealistic case of 20 generations per year), the upper bound of the 95% CI for Kiribina was always defined, and ranged from ~1,300–14,000. Ecological data concerning the differences between these chromosomal forms is very incomplete, owing to the lack of molecular markers and the consequent necessity of determining karyotype from semi-gravid females by laborious cytogenetic methods. Preliminary indications are that Kiribina may prefer anthropogenic breeding sites such as rice fields, whereas Folonzo predominates in association with more natural and permanent habitats such as marshes and swamps. If confirmed by follow-up studies, the different breeding habitats may help explain differences in Ne, because rice fields do not produce anophelines continuously. Rice fields are flooded in June/July and harvested in October/November, with a second crop sometimes grown between March and June. Even during these intervals, the ability of An. funestus to exploit rice fields depends upon the stage of rice growth [32].

Estimates of Ne produced by the F-statistic and ML methods were very similar, except for the two-year sampling interval for the Folonzo form, where Ne values estimated from F were several-fold larger than those from the ML method and were inconsistent with both of the corresponding one-year sampling intervals. Unlike the ML method, the F-statistic method is known to have an upward bias if rare alleles are present, and it is possible that the longer sampling interval resulted in more rare alleles that led to an overestimate of Ne [22]. Because of the inconsistent performance of the F-statistic method and the greater precision of the ML method [8], ML estimates of Ne should be given greater weight. However, it should be noted that both methods assume the absence of mutation, selection and migration. While the first assumption is not unreasonable over short sampling intervals, and direct selection on microsatellite markers seems unlikely, the assumption of an isolated population without immigration may not be [24]. In the short-term, ignoring immigration, if it actually exists, leads to underestimates of Ne [8]. Both methods also assume constant population size and discrete generations, which are not realistic for An. funestus. Fluctuating population size with overlapping generations can lead to overestimates of Ne [33].

Long-term Ne

Unlike short-term Ne that reflects recent effects on genetic variation in a focal population, long-term Ne reflects evolutionary forces and demographic processes over a much greater geographic and historical frame (on the order of Ne generations; [8]), thus it approaches the effective size of the entire species. Long-term Ne was estimated from genetic diversity based on microsatellites and mtDNA sequences from samples collected in 2002 (Table 3). Microsatellite data from other years gave very similar estimates and, therefore, are not shown. As expected, long-term estimates are 2 orders of magnitude larger than the short-term Ne values (105 versus 103). Moreover, the estimates from microsatellites and mtDNA are roughly in accord. Departing from recent convention for estimates of Ne in other Afro-tropical vector species, 9.3 × 10-6 was assumed as the mutation rate for microsatellites, an order of magnitude slower than the upper-bound estimate for An. gambiae [12] but consistent with measurements from D. melanogaster [28]. The mtDNA mutation rate was assumed to be 5.7 × 10-8. The uncertainties in these mutation rates as applied to An. funestus mean that the long-term Ne values lack precision. It also should be noted that because of recent population expansion, An. funestus populations west of the Rift Valley (including those from Burkina Faso) violate the assumption of MDE under which long-term Ne is derived [3,6]; the departure from MDE is suggested by divergent estimates of θ from S or π seen in Table 3. The degree to which the long-term estimates of Ne are biased by population expansion is not known. However, long-term Ne estimates from Burkina Faso can be compared to long-term estimates from countries east of the Rift Valley where the evidence does not support a population expansion [3]. Using values of gene diversity derived from the data of Michel et al [3], who employed 10 of the 16 microsatellite loci used in the present study (only two on 3R), long-term Ne values from Malawi or Tanzania samples were smaller by about one-third (59,340 and 49,000) relative to values from Burkina Faso based on the same 10 loci (Kiribina: 160,940; Folonzo: 169,700).

Table 3

Long-term Ne based on genetic diversity estimates (θ) from 2002

	Folonzo		Kiribina
Locus	θ	Ne	θ	Ne
Microsatellites
All loci (16)	0.79a	250,590	0.77a	202,030
Non-3R loci (11)	0.76a	174,860	0.75a	174,420
3R loci (5)	0.85a	412,880	0.81a	262,780
mtDNA	0.00534b	93,684	0.00563b	98,771
	0.01939c	340,175	0.01774c	311,228

aθF [25]

bMean pairwise nucleotide difference per base, π [29]

cSegregating sites per base, S [29]

Comparison between Folonzo and Kiribina Ne values is not hindered by uncertainty in mutation rate, as this rate should be the same for both forms at a given set of loci. The Ne estimates from mtDNA are essentially identical between the forms. Across all 16 microsatellite loci, it appears that Ne for Folonzo is slightly larger than that for Kiribina. However, partitioning the loci into two groups – those residing on chromosome 3R or those that map elsewhere in the genome – reveals that the difference between forms can be explained by loci on 3R, the arm that carries the two main chromosomal inversions involved in distinguishing the forms [4,5]. There is higher diversity on this arm in Folonzo than in Kiribina. This is not altogether surprising, as the deterministic algorithm which defines these forms allows more polymorphism in Folonzo with respect to alternative arrangements on 3R. However, this is not the entire story. The relative level of genetic diversity at the 5 loci on 3R versus the 11 loci outside 3R is significantly higher in both forms (Mann-Whitney test: Folonzo U = 43.5, P < 0.03; Kiribina U = 49.5, P < 0.01). The reason(s) for this pattern are not clear from the present data, but higher genetic diversity provides more scope to natural selection and further investigation of chromosome 3R seems warranted.

Conclusion

The present estimates of effective population size for An. funestus, though approximate, preclude strong seasonal bottlenecks. Based on microsatellite variation between temporal samples of An. funestus in Burkina Faso, short-term estimates of Ne were on the order of 103 and were consistently smaller for the Kiribina than the Folonzo chromosomal form, possibly related to the preference of Kiribina for breeding in rice fields that are not in continuous production. Long-term estimates of Ne were consistent between classes of marker (microsatellites or mtDNA), and were two orders of magnitude larger than short-term estimates. Although long-term Ne did not differ between chromosomal forms, in each case there was significantly higher genetic diversity on the chromosome arm (3R) that carries the principal inversions defining the two taxa, a finding that merits further study.

As Ne reflects the strength of genetic drift, it is a central parameter in descriptions of population genetic structure. The Ne estimates derived here provide some insight into the complex population structure of An. funestus in Burkina Faso where the two chromosomal forms coexist. However, the geographic extent and nature of these forms has not been well characterized outside of Burkina Faso. In general, An. funestus is a heterogeneous species that occupies diverse habitats across Africa, where its population structure and history may differ. A more complete understanding of the forces that structure genetic variation in An. funestus will depend upon additional studies in other parts of its range.

Authors' contributions

APM participated in field collections, conducted the microsatellite genotyping, performed the data analysis, and drafted the manuscript. OG performed karyotype analysis and mtDNA sequencing and analysis. WG participated in collection efforts and performed karyotype analysis. N'FS, CC and NJB participated in the design and coordination of the study. NJB participated in drafting the manuscript. All authors read and approved the final manuscript.

Additional file 1

Microsatellite variation at 16 loci from two chromosomal forms of An. funestus during three consecutive sampling years in Burkina Faso, West Africa. The data provided in this table represent polymorphism statistics estimated at each of 16 microsatellite loci from population samples of An. funestus collected in Burkina Faso.

Click here for file