Population structure analyses and demographic history of the malaria vector Anopheles albimanus from the Caribbean and the Pacific regions of Colombia.

BACKGROUND: Anopheles albimanus is an important malaria vector in some areas throughout its distribution in the Caribbean and the Pacific regions of Colombia, covering three biogeographic zones of the neotropical region, Maracaibo, Magdalena and Chocó.
METHODS: This study was conducted to estimate intra-population genetic diversity, genetic differentiation and demographic history of An. albimanus populations because knowledge of vector population structure is a useful tool to guide malaria control programmes. Analyses were based on mtDNA COI gene sequences and four microsatellite loci of individuals collected in eight populations from the Caribbean and the Pacific regions of Colombia.
RESULTS: Two distinctive groups were consistently detected corresponding to COI haplotypes from each region. A star-shaped statistical parsimony network, significant and unimodal mismatch distribution, and significant negative neutrality tests together suggest a past demographic expansion or a selective sweep in An. albimanus from the Caribbean coast approximately 21,994 years ago during the late Pleistocene. Overall moderate to low genetic differentiation was observed between populations within each region. However, a significant level of differentiation among the populations closer to Buenaventura in the Pacific region was observed. The isolation by distance model best explained genetic differentiation among the Caribbean region localities: Los Achiotes, Santa Rosa de Lima and Moñitos, but it could not explain the genetic differentiation observed between Turbo (Magdalena providence), and the Pacific region localities (Nuquí, Buenaventura, Tumaco). The patterns of differentiation in the populations from the different biogeographic provinces could not be entirely attributed to isolation by distance.
CONCLUSION: The data provide evidence for limited past gene flow between the Caribbean and the Pacific regions, as estimated by mtDNA sequences and current gene flow patterns among An. albimanus populations as measured by MS loci which may be mainly influenced by semi-permeable natural barriers in each biogeographical region that lead to the genetic differences and effective population sizes detected. The relatively high genetic differentiation in the port city of Buenaventura may be the result of specific ecological conditions, human migration and activities and/or differences in effective population sizes. This knowledge could serve to evaluate and coordinate vector control strategies in these regions of Colombia.

Background

In Colombia, there are up to 100,000 malaria cases reported annually, approximately 10-20% of the cases in the Americas [1,2]. Anopheles albimanus is considered an important malaria vector in rural and periurban areas throughout its distribution in Colombia, along the Caribbean and Pacific Coasts and on San Andres Island [3]. A recent study showed that after 20 years it continues to be the predominant species in these coastal regions, accounting for 61 and 99% of the total capture of adult mosquitoes of the Caribbean and Pacific regions, respectively [4]. These regions represent substantial topographic, climatic and vegetation contrasts. In particular, the Caribbean coast is drier and has higher temperatures, while the Pacific coast is one of the rainiest regions globally, with high relative humidity and precipitation levels exceeding 5,000 mm/year [5]. The distribution of An. albimanus in Colombia is particularly remarkable because it includes coastal areas from three different biogeographic provinces, Maracaibo, Chocó and Magdalena in the neotropical Caribbean and Pacific regions of Colombia as described by Morrone [6].

The Pacific and the Caribbean regions are characterized by different levels of malaria transmission [1]. The Pacific region presents moderate to high transmission, while the Caribbean coast has presented low numbers of malaria cases and occasional epidemics [7]. Anopheles albimanus has been incriminated as a vector of Plasmodium falciparum and Plasmodium vivax, VK210 and VK247, in the Pacific region [4,8,9]. Despite its abundance in the Caribbean region and the high number of mosquitoes evaluated by Gutierrez et al [4], no infected individuals were detected, even though it was incriminated epidemiologically during a malaria outbreak in Guajira department in which An. albimanus was the predominant species [7]. This species is geographically variable in its infectivity by Plasmodium species throughout its distribution in the Americas [10-13].

Multiple factors could explain why An. albimanus has not been recently incriminated as a malaria vector in the Caribbean region: (i) because the number of human malaria cases reported from this region is lower than from the Pacific [1], more focused vector incrimination studies in this area are needed, (ii) Caribbean populations of An. albimanus could be more zoophilic because of extensive cattle ranching [3,5], (iii) different Anopheles species could be involved in malaria transmission at a local or regional level in the Caribbean coast [4,14,15], (iv) local adaptation between mosquito host and parasite could exist, as demonstrated in southern Mexico where reciprocal selection has led to local adaptation of P. vivax, and parasite populations are most compatible with their sympatric mosquito species (e.g., An. albimanus versus Anopheles pseudopunctipennis) [16], (v) the geographic, ecological and climatic differences between the regions could promote population differentiation of An. albimanus in Colombia and/or, (vi) demographic events could have influenced the current An. albimanus distribution and its role as a malaria vector in Colombia.

Interactions between changes in climatic conditions over large time scales with geographic features, sea level changes and contemporary factors such as human migration have influenced the distribution and diversification of different species [17-22]. For example, it has been suggested that the genetic differences detected in Anopheles darlingi, another important neotropical malaria vector, were affected by climate change at the end of the Pleistocene [23], Brazilian Amazon biogeography [24], and geographic barriers [25,26]. Similarly, studies of An. albimanus from Central and northern South America reported high variation in the intergenic spacers of nuclear ribosomal DNA [27] and differentiation between Central and South America using microsatellite loci and a mitochondrial DNA (mtDNA) marker [28,29]. Within Colombia, isoenzyme analyses revealed higher variability in Caribbean populations of An. albimanus than in Pacific populations and some loci showed significant allele frequency differences between these regions. However, cytogenetic data showed the same chromosomal banding patterns in all populations in the two regions, suggesting con-specificity in Colombia [30].

Several recent mtDNA studies have provided detailed accounts of species population structure and history [23,26,31-33], while studies based on highly polymorphic microsatellite loci have shown great potential for the evaluation of gene flow between populations of anophelines at a finer geographical and evolutionary scale [24,25,34-36]. Such studies have potential applied benefits, guiding and informing mosquito control strategies, such as the release of genetically modified mosquitoes refractory to the malaria parasite and the spread of insecticide resistance genes [37,38]. However, no mtDNA analyses or microsatellite loci information about malaria vector populations are available in Colombia. The taxonomic status of An. albimanus as a single lineage as defined by De Queiroz [39] is not disputed; however, additional detailed studies of inter-population genetic diversity could provide essential information about both historical and current distributions and processes affecting this species. Therefore, a focused population genetics and micro-evolutionary study of this important vector was conducted, using mtDNA cytochrome oxidase subunit I (COI) gene and four microsatellite markers, to address the hypothesis that historical and contemporary ecological and demographic processes may have led to some level of population differentiation of An. albimanus in the Caribbean and Pacific regions of Colombia. These data provide evidence for restricted past gene flow between the Caribbean and Pacific regions of Colombia as estimated by COI analysis, and contemporary gene flow among the regions as measured by four microsatellite markers (MS).

Methods

Sampling strategy

Adult An. albimanus specimens were collected from March 2005 to November 2007 in seven localities on the Caribbean and the Pacific coasts of Colombia (Figure 1, Table 1): Los Achiotes (ACH), Moñitos (MON), Santa Rosa de Lima (SRL), Nuquí (NUQ), Pizarro (PIZ), Buenaventura (BUE) and Tumaco (TUM). For details on sampling strategy and species identification see Gutiérrez et al [4]. Additional specimens collected in 2007 (Gutierrez et al, unpublished data) from Turbo (TUR), Antioquia, were also included (Figure 1, Table 1). Human landing catches of adults were conducted under an informed consent agreement using a protocol and collection procedures that were reviewed and approved by a University of Antioquia (SIU-UdeA) Institutional Review Board. ACH, SRL and MON are in the Maracaibo biogeographic province, NUQ, PIZ, BUE, and TUM are in the Chocó biogeographic province and TUR is in the Magdalena biogeographic province [6]. For COI sequencing samples from all sites were used, with individuals from Alto Guandipa in Mosquera (MTU) and La Ensenada in Santa Bárbara (STU), both sites in Tumaco, analysed as a single population (indicated as TUM; Figure 1A). For MS analysis, PIZ was excluded because there were not enough specimens with successful PCR amplification; STU and MTU from Tumaco were analysed separately because they are ~50 km apart (Figure 1B). Approximately 50% of the mosquitoes (randomly selected) included in this study were confirmed as An. albimanus using an ITS2-PCR-RFLP based assay [40,41].

Figure 1

A) Distribution of collection localities of . Altitude (meters) is based on a gray scale, and the distance scale is in km, B) Vegetation conditions of the Caribbean and the Pacific regions of Colombia, C) Temperature zones (according to Thornthwaite System) and moisture index, D) Average multi- annual precipitation (mm), E) Average multi- annual temperature (°C). Figures B to E were adapted from IGAC, with permission [5].

Table 1

Anopheles albimanus collection data

Department	Locality (Abbreviation)Collection site*	Location(longitude/latitude)	Collection date(month/year)
Caribbean Region
Magdalena	Los Achiotes (ACH)	11°15' N, 73°36' W	8/2005, 2, 6/2006
Bolívar	Santa Rosa de Lima (SRL)		7-10/2005, 2-3, 6/2006
	Cienaga	10°26' N, 75°21' W
	Hatillo	10°25' N, 75°22' W
Córdoba	Moñitos (MON)	9°15' N, 76°06' W	8-9/2005, 6-8, 11/2006
	Tierra Blanca	9°13' N, 76°08' W
	Santander de la Cruz	9°11' N, 76°10' W
Antioquia	Turbo (TUR)		11/2007
	Yarumal	8°07' N, 76°44' W
	Camerun	8°08' N, 76°43' W
Pacific Region
Chocó	Nuquí (NUQ)	5°42' N, 77°16' W	3-6, 8-9, 11/2005, 2-3, 6/2006
	Panguí	5°39' N, 77°18' W
	Panguí viejo	5°40' N, 77°17' W
	Pizarro (PIZ)	4°57' N, 77°21' W	6/2005
Valle del Cauca	Buenaventura (BUE)		6, 10/2005, 2, 8/2006
	Barrio La unión	3°51' N, 77°0' W
	La Barra	3°57' N, 77°22'W
	Puerto España	4°02' N, 77°26' W
Nariño	Tumaco (TUM)		12/2005, 10/2006
	Alto Guandipa	2°29' N, 78°26' W
	La Ensenada	2°27' N, 77°58' W

* For each locality the name of the collection site is given. Collectors: Naranjo, N.; Gutiérrez, L.; Cienfuegos, A.; Córdoba, L.; Solano, U.; Pinto, J.; Aviles, P.; Pozo, C.; García, L.; Quiñones, J. The specific collection sites are 2-11 km apart; two collection sites from TUM are ~50 km apart: Alto Guandipa, Mosquera (MTU) and La Ensenada, Santa Bárbara (STU).

DNA extraction, COI gene amplification and sequencing

DNA was extracted from individual mosquito abdomens following a salt precipitation protocol [42]. The DNA was then resuspended in 25 μL of TE buffer (10 mM Tris, 1.0 mM EDTA) and stored at -20°C. A 1,300-bp segment of mtDNA COI gene was amplified by PCR using primers UEA3 5'-TAT AGC ATT CCC ACG AAT AAA TAA-3' and UEA10 5'-TCC AAT GCA CTA ATC TGC CAT ATT A-3' [43]. The PCR was performed in 25 μL containing 2 μL of DNA, 0.2 μM each primers, 1× reaction buffer, 1.5 mM MgCl2, 0.1 mM each dNTP, and 1 U Taq polymerase (Bioline, London, UK). The cycling conditions were: 5 min denaturation at 94°C, followed by 36 cycles of 1 min denaturation at 94°C, 1 min annealing at 50.2°C, and 1 min 15 seconds extension at 72°C, ending with a final extension at 72°C for 7 min. The PCR products were purified using the Wizard SV Gel and PCR Clean-Up System (Promega, Madison, WI), following the protocol recommended by the manufacturer. Both strands of the purified PCR products were sequenced, with a region of ~600-bp overlap, for 220 individuals. Only DNA sequence segments with higher than 30 Phred values [44,45] were used in the analyses. The forward and reverse chromatograms were manually corrected using the electropherogram viewer Chromas Lite© [46]. Sequences from each mosquito were assembled by pairwise alignment using BioEdit Sequence Alignment Editor [47], then multiple alignment was performed using ClustalX [48]. Unique haplotypes were determined using DAMBE, version 4.5.68 [49]; identical sequences were considered to be a single haplotype.

Microsatellite genotyping

The An. albimanus mosquitoes genotyped from each locality were selected randomly. Four di-nucleotide microsatellite (MS) loci described for An. albimanus by Molina-Cruz et al [29], were genotyped: 6-41, 1-90, 2-14 and 2-25. Each locus was amplified by PCR performed in a 25 μL volume containing 1 μL of DNA, 0.25 μM of each primer, 1× reaction buffer, 2.5 mM MgCl2 (6-41, 1-90, 2-14) or 2.0 mM MgCl2(2-25), 0.2 mM for each dNTP, and 1 U Taq polymerase (Bioline, London, UK). PCR additives were used as follows: 5% DMSO for reactions 1-90 and 2-25 and 0.5 μg/ml BSA for reaction 2-14. For all reactions, the cycling conditions were: 5 min denaturation at 95°C, followed by 35 cycles of 30 sec denaturation at 94°C, 20 sec annealing at 62°C (6-41), 50°C (1-90), 55°C (2-14), 58°C (2-25), and 30 sec extension at 72°C, ending with a final extension at 72°C for 10 min. Amplified fragments were separated by electrophoresis on DNA denaturing 6% polyacrylamide sequencing gels, and MS alleles were visualized by silver staining. To estimate allele sizes the length of bands was compared to a 10 bp DNA ladder (Invitrogen Inc., Carlsbad, CA, USA) over the migration/size of each MS allele using Quantity One® software (Biorad Laboratories, Hercules, CA, USA). The most frequent homozygote MS alleles per locus were sequenced to provide reference sizes to estimate the number of repeats of other alleles. Allele sequences are available in GenBank [GenBank: FJ785408-FJ785419].

Data analyses

Descriptive statistics

Indices of population and overall genetic diversity for An. albimanus were determined per locus using both haplotype frequencies of mtDNA COI gene and allelic frequencies of MS loci. Haplotype and nucleotide diversities were generated using DnaSP version 4.50.2 [50] and Arlequin version 3.11 software [51]. Nucleotide composition and patterns of nucleotide substitution were characterized using MEGA version 4.0 [52]. For the complete COI data set, an appropriate model of nucleotide substitutions was determined using the programme Modeltest 3.8 [53,54].

MICRO-CHECKER 2.2.3 software [55] was used to detect potential errors that may occur at each MS locus during genotyping or the interpretation of data such as null alleles, stuttering and large allele dropout. Allele and genotype frequencies of the amplified alleles were compared and adjusted if necessary for population genetic analysis. Number of alleles (Na), expected heterozygosity (H), observed heterozygosity (H), allele richness (Rs) and Hardy-Weinberg Equilibrium (HWE) were estimated for MS loci using GenAlEx version 6.1 [56] and FSTAT v 2.9.3.2 [57]. Statistical significance for HWE and linkage disequilibrium (LD) for each pair of loci was assessed using exact probability tests available in GENEPOP version 4.0 [58]. Whenever multiple comparisons were carried out simultaneously, the sequential Bonferroni procedure [59] was applied.

Population differentiation test

Using COI and MS data, genetic differentiation between populations was measured by estimating the fixation index FST, using Arlequin software, and the significance was tested by permutation tests (10,000 replicates). Inbreeding coefficient (F) for each locus and overall loci by population [60] and the number of migrants per generation (Nm) between localities were also calculated in Arlequin. Effective population size (Ne) was estimated for each population under two methods, heterozygosity excess (HE) and linkage disequilibrium (LD) using NeEstimator software version 1.3 [61]. Geographical coordinates and distances of each sampling location were obtained using Google earth® software [62], and the programme SAMOVA 12.02, Spatial Analysis of Molecular Variance, was used to maximize the proportion of total genetic variance due to differences between groups of populations as well as identifying possible genetic barriers between them, without prior information of the sampling locations as is necessary for AMOVA [63]. Analysis of Molecular Variance (AMOVA) was used to examine population variation among, within, and between collection sites. The probability that a sampled individual belonged to each reference population was estimated using assignment statistics in Geneclass 2.0 [64]. The Bayesian method of Rannala and Mountain [65] was selected as the computation criterion and the re-sampling algorithm; the method was performed with a minimum number of simulated individuals of 10,000 and a type I error of 0.01. In addition, the BAPS 5 programme (Bayesian Analysis of Population Structure) was used for Bayesian inference of genetic structure and admixture analyses performed for the An. albimanus populations [66]. A dendogram based on MS genetic distances was constructed using the Unweighted Pair Group Method Arithmetic average (UPGMA) cluster analysis in TFPGA programme version 1.3 [67] to test the genetic relationship among different populations. The correlation between genetic distance estimated from MS and geographical distances, assuming isolation by distance, was evaluated by the regression FST/(1-FST) on the natural logarithm (Ln) of pairwise geographical distances between collection sites. Significance of the correlation coefficient was assessed applying the Mantel test (10,000 randomizations) [68] using the Isolation By Distance Web Service version 3.15 [69].

Demographic inference and neutrality tests

DnaSP and Arlequin software were used to test the hypothesis that all mutations on COI gene are selectively neutral [70], employing Tajima's D [71] and Fu's F[72]. Confidence intervals were tested by 10,000 coalescent simulations. Analyses for constant size and for growing populations were carried out to determine the distribution of the observed pairwise nucleotide site differences or mismatch distribution and the expected values (at equilibrium for no recombination) in a stable population or in growing and declining populations. Statistically significant differences between observed and model distributions were evaluated with the Sum of Squared deviation (SSD) and the raggedness statistic (r) [73] to reject the hypothesis of demographic expansion. Time since the population expansion can be estimated from t = τ/2μ, where τ (tau) is the date of the growth or decline measured in units of mutational time [τ = 2μt; t is the time in generations, and μ is the mutation rate per site (sequence size) and per generation] [74,75]. To analyse patterns of An. albimanus population history, haplotype networks were estimated by a parsimony-based method, which calculates the maximum number of mutational connections between pairs of sequences by the 95% parsimony criterion using the TCS computer software [76]. To obtain the most likely connection between two haplotypes and resolve some ambiguous loops in the network, several recommendations were followed [77-79]. Also, data were analysed using a neighbour-net network, which constructs split networks from inferred distance matrices, in the computer programme SplitsTree4, version 4.10 [80] and a neighbour-joining (NJ) algorithm was used to generate a tree in PAUP Version 4.0 [81] based on the genetic distance between haplotypes. For the latter method some analyses including an outgroup [82] were performed with COI gene sequences of An. darlingi collected from Colombia (unpublished data). The NJ tree and neighbour-net network analyses were performed with 1,000 bootstrap replications under the probability model identified using Modeltest. MS data were also used to estimate demographic processes such as recent population bottleneck and/or expansion, and heterozygosity tests were used to analyse deviations from Mutation-Drift Equilibrium (MDE) for each sample across all loci. At selectively neutral loci, the expected heterozygosity calculated from allele frequencies data (He) assuming HWE, and from the number of alleles and sample sizes (Heq), assuming a population at MDE, are expected not to be significantly different. Thus, if a significant number of loci show He>Heq, this indicates that the population recently experienced a bottleneck, conversely, a He83].

Results

Intra-population genetic diversity

A 1,058-bp sequence of the An. albimanus COI gene corresponding to positions 1,802-2,859 of Anopheles gambiae s.s. mitochondrion complete genome and aligned with RefSeq NC 002084 [84] was analysed for 220 An. albimanus mosquitoes from eight localities (Figure 1A, Table 2). All 112 haplotype sequences are available in GenBank under the accession numbers: FJ015158-FJ015269. Sequences contained no missing data such as ambiguous base pairs and alignment revealed 107 variable sites, 67 of which were parsimony informative. A total of 115 nucleotide substitutions (with a majority predicted to introduce synonymous amino acid changes), 101 transitions and 14 transversions were observed. The best-fit DNA substitution model selected by the hierarchical likelihood ratio test (hLRT) was TrN+I+G [85], with invariable sites (I = 0.7675) and Gamma distribution shape parameter (G: 0.9555). Using Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC), the best-fit model was the Transitional model: TIM+I+G [86], with invariable sites (I = 0.7664) and Gamma distribution shape parameter (G: 0.8906), A-T composition of the COI sequences was 39.7% and 29.5%, respectively, similar to that found in other Anopheles species [23], and the above parameters were adjusted to construct the split networks from inferred distance matrices and the neighbour-joining analyses.

Table 2

Description of shared COI haplotypes and statistics of genetic polymorphisms for An. albimanus

Population	Haplotypes	n	S	h	Hd (SD)	Pi (SD)
ACH	A(3), C(9), F(2),	24	18	13	0.8551 +/- 0.065	0.0028 +/- 0.002
SRL	A(10), C(2), F(2), G(1), N(1),	30	26	19	0.8920 +/- 0.052	0.0031 +/- 0.002
MON	A(8), C(1), F(1), G(2), O(1), P(1)	29	27	20	0.9261 +/- 0.041	0.0030 +/- 0.002
TUR	A(6), C(1), G(1), H(2), M(1), N(1), O(1), P(1)	29	28	22	0.9581 +/- 0.028	0.0038 +/- 0.002
Caribbean Region	A(27), C(13), F(5), G(4), H(2), M(1), N(2), O(2), P(2)	112	58	61	0.9274 +/- 0.018	0.0035 +/- 0.002
NUQ	B(8), D(4), E(2), H(2), K(2), L(1)	30	50	15	0.9103 +/- 0.037	0.0103 +/- 0.005
PIZ	B(5), D(4), E(1), I(3), J(2), L(1)	24	42	14	0.9275 +/- 0.033	0.0097 +/- 0.005
BUE	B(11), E(1), I(1), J(1), M(1)	30	42	15	0.8483 +/- 0.056	0.0088 +/- 0.005
TUM	E(1), K(1),	24	33	19	0.9710 +/- 0.024	0.0078 +/- 0.004
Pacific Region	B(24), D(8), E(5), H(2), I(4), J(3), K(3), L(2), M(1)	108	72	53	0.9396 +/- 0.016	0.0098 +/- 0.005
All localities	-	220	107	112	0.9666 +/- 0.006	0.0152 +/- 0.007

n, number of mosquitoes sequenced; S, number of segregating sites; h, number of haplotypes; Hd, haplotype diversity; Pi, nucleotide diversity; (SD), Standard deviation. Number in parentheses represents the frequency of each haplotype by region or locality. Underlined haplotypes are shared between Caribbean and Pacific populations.

Although diversity values were low overall, Pacific region populations had higher nucleotide diversities (0.0098, range 0.0078-0.0152) than Caribbean populations (0.0035, range 0.0028-0.0038), and haplotype diversity values were similar. There were 112 unique haplotypes (Table 2), 16 shared between different localities (named in order of frequency using capital letters), with the remainder exclusive to a particular geographic site (Figure 2, Table 2). In concordance with the parsimony-based network, the distance-based haplotypes networks (data not shown) and neighbour-joining tree (Figure 3) illustrated mainly two distinctive groups, corresponding to samples from the Caribbean or the Pacific regions.

Figure 2

Parsimony-based haplotype networks of 112 haplotypes from 1058 bp of the . The color of the circles and numbers in parentheses depict the origin and frequency of each haplotype, respectively. Each black bar represents one mutational step; the dashed line and m.h represent missing haplotypes. Unique haplotypes from Caribbean populations are represented in the network as Los Achiotes (A), Santa Rosa de Lima (SR), Moñitos (M) and Turbo (Y and C), and haplotypes from Pacific populations as Nuquí (N), Pizarro (P), Buenaventura (B) and Tumaco (T).

Figure 3

Neighbour-Joining tree based on 112 haplotypes from 1058 bp of the . The percentage of trees in which the associated haplotypes clustered together in the bootstrap test (1,000 replicates) is shown next to the branches, retaining only groups with frequency ≥75%. Shared haplotypes (M, H) from the Caribbean and Pacific populations are represented in the tree.

For MS analyses, 280 An. albimanus were genotyped and the four microsatellite loci were polymorphic in all collection sites studied (Figure 1B). The number of alleles per locus ranged from three to twelve: locus 6-41 showed the lowest value and 1-90 the highest (Table 3). Allelic richness (Rs) per locus ranged between three in MON and MTU (6-41) to eleven, also, in MTU (2-25). STU showed the lowest average allelic richness (6.125), while BUE showed the highest (8.144) compared to all other localities. The expected heterozygosity (He) across all samples ranged from 0.503 (ACH) to 0.87 (MTU), while the average expected heterozygosity ranged from 0.681 (STU) to 0.784 (TUR). Of 32 tests performed for HWE, 12 (37.5%) remained significant after the sequential Bonferroni correction (p < 0.01); all of them were associated with heterozygote deficits, in each locus as follows: 6-41 (STU), 1-90 (TUR and BUE), 2-14 (ACH, MON, NUQ and BUE) and 2-25 (ACH, SRL, TUR, NUQ and BUE). Of 48 tests conducted, no locus was at linkage disequilibrium after the sequential Bonferroni correction (p < 0.05) (Table 3). The frequency of null alleles at each locus was compared; however, there were no significant changes in comparison with the initial results. Under LD the average Ne calculated for all populations was 308 individuals (CI 215-493) and the average Ne for all Caribbean populations (including TUR) was an infinite number of individuals (CI 1,613- ∞), while for all Pacific populations it was 256 individuals (CI 136- ∞). SRL showed the lowest value (231 individuals), and TUR, NUQ, STU and MTU showed an infinite result (∞). Under the HE model all of populations showed infinite values (See Additional file 1: Estimates of effective population size (Ne) and heterozygosity tests based on MS data for An. albimanus).

Table 3

Genetic variability estimates at MS loci in An. albimanus from three biogeographic provinces of Colombia

Locus		ACH	SRL	MON	TUR	NUQ	BUE	STU	MTU	All samples
		Maracaibo/Caribbean Region			Magdalena	Chocó/Pacific Region
	N	40	40	40	40	40	40	21	19	280
6-41	Na	4	4	3	4	4	5	4	3	7
	Rs	3.473	3.475	3.000	3.475	3.921	4.334	3.898	3.000	3.675
	He	0.503	0.617	0.627	0.662	0.583	0.623	0.484	0.555	0.659
	Ho	0.325*	0.425*	0.625	0.575	0.525	0.450	0.381**	0.263*	0.446
	FIS	0.365	0.323	0.016	0.143	0.113	0.289	0.236	0.545
1-90	Na	8	12	8	13	9	12	7	7	16
	Rs	7.075	9.702	7.681	10.63	8.077	9.526	6.904	7.000	8.866
	He	0.742	0.838	0.813	0.848	0.805	0.819	0.746	0.742	0.813
	Ho	0.650	0.800	0.7	0.800***	0.775	0.500***	0.714	0.632	0.683
	FIS	0.136	0.058	0.151	0.069	0.05	0.4	0.067	0.176
2-14	Na	8	9	8	11	10	10	6	8	14
	Rs	7.321	7.886	6.132	8.469	8.512	8.902	5.898	8.000	8.606
	He	0.811	0.807	0.674	0.792	0.802	0.830	0.717	0.794	0.812
	Ho	0.700***	0.575*	0.600***	0.825	0.500***	0.525***	0.857	0.579	0.604
	FIS	0.149	0.299	0.122	-0.029	0.387	0.379	-0.173	0.295
2-25	Na	10	9	10	9	9	11	8	11	12
	Rs	8.373	7.799	8.371	7.991	7.772	9.814	7.802	11.000	8.926
	He	0.832	0.829	0.833	0.834	0.802	0.860	0.776	0.871	0.858
	Ho	0.475***	0.600**	0.725	0.750**	0.375***	0.500***	0.571	0.789	0.558
	FIS	0.439	0.288	0.142	0.114	0.541	0.429	0.286	0.121
All loci	Na	7.5	8.5	7.25	9.25	8	9.5	6.25	7.25	7.94
	Rs	6.56	7.215	6.296	7.64	7.071	8.144	6.125	7.25	7.037
	He	0.722	0.773	0.737	0.784	0.748	0.783	0.681	0.741	0.746
	Ho	0.537	0.600	0.663	0.738	0.544	0.494	0.631	0.566	0.596

N: Number of mosquitoes, Na: number of alleles, Rs: Allele richness; H: expected heterozygosity, H: observed heterozygosity, F: inbreeding coefficient. All loci/samples: mean values over loci or collecting sites. Asterisks are showing significant heterozygote deficits according to exact tests against Hardy-Weinberg proportions after the sequential Bonferroni correction; *p < 0.05, **p < 0.01, ***p < 0.001.

Genetic structure and differentiation

In the AMOVA analyses using haplotype frequencies with the eight populations separated into two groups, including TUR in the Caribbean region, 5.85% (p < 0.05) of the total variance was explained at the among regions level. When populations were separated into two groups, including TUR in the Pacific region, only 3.86% (p < 0.05) of the total variance was due to the variation among regions (Table 4).

Table 4

Analysis of Molecular Variance (AMOVA) using COI and MS data in An. albimanus from Colombia

Source of variation	Variance components	Percentage of variation	FST
MS Two groupsAmong the Caribbean and the Pacific regions(including TUR in the Caribbean region)	0.04944	3.08	0.05040*
Among populations within regions	0.03140	1.96
Within populations	1.52313	94.96
MS Two groupsAmong the Caribbean and the Pacific regions(including TUR in the Pacific region)	0.05419	3.37	0.05162*
Among populations within regions	0.02871	1.79
Within populations	1.52313	94.84
COI Two groupsAmong the Caribbean and the Pacific regions(including TUR in the Caribbean region)	0.02925	5.85	0.08934*
Among populations within regions	0.01539	3.08
Within populations	0.45509	91.07
COI Two groupsAmong the Caribbean and the Pacific regions(including TUR in the Pacific region)	0.01916	3.86	0.08262*
Among populations within regions	0.02183	4.40
Within populations	0.45509	91.74

* p < 0.05; significant level based on 1,023 permutations

The Fvalue estimated from COI sequences for the comparison between the Caribbean and Pacific regions was 0.07 (p < 0.05), and most of the pairwise comparisons of Fbetween localities were significant; however, no significant Fvalues at p < 0.01 were observed in comparisons between SRL, MON and TUR from the Caribbean coast, or between NUQ, PIZ and BUE from the Pacific coast (Table 5). Estimates of pairwise genetic differentiation (F) and gene flow (Nm) among populations using MS data are shown in Table 6. Overall significant Fvalues were observed between population belonging to the same region, Caribbean or Pacific. Of 28 differentiation tests, 18 were significant (p < 0.01). Comparisons between ACH and STU showed the highest degree of differentiation. An. albimanus samples collected from TUR, located in the middle of the collection range, showed lower differentiation with its nearest localities: NUQ (Pacific) and MON and SRL (Caribbean). Nm estimates among An. albimanus populations ranged from 4 to 59 individuals, presenting an average value of 21 individuals.

Table 5

F-Statistics based on pairwise estimates of COI haplotype frequencies of An. albimanus from 8 localities in Colombia

Populations		ACH	SRL	MON	TUR	NUQ	PIZ	BUE	TUM
ACH		___
SRL		0.05820**	___
MON		0.06177**	-0.00873	___
TUR		0.05614**	0.00155	-0.00548	___
NUQ		0.11671**	0.09885**	0.08179**	0.06150**	___
PIZ		0.10870**	0.09058**	0.07319**	0.05698**	-0.00089	___
BUE		0.14841**	0.12989**	0.11293**	0.09591**	0.02299	0.03041*	___
TUM		0.08696**	0.06914**	0.05171**	0.03550**	0.05450**	0.04907**	0.09020**	___

**p < 0.01; * p < 0.05

Table 6

Estimates of pairwise genetic differentiation (F) and gene flow [48] among populations of An. albimanus based on MS data

Populations	ACH	SRL	MON	TUR	NUQ	BUE	STU	MTU
ACH	__	59	31	13	6	9	4	7
SRL	0.00834	__	56	38	8	11	5	10
MON	0.01602*	0.00889	__	26	10	15	5	11
TUR	0.03596***	0.01286*	0.01880**	__	47	22	12	24
NUQ	0.08121***	0.05584***	0.04891***	0.01046	__	18	41	35
BUE	0.05404***	0.04186***	0.03227***	0.02211**	0.02697**	__	10	35
STU	0.11405***	0.09112***	0.09023***	0.04015***	0.01188	0.04721***	__	15
MTU	0.06623***	0.04589**	0.04233***	0.02050*	0.01439	0.01421	0.03274*	__

Above diagonal: Nvalues; below diagonal: Fvalues; * p ≤ 0.05; ** p ≤ 0.01; *** p ≤ 0.001

The SAMOVA was used to estimate heterogeneity among populations from the Caribbean and Pacific regions. This test maximized the proportion of total genetic variance between populations at six groups: 1) ACH, SRL and MON, 2) TUR, 3) NUQ, 4) BUE, 5) STU, 6) MTU (Figure 4A). In agreement with SAMOVA data, the UPGMA dendogram based on genetic distances showed two distinctive groups in which ACH, SRL and MON constituted a cluster (82% support) and TUR, NUQ, BUE, STU and MTU were included in a second one at 61% support (Figure 4B). Given the number of MS loci tested in this study and in the low support values for the UPGMA dendogram, the possibility that data on the apparent relationship between populations from the Caribbean and Pacific regions may be biased cannot be excluded. However, BAPS clustering was also congruent with the results obtained by SAMOVA and UPGMA analyses, which proposed two groups without admixture among them, except for TUR which showed signs of admixture with both. AMOVA was conducted with both the COI and MS data to test the Isthmus of Panama as a putative barrier between the Caribbean and Pacific regions, including and excluding TUR. The AMOVA corresponding to two groups: the Caribbean and the Pacific samples (excluding TUR), showed the highest percentage of the total variance and Fvalue, 4.3% (p < 0.05) and 0.06049 (p < 0.05), respectively (Table 4). Results of the assignment statistics showed that on average 28.6% (80 of 280) of the individuals were correctly assigned to their original reference site. Samples from STU and BUE presented the highest proportion of correctly assigned individuals (38%). In general, mis-assignments occurred between samples of the same region, either the Caribbean or the Pacific (Table 7).

Figure 4

Cluster analysis based on MS data. A: Map of the sampling points and the barriers between the groups of populations defined by SAMOVA. B: UPGMA dendogram based on Nei's genetic distance for the eight populations of An. albimanus.

Table 7

Data of assignment tests based on MS among samples of An. albimanus from Colombia

Populations	ACH	SRL	MON	TUR	NUQ	BUE	STU	MTU
ACH	0.25*	0.20	0.18	0.08	0.05	0.08	0.03	0.15
SRL	0.20	0.20*	0.20	0.13	0.10	0.05	-	0.13
MON	0.18	0.10	0.35*	0.05	0.13	0.10	0.03	0.08
TUR	0.03	0.18	0.10	0.33*	0.18	0.05	0.08	0.08
NUQ	0.05	0.03	0.10	0.23	0.25*	0.05	0.15	0.15
BUE	0.10	0.03	0.08	0.05	0.10	0.38*	0.15	0.13
STU	-	-	-	-	0.33	0.10	0.38*	0.19
MTU	0.11	0.16	0.11	-	0.11	0.05	0.37	0.11*

Values are proportions of mosquitoes from the original reference site (rows) assigned to each population site (columns), *Proportions of samples correctly assigned.

A correlation test between genetic distance based on MS loci, measured by FST/(1- FST), and the geographic distance (Ln) was statistically significant, although with low values, for the whole dataset (Mantel tests: R= 0.36, p = 0.003). The IBD model best explains genetic differentiation among populations from the Caribbean region: ACH, SRL and MON (excluding TUR), but it could not explain the low differentiation observed within the Pacific region populations (NUQ, BUE, STU and MTU). No statistically significant correlation was detected when tests of IBD were carried out separately for the Caribbean (including TUR) or the Pacific region (with or without TUR).

Demographic inference and neutrality tests

The most common haplotypes were A (n = 27), exclusively Caribbean, and B (n = 24), exclusively Pacific. Three localities (TUR, NUQ and BUE) share haplotypes H and M from both regions (Figure 2, Table 2). Haplotypes from the two regions differed by more than 13 mutational steps. Interestingly, the shared haplotypes (H, M) and the Pacific haplotypes J, P38, B357, B857 and N2865 clustered with the Caribbean region (Figure 2). Haplotypes A and B were the most common interior haplotypes and so are most likely to be ancestral [87]. The majority of haplotypes in both regions were tip alleles, considered to be more recently derived and geographically restricted [77,87]. In addition, private haplotypes located peripherally are suggestive of a demographic expansion [72,75] within the Caribbean region with subsequent limited gene flow between the two regions. In general, the parsimony-based network (Figure 2), distance-based haplotype networks and neighbour-joining tree (Figure 3) consistently show two distinctive groups, with a large star-shaped Caribbean network surrounding haplotype A, a smaller network around haplotype C, and a small Pacific network surrounding haplotype B. In the NJ tree, the Caribbean haplotypes were supported at 86% and the Pacific haplotypes at 98% (Figure 3), and the bootstrap support level was 99.9% for the branch that connects Caribbean and Pacific haplotypes in the An. albimanus neighbour-net network.

Both neutrality tests found the Caribbean populations have significant negative values, and one test was significant for the Pacific populations (see additional file 2: Results of neutrality tests based on COI sequences of An. albimanus from Colombia), which suggested a possible past demographic expansion event. The distribution of the number of differences (mismatches) between pairs of DNA sequences from Caribbean populations demonstrated the expected unimodal distribution for Caribbean all populations together and also for each population. The expected distribution did not differ significantly from the sudden-expansion model. The distribution for each locality from the Pacific region and for the grouped Pacific region showed a multimodal distribution typical of populations at equilibrium. Although the r (0.009; p = 0.91) and SDD (0.002; p = 0.75) values were not significant for Caribbean populations, the estimated values were small, also supporting a population expansion. As in previous Anopheles studies [23,31], the Drosophila (D. melanogaster and D. yakuba) mutation rate values of 10-8/site/year [88] and 10 generations/year [89] were assumed to estimate the time for the Caribbean population expansion. For An. albimanus, using τ = 4.654, this is approximately 21,994 years ago (95% CI, 8,969-34,347) during the late Pleistocene.

The heterozygote test performed under SMM, IAM and TPM using MS data showed different results. Using IAM, heterozygote excess was detected for all of the populations tested; however, none were statistically significant. Under SMM four populations showed heterozygote deficits: STU, BUE, NUQ and ACH; only STU was statistically significant. For TPM three populations presented heterozygote excess: MON, SRL and NUQ; only NUQ was statistically significant. STU showed heterozygote deficits that were not statistically significant.

Discussion

Distinctive evolutionary genetics and demographic history evidence

Nucleotide diversity values for An. albimanus COI sequences, including synonymous and non-synonymous changes, were similar to those found for An. darlingi COI sequences from Central and South America [23], An. gambiae [84] and other species within the Insecta class. In particular, all are characterized by having an adenine-thymine rich mitochondrial genome. Anopheles albimanus COI sequences presented a combined frequency of ~70% AT. The results show that An. albimanus populations from the Caribbean and Pacific regions in Colombia have moderately high genetic diversity, in contrast to the lower diversity in Pacific localities detected using isozymes [30]. This difference is likely a result of the greater sensitivity of DNA markers to detect genetic variability, since DNA polymorphic sites are not necessarily seen at the protein level [90]. Nucleotide diversity was higher in most of the Pacific populations; this may indicate that different An. albimanus maternal lineages invaded these two Colombian regions at different times and (or) that the Pacific populations are older relative to the Caribbean ones, perhaps most closely related to populations from Cuba and Central America, which had similar values using ND5 sequence [29].

Panmictic populations of An. albimanus (within ~665 km) had been observed in Central America, Costa Rica and Panama, Colombia and Venezuela. It appears that An. albimanus populations have had different periods of isolation followed by secondary contact throughout the species range in America [27-29]. The data showed some haplotypes from Pacific populations clustering with the Caribbean group (Figure 2), possibly a genetic signature of a panmictic gene pool that existed before the late Pleistocene. Multiple factors may be responsible for the observed distribution and frequency of haplotypes from the Caribbean and Pacific regions in Colombia. Most haplotypes are not shared between regions, perhaps because of the distance between populations (>200 km), distinctive demographic history, human impact (insecticide use could have led to local haplotype extinctions) and/or distinctive ecological conditions in each region. Previous mtDNA analysis of An. albimanus did not detect a significant correlation between genetic and geographic distance [29]. Although, genetic differences and geographical distances were not directly compared because the populations were not in MDE, distance did not appear to be a significant factor affecting the data.

The neutral model is rejected for the Caribbean populations, indicating possible previous population expansion and/or natural selection. In general, all tests illustrate two distinctive groups, corresponding to haplotypes from the Caribbean and from the Pacific regions. The data did not show significant differences between An. albimanus populations from TUR, located in the Magdalena province as described by Morrone [6], and the other populations from the Caribbean coast. The relatively strong mtDNA support for expansion in the Caribbean An. albimanus populations may indicate different demographic histories for this species in these two regions of Colombia, and the AMOVA support for variance resulting from these two groups is significant, albeit low. Population expansion for An. albimanus from the Caribbean coast in Colombia is estimated to the late Pleistocene, similar to An. darlingi [23].

Contemporary population structure

Microsatellite loci used in this study have not been physically mapped in An. albimanus polytene chromosomes, consequently, their location on chromosome inversions is unknown and neutrality cannot be assumed [91]. Chromosomal inversions have not been detected in An. albimanus and it is conspecific along its distribution [30,92]; however, the microsatellite loci analysed in this study were highly polymorphic (Table 3), showing their usefulness for evaluating the population structure of An. albimanus in Colombia.

Significant departures from HW equilibrium were detected associated with heterozygote deficits (Table 3), similar to that reported in other anopheline microsatellite studies [25,29,34,93,94]. Heterozygote deficits are attributed to either significant inbreeding, Wahlund effect, natural selection or the presence of null alleles [24,95,96]. Inbreeding as the possible cause of heterozygote deficits is not considered due to the fact that it affects all loci equally, which is not compatible with the heterogeneity detected in this study (Table 3). The Wahlund effect refers to reduction of heterozygosity in a population caused by subpopulation structure [97]. In the data, SAMOVA, UPGMA and BAPS based on MS loci identified some degree of population subdivision, mainly among the Caribbean and the Pacific regions, and within the Pacific region, thus a part of the heterozygote deficits detected could be the result of Wahlund effect. High levels of heterozygote deficits could also be the result of null alleles, as a result of accumulation of mutations in the primer binding sites [34]. In this study, there was failure of amplification of different loci for some specimens; however poor DNA quality was discarded as a cause, similarly, Molina-Cruz et al [29], observed similar cases of amplification failure. In summary the most likely causes of heterozygote deficits in the present study were null alleles and Wahlund effect.

The moderate level of population differentiation detected (F) was not observed with isozyme markers [30]. In theory, the analysis of a low number of loci could potentially increase population differentiation values (F), but at the same time, the Fis the estimate presenting the lowest differentiation bias [98]. The SAMOVA, UPGMA and BAPS detected two distinctive clusters: ACH, SRL, MON and TUR, NUQ, BUE, STU and MTU. Unlike the COI data, the TUR population clustered with the Pacific region populations. Low admixture was detected between the cluster represented by ACH, SRL and MON from the Caribbean region and the other five populations. Thus, the largest genetic differentiation was observed in comparisons among the Caribbean and the Pacific region, and SAMOVA analyses placed TUR, NUQ, BUE, STU and MTU in five different groups perhaps suggesting a specific barrier that reduces gene flow between these populations.

Data could suggest that features related to the three different biogeographic provinces, Maracaibo, Magdalena and Chocó [6], such as vegetation, weather, moisture, precipitation and temperature (Figure 1B-E) [5], constitute a semi-permeable barrier that reduces gene flow of An. albimanus populations at the inter-regional or the inter-biogeographic provinces level. Also the differences in Ndetected among the populations may have contributed to the genetic differentiation of these regions [99]. Because some differentiation was observed between populations in the Chocó province (Pacific region) it is possible that other processes besides ecological conditions may affect free gene flow between An. albimanus populations.

The level of differentiation (F), among the Caribbean and the Pacific region of Colombia is comparable to that reported among countries in South America and between Cuba and continental An. albimanus populations (F= 0.057 and 0.059, respectively) [29]. Similar levels of differentiation were also found among Anopheles nuneztovari s.l from Cordoba, Norte de Santander and Valle in the east, west and southeast of Colombia (F= 0.024-0.06), using RAPDs (Randomly Amplified Polymorphic DNA) [100]. In this study, the authors suggested IBD could explain their results. In a study of An. darlingi populations from Cordoba, Meta and Chocó in the northwest and west of Colombia, using AFLPs (Amplified Fragment Length Polymorphism) and RAPDs (RAPD F= 0.084, AFLP F= 0.229), results suggested that the observed genetic differences could be the result of the biogeographic characteristics of each particular region [101]. The data showed that in addition to the differences in the demographic history of each region, the presence of semi-permeable biogeographic barriers, could contribute to the differentiation observed using MS.

In general, in the Caribbean region, there was low genetic differentiation, partly explained by the high Ne, that could have increased gene flow and decreased population genetic structure [24]. IBD was the model that best explained differentiation among ACH, SRL and MON; however, if TUR was included in the analysis, the resulting data did not fit this model. TUR, the most genetically distant Caribbean region population, (Figure 4B), could be influenced by its location in Magdalena biogeographic province [6], and by differences in the effective population sizes. Similarly, low non-significant differentiation was observed among the Pacific region populations between NUQ-STU (362 Km), NUQ-MTU (371 Km) and BUE-MTU (217 Km) (Table 5) and the level of differentiation observed was not congruent with IBD.

A significant level of differentiation among the populations closer to BUE was detected, probably influenced by ecological characteristics [5] that could reduce gene flow with NUQ and STU (Figure 1B-E), in addition to their differences in effective population sizes. Ecologic and climatic variation with appropriate conditions for mosquito development have been recognized among the main causes for peaks of mosquito abundance and subsequent peaks of malaria cases [102]. Therefore, further studies would be essential in BUE in comparison with others localities from the Pacific region. Even though there was restricted gene flow with respect to adjacent populations, BUE presented the highest genetic diversity, characterized by a high proportion of private alleles in low frequency. This could be explained by the characteristics of BUE, the main Colombian port city, where human migration and activities may promote gene flow among mosquito populations with different genetic pools, thereby increasing variability, as has been reported for An. gambiae [35,103,104]. The levels of differentiation observed between MTU and STU could also be due to ecological differences [5]; however, an overestimation of differentiation (F), is not discarded given the low sample number [97,105].

Most population structure analyses assume equilibrium, however, in several studies on anopheline vectors this has been violated, for example for An. gambiae and Anopheles arabiensis [106], Anopheles dirus [94], Anopheles moucheti [36] and An. darlingi [34]. In this study, the heterozygosity test did not detect a significant departure from equilibrium for most of the populations. Nevertheless, under the TPM model, there was evidence of a bottleneck for NUQ which was not sustained by the SMM and IAM models and in disagreement with the effective population size observed, perhaps because of the mutational model and low number of analysed loci [83]. Under the SMM model, the STU population appears to be expanding, consistent with the effective population size.

The Ne calculated under the HE model assumes equilibrium [107]. Some of the Caribbean and Pacific populations were at H-W disequilibrium for some loci; it is possible that the departure from equilibrium influenced the values detected with HE, nevertheless, the confidence intervals obtained with LD agreed with those of HE. In general, high Ne values were obtained for all localities and on average, the effective population size (Ne = 308 individuals) was three fold higher than in the study with An. albimanus from Central and South America (Ne = 96 individuals) [29]. For TUR and the Pacific region populations, except BUE, an infinite Ne was obtained under both models. According to Donnelly et al. [99] heterogeneity in the Ne could contribute to the levels of population genetic differentiation seen in this study. Similarly, Ne heterogeneity has been reported for An. darlingi [24,34] and within the An. gambiae complex [108-111]. In An. gambiae, the Ne heterogeneity observed in the savana cytotype was demonstrated to be region specific [111]. In the present study, high Ne might be associated with a wide availability of breeding places that are suitable for An. albimanus, such as lakes, ponds, mine excavations and fish ponds [3,112]. The high precipitation levels characteristic of the Pacific region [5], could contribute significantly to the availability of a variety of breeding sites for An. albimanus.

Conclusion

The genetic diversity of Caribbean and Pacific populations of An. albimanus detected using COI gene sequences appears to be mostly influenced by demographic events in An. albimanus populations during the late Pleistocene. The data indicated little to moderate genetic differentiation among An. albimanus populations from the Caribbean and the Pacific regions. Current gene flow patterns may be mainly influenced by semi-permeable natural barriers in each biogeographical region that lead to the genetic differences and effective population sizes found in this study.

Current hydrological and climatic variation in Colombia is a factor associated with malaria outbreaks for multiple localities throughout endemic regions of the country [113]. Therefore additional studies could test the influence of specific ecological and climatic conditions on the genetic differentiation in An. albimanus populations from the Caribbean and Pacific regions. It would also be of interest to examine additional populations of other Anopheles species to evaluate the generality of these patterns and to test the hypothesis that Caribbean expansion occurs among the Anopheles species whose distributions include these Colombian geographic areas.

In summary, this study showed current patterns of gene flow among the different An. albimanus populations from Colombia. This knowledge could be applied to existing vector control strategies since data on vector population genetics allow the inference of the spread of genes important for insecticide resistance or refractoriness to malaria parasites. In addition, this information contributes to the understanding of the epidemiology and the dynamics of disease transmission [99]. Moreover, effective population sizes calculated in this study may be used to evaluate local control measure effectiveness.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

LAG contributed in the design of the study, extracted DNA, performed all COI molecular procedures and analysis, co-participated in MS genotyping and analysis and drafted the manuscript. NJN contributed in the design of the study, was involved in field surveys, conducted procedures for microsatellite genotyping and analysis and helped in manuscript preparation. AVC co-participated in the design of the study and performed all molecular identification assays. CEM, SL and JEC co-participated in the design of the study and helped drafting the manuscript. MMC participated in the design of the study, data analysis, draft of the manuscript, general supervision of the research group and funding acquisition. All authors read and approved the final document.

Additional file 1

Estimates of effective population size (. The table provided include the estimated values of effective population size (Ne) and heterozygosity tests based on MS data for An. albimanus from eight sites from the Caribbean and the Pacific regions of Colombia.

Click here for file

Additional file 2

Results of neutrality tests based on . The table provided describe the estimated values of neutrality tests based on COI sequences of An. albimanus from eight sites from the Caribbean and the Pacific regions of Colombia.

Click here for file