Distinct population structure for co-occurring Anopheles goeldii and Anopheles triannulatus in Amazonian Brazil.

To evaluate whether environmental heterogeneity contributes to the genetic heterogeneity in Anopheles triannulatus, larval habitat characteristics across the Brazilian states of Roraima and Pará and genetic sequences were examined. A comparison with Anopheles goeldii was utilised to determine whether high genetic diversity was unique to An. triannulatus. Student t test and analysis of variance found no differences in habitat characteristics between the species. Analysis of population structure of An. triannulatus and An. goeldii revealed distinct demographic histories in a largely overlapping geographic range. Cytochrome oxidase I sequence parsimony networks found geographic clustering for both species; however nuclear marker networks depicted An. triannulatus with a more complex history of fragmentation, secondary contact and recent divergence. Evidence of Pleistocene expansions suggests both species are more likely to be genetically structured by geographic and ecological barriers than demography. We hypothesise that niche partitioning is a driving force for diversity, particularly in An. triannulatus.

Malaria remains one of the most important tropical diseases in the world with an estimated 216 million cases in 2010 (WHO 2011). Brazil has the largest incidence of disease and malaria related deaths in South America, averaging 300 thousand cases annually (MS 2008, Barreto et al. 2009). Nearly all cases occur in the Amazon Basin (Akhavan et al. 1999, Moutinho et al. 2011), where the relatively low levels of transmission have been difficult to control (da Silva-Nunes et al. 2011). The Amazon Basin is one of the world's most important bioregions, with a variety of habitats suitable for Anopheles species, many of which are malaria vectors.

Most primary vectors are widespread because of excellent colonising abilities and adaptation to variable environmental conditions (Costantini et al. 2009, Cohuet et al. 2010, Loaiza et al. 2012). Though not major malaria vectors in South America, Anopheles goeldii and Anopheles triannulatus are locally and regionally important (Rubio-Palis 1994, Póvoa et al. 2001, Galardo et al. 2007) with large geographic distributions and a moderate to wide habitat range (unpublished observations). Both species can colonise and increase abundance in altered or temporary environments (Tadei & Dutary-Thatcher 2000). Overlapping species distributions and mostly generalist behaviour (McKeon et al. 2013), in addition to the hypothesis of competition towards exclusion, noted by Tadei et al. (1998) make An. goeldii and An. triannulatus ideal comparisons.

An. goeldii (Rozeboom & Gabaldón), often mistakenly identified as Anopheles nuneztovari (Sallum et al. 2008 ), has been resurrected from synonymy (Calado et al. 2008) and is part of the Nuneztovari complex (Montoya-Lerma et al. 2011). An. goeldii , thought to be restricted to the Amazon Basin (Calado et al. 2008), may extend further north in sympatry with An. nuneztovari (Scarpassa & Conn 2011). The close genetic association between An. nuneztovari , an important malaria vector in Colombia and Venezuela (Bourke et al. 2010), and An. goeldii , raises concerns over the putative vector status of the latter and its habitat specificity. Thus, studies involving the population structure and ecological distribution of An. goeldii are useful in determining the extent of its distribution and divergence.

An. triannulatus (Neiva & Pinto) has been reported from Nicaragua to Argentina (Faran 1980). Originally described as a polymorphic taxon (Faran & Linthicum 1981), the intraspecific variation was attributed to adaptation to different habitats (Forattini 1962, Faran 1980, Rosa-Freitas et al. 1998). More recently, it has been recognised as a complex of at least three species (Silva-do-Nascimento & Lourenço-de-Oliveira 2002, Silva-do-Nascimento et al. 2006, 2011). Recent studies showing anthropophilic feeding and biting activity both indoors and out (Brochero et al. 2006, Gutiérrez et al. 2009), support vector incrimination in some areas of Latin America (Silva-do-Nascimento et al. 2006, de Barros et al. 2007).

As factors affecting speciation and population differentiation evolve at different rates (de Queiroz 2007), local and regional analysis of ecology and population structure may identify behavioural or ecological differences that complement divergence studies. Use of nuclear markers in addition to mitochondrial markers reduces the process error associated with coalescence (Zink & Barrowclough 2008) and has led to an informal set of genes sequenced frequently for phylogenetic studies including the mitochondrial cytochrome oxidase I ( COI ) (Avise 2000, Molina-Cruz et al. 2004, Sunghoon et al. 2011), nuclear white gene (Besansky & Fahey 1997 , Reidenbach et al. 2009) and internal transcribed spacer (ITS)2 gene family (Zapata et al. 2007, Shultz & Wolf 2009, Wiemers et al. 2009). Within the mitochondrial genome the more conservative “folmer region” has become the standard barcode fragment with a 3% divergence threshold for speciation events (Hebert et al. 2003).

The diversity of Neotropical anophelines has been attributed to the ability of dipterans to readily adapt to and utilise a broad variety of ecological niches (Grimaldi & Engel 2005) with more generalised populations being more ecologically heterogeneous (Bolnick et al. 2007). Thus, the genetic heterogeneity seen in the Triannulatus and Nuneztovari complexes may be explained, in part, by environmental heterogeneity (Gram & Sork 2001), including local or microgeographic adaptation, past and present ecological barriers and/or demographic events (Brouat et al. 2004). Additionally, information on the dynamics of larval anophelines is limited (Shililu et al. 2003) and may provide insights into population differences, thereby complementing adult studies.

The purpose of this study was to examine the distribution, abundance and genetics of largely generalist and sympatric species from 76 larval sites to address the following questions: (i) what is the level of habitat differentiation between species, (ii) do patterns of diversity in An. goeldii and An. triannulatus indicate a common cause and (iii) can genetic differentiation be explained by demographic phenomena and (or) geographic boundaries?

MATERIALS AND METHODS

Mosquito collection - Fourth instar anopheline larvae were collected from 13 localities, including eight in the state of Roraima (RR) in northern Amazonian Brazil and five in the state of Pará (PA) in eastern Amazonian Brazil (Fig. 1, Supplementary data). Sampling was done once for each randomly selected larval site. If the site was positive for anophelines, collections continued for an hour to measure relative species abundance. The study was conducted over three years (2009-2011) during malaria transmission season in each state. The larval sites in RR were located between Amajari and Ecuador, roughly 500 km along the RR-174, in 2009 and 2011, and in PA between Mojú and Marabá, approximately 560 km along PA-475/263/150 in 2010. The collection protocol was approved by the Evandro Chagas Institute, Pará State Ethical Committee and the New York State Department of Health, Institutional Review Board. Physical characteristics of the larval habitats were recorded. These included biotic factors such as temperature, alkalinity, pH, conductivity, salinity and turbidity. Additionally, canopy coverage was estimated by spherical densitometer (Minakawa et al. 2005).

Fig. 1

map of Brazil indicating the location of the 13 localities and species distributions. See Supplementary data for more information. PA: state of Pará; RR: state of Roraima.

Larvae were identified morphologically as An. triannulatus s.l. and An. nuneztovari using the key of Deane et al. (1946) and confirmed by species specific ITS2 restriction fragment length polymorphism (Zapata et al. 2007). Morphological keys do not reflect the recent differences between the Colombian and Venezuelan An. nuneztovari and the Brazilian An. goeldii , samples, therefore samples were sequenced and compared to those from Calado et al. (2008) using a Bayesian inference (BI) phylogenetic approach (data not shown). Total genomic DNA was extracted using the DNeasy Blue Tissue Kit (Qiagen, CA, USA) and maintained at Griffin Laboratory at -80°C.

Statistical analysis - A paired sample t test was used to compare An. triannulatus and An. goeldii larval occurrence across all sites. Subsequent analysis of similar larval sites based upon the habitat type and type of water body sampled (pond, pool, swamp etc.), examined variation in larval densities using two-way analysis of variance (ANOVA). Additional analyses examining differences in habitat ecology, based on biotic factors, surrounding vegetation and canopy coverage were examined to test the hypothesis that environmental heterogeneity is positively correlated with genetic diversity.

Amplification and sequencing - A representative sample of one-50 mosquitoes per locality was selected for DNA amplification and sequencing of the COI . Subsets of up to 10 samples per locality for the white gene and up to five for ITS2 were amplified and sequenced. The 658-bp barcode fragment of the COI gene was amplified using the forward primer LCO1490 and the reverse primer HCO2198 (Folmer et al. 1994). Individual polymerase chain reaction (PCR) reactions were performed using Ready-To-Go RT-PCR Beads (Amersham Pharmacia/Biotech, NJ, USA) or PCR Supermix High Fidelity (GoTaq, Promega, CA) run on a Bio-Rad PTC-100 or 200 series thermal cycler (Bio-Rad Inc), using the conditions stipulated in Hebert et al. (2003) . The Applied Genomics Technology Core (Wadsworth Center) carried out the sequencing on an ABI PRISM 3700 automated DNA sequencer. The forward and reverse COI sequences were aligned using Sequencher 3.0 (Gene Codes Corps, MI, USA), grouped together by sight and trimmed in MEGA v. 3.1 (Kumar et al. 2004) creating a 600-bp fragment. In addition, a subset of three-five An. triannulatus samples from each locality was sequenced following the protocol of Pedro and Sallum (2009), to yield a 449-bp fragment of the 3' COI . High levels of intraspecific variation, deep structure and complex phylogenies based on the 3' end of the COI gene (Moreno et al. 2013) are not congruent with recently described species An. triannulatus , Anopheles haplophylus and An. triannulatus “C” (Silva-do-Nascimento & Lourenço-de-Oliveira 2002, Silva-do-Nascimento et al. 2006, 2011). The barcode region is more variable than the 3' end of the COI (Ruiz-Lopez et al. 2012) and able to improve resolution at deeper nodes (Saunders & Le Gall 2010). Therefore, analyses herein focused on the barcode fragment and included only a subset based on the 3' end to allow systematic comparisons of An. triannulatus populations across Brazil.

Nuclear genes were included to investigate the history of the nuclear genome and provide independent evidence of population divergence (Loiaza et al. 2012). Additionally, nuclear genes do not suffer the short comings of mitochondrial markers (haploid, small genome size, short coalescent time, maternal inheritance) ( Ballard & Whitlock 2004) and are not as biased in base composition (Brower & De Salle 1998). An 800-bp fragment of the white gene was amplified using the W2R and WF primers, with the PCR conditions as reported in Mirabello and Conn (2008) . A 500-bp fragment of the ribosomal ITS2 was amplified using the primers 18S and 28S following the parameters in Li and Wilkerson (2005) . The PCR products were cleaned, sequenced and aligned creating a 570-bp fragment of the white gene including the intron. The final product size for ITS2 was 433 bp and 371 bp for An. triannulatus and An. goeldii , respectively. A single sequence per individual/gene was included in all subsequent analyses. Haplotypes of heterozygous individuals were inferred for each set of gene sequences using PHASE 2.1 software (Stephens et al. 2001, Stephens & Donnelly 2003, White et al. 2009). Unique sequences for all markers are available in GenBank (accessions KC167670-KC167816).

Genealogical relationship - Genealogical trees for each species were estimated for the concatenated ( COI + white ) data set using BI analysis performed with MrBayes v. 3.1 (Huelsenbeck & Ronquist 2001, Ronquist & Huelsenbeck 2003). Data were partitioned by gene using the model of nucleotide substitution (HKY+I and TrN+G for An. triannulatus and TrN+I+G and TPM2uf+G for An. goeldii ) that best fit the white gene and COI , respectively, determined by jModelTest (Posada 2008). ITS2 sequences were excluded from these analyses because of the limited sample size. The settings were two simultaneous, independent runs of the Markov Chain Monte Carlo for two million generations, sampling every 1,000 generations with a “burnin” of 25%. The outgroup Anophe-les marajoara (DQ076221 and AY956296 deposited in GenBank) was chosen based upon its phylogenetic position in Sallum et al. (2000) .

Genetic variation - Genetic structure of multiple populations based on statistical parsimony groupings were examined by analysis of molecular variance (AMOVA) using Arlequin v. 3.1.1 (Excoffier et al. 2005) and spatial analysis of molecular variance (SAMOVA), v. 1.0 ( Dupanloup et al. 2002). Both algorithms were used to cluster the 600-bp barcode fragment data into genetically and in the case of SAMOVA, geographically, homogeneous populations by generating F statistics ( F , F , F ). Each species was analysed separately, An. triannulatus into groups of K = 2-16 and An. goeldii into groups of K = 2-20, with 1,000 simulated annealing steps from each of 100 sets of initial starting conditions. Groups were based on individual larval sites. A second analysis was undertaken for An. triannulatus using the 449-bp fragment haplotypes of Pedro and Sallum (2009) (GU445849-899 and AF417702) combined with the samples in the present study, to test their hypothesis of geographical sub-structuring within the species based on five populations across Brazil.

Population structure and demographic history - A statistical parsimony network (SPN) estimated genealogical relationships among haplotypes with a 95% identity for each marker using TCS 1.13 (Clement et al. 2000), except An. goeldii ITS2 sequences, which were examined at 94% sequence similarity. Homoplasy in all networks was resolved using the algorithm estimation rules in Crandall and Templeton (1993) . Estimates of time to coalescence were calculated for the COI fragment only and compared using θ values (Watterson 1975).

The differentiation and polymorphism statistics for COI and white gene sequences by species and locality were computed in DnaSP, v. 5.0 (Librado & Rozas 2009). Because of the high proportion of shared haplotypes among habitats within a single locality, samples from an individual town/city were pooled to more accurately represent diversity statistics. Furthermore, haplotype distributions across four ecoregions were examined for each of the seven environmental variables to test the hypothesis that climatic differences sustain (and drive) diversity (Hoorn et al. 2010). To test this hypothesis, we considered both the availability of breeding sites and haplotype diversity as measures of gene flow. Individual sites were grouped within ecoregions, based on descriptions from Rosa-Freitas et al. (2007) .

The hypothesis of strict neutrality for each species was examined using DnaSP, with the statistics D (Tajima 1989) and D and F (Fu 1993) evaluating the likelihood of background selection and Fu's F (Fu 1997) and R (Ramos-Onsins & Rozas 2002) examining evidence of population expansion. The mismatch distributions (simulated in Arlequin) were conducted to confirm demographic expansion using sudden and spatial models of expansion. Statistically significant differences between observed and simulated distributions were evaluated with the sum of square deviations to reject the hypothesis of demographic expansion (Zarza et al. 2008). Dates of population expansion were estimated for the COI gene with the formula T = τ/2 u (Rogers & Harpending 1992) using 10 generations per year and the mutation rate of Drosophila estimated at 1 × 10−8 substitutions per site per year (Walton et al. 2000, Yang et al. 2011).

RESULTS

Site diversity and larval abundance - Larvae were predominantly collected from the edges of naturally occurring ponds, comprising 49.19% (n = 200) of the total Anopheles identified. An. goeldii and An. triannulatus larvae were relatively uncommon in collections from seepages, stream margins and swamps, where they comprised only 10.9% of the total. Both species were collected from the same site 17 times (Supplementary data). Twelve sites contained An. triannulatus and not An. goeldii , whereas 13 sites contained An. goeldii , but not An. triannulatus . A large proportion, 68.97% (n = 140), of An. triannulatus larvae were sampled from ponds. An. goeldii larvae appeared equally likely to inhabit artificial pools (37.93%, n = 77), ponds (29.55%, n = 60) and temporary ditches (25.12%, n = 51). No correlations between abundance and habitat type for either species were detected by ANOVA (data not shown).

Habitat characteristics of An. triannulatus and An. goeldii were similar: the 76 sites were not differentiated by larval abundance ( t = -0.2980, df = 42, p = 0.3836). ANOVA detected a significant difference in turbidity between sites containing only An. goeldii and those negative for either species ( F = 4.89, p = 0.032) that was not supported by a t test ( t = 1.179, df = 16, p = 0.055). There are no other significant differences between environmental variables for sites unique for either species and sites either positive or negative for both (Table I ).

TABLE I

Characteristics (± 95% confidence intervals) of habitats with Anopheles triannulatus only, Anopheles goeldii only, both species or neither species present

	An. triannulatus	An. goeldii	Both	Neither
Collections (n)	12	13	17	34
pH	5.5 ± 0.00	5.58 ± 0.15	5.68 ± 0.24	5.62 ± 0.16
Alkalinity (ppm)	35 ± 29.40	33.85 ± 27.50	47.06 ± 27.91	52.35 ± 42.58
Temperature (°C)	26.29 ± 1.37	27.44 ± 1.88	26.79 ± 1.10	26.48 ± 0.95
Conductivity (μS/cm)	35.68 ± 16.33	39.49 ± 19.68	34.01 ± 12.93	26.04 ± 9.76
Salinity (ppm)	18.09 ± 8.13	19.92 ± 10.31	25.11 ± 17.34	13.24 ± 4.90
Turbidity (JTUs)	21.67 ± 15.93	28.46 ± 16.88	17.50 ± 6.88	13.09 ± 5.36
Canopy coverage (%)	11.28 ± 14.20	13.74 ± 13.43	13.77 ± 11.38	10.55 ± 6.33

Genealogical relationship - The concatenated ( COI + white gene with intron) BI tree for An. triannulatus , with 50 parsimoniously informative sites, indicated two strongly supported (97 and 98 posterior probabilities) distinct populations of An. triannulatus (Fig. 2 ). The two populations correspond to the Brazilian states, with the exception of samples from Mojú and Goianésia (9 and 10, respectively), that were included in the northern group (RR localities). The close relationship between the RR samples and localities Mojú and Goianésia was also evident in the SPN. This may indicate a historical boundary separating north and south populations, with the northern limit occurring just south of the Amazon River and including parts of Goianésia and Mojú. The more derived population 1 corresponds to localities in the Amazon Basin, near the hypothesised population origin (Pedro & Sallum 2009), with the majority of population 2 occurring along one of the proposed routes of lineage expansion. Within grouping 2, smaller subdivisions, related to geography and corresponding to populations 2 and 3, were apparent between haplotypes 9 and 3 that are separated by six mutations and two haplotypes 4 and 1, separated by two mutations (Fig. 3).

Fig. 2

concatenated (cytochrome oxidase I + white gene) Bayesian inference trees, partitioned by gene using the model of nucleotide substitution that best fit the data as determined by jModelTest. blue: samples from Goianésia, Tucuruí, Jacundá and Marabá [state of Pará (PA)] corresponding to An. triannulatus population 3; green: Mojú and Goianésia samples (PA) in addition to An. triannulatus population 2; orange: localities in the state of Roraima and An. triannulatus population 1. Asterisk means Anopheles marajoara outgroups.

Fig. 3

statistical parsimony network of cytochrome oxidase I sequences at 95% support. black: localities in the state of Roraima; grey: Mojú and Goianésia samples from the state of Pará (PA); n: number of sequences; white: samples from Goianésia, Tucuruí, Jacundá and Marabá (PA). For An. triannulatus , black refers to population 1, grey (with the addition of haplotype 14) to population 2 and white to population 3.

The concatenated ( COI + white gene with intron) BI tree for An. goeldii , with 32 parsimoniously informative sites, indicated two moderately supported (72 and 68 posterior probabilities) populations which correspond to the Brazilian states (Fig. 2). Unlike An. triannulatus , there does not appear to be any sub-structuring in the samples from PA or support for a separate Goianésia/Mojú population.

Genetic variation - AMOVA indicated consistently higher F and F fixation indices across genes and various COI groupings for An. triannulatus , i.e., genetic structuring among and within localities rather than between groups (Table II). However, there was strong genetic differentiation (41.01%) also seen in the COI grouping between states, which correlates well with the groups defined by the Bayesian tree.

TABLE II

Analysis of molecular variance calculations for inter and intravariation among species and populations for various markers

		Anopheles triannulatus		Anopheles goeldii
Gene	Source of variation	Variance (%)	Fixation index	Variance (%)	Fixation index
COI	Between RR and PA	41.01	FCT = 0.41a	45.12	FCT = 0.45b
	Among localities within states	29.23	FSC = 0.50b	5.29	FSC = 0.10b
	Within a given locality	29.76	FST = 0.70b	49.59	FST = 0.50b
white	Between RR and PA	-3.01	FCT = -0.03c	89.46	FCT = 0.89a
	Among localities within states	58.62	FSC = 0.57b	0.42	FSC = 0.04c
	Within a given locality	44.38	FST = 0.56b	10.12	FST = 0.90b
COI	Between SPN populations	38.85	FCT = 0.39c	-	-
	Among localities	33.92	FSC = 0.55b	-	-
	Within a given locality	27.23	FST = 0.73b	-	-

a : p < 0.05; b : p < 0.001; c : not significant; COI : cytochrome oxidase I; PA: state of Pará; RR: state of Roraima; SPN: statistical parsimony network derived groupings.

The small genetic differentiation of An. goeldii between localities suggests panmixia with no significant barriers to gene flow within individual states. AMOVA indicated a large proportion of variation between the two states and among samples from different breeding site types in a given locality ( Table II). Genetic structuring in the latter may indicate environmental differences and niche partitioning as drivers of diversity.

SAMOVA analysis for both species found no evidence of a geographic barrier. The addition of the An. triannulatus 449-bp fragment of the COI gene (n = 35) to all samples from Pedro and Sallum (2009) was unable to replicate the five groupings originally described across Brazil.

Population structure and demographic history - Analysis of mitochondrial data from 129 An. triannulatus and 133 An. goeldii specimens (Fig. 3) led to the discovery of two very different demographic histories between these sympatric species, which often co-occur within the same larval habitats. The time of the beginning of lineage divergence was calculated following Da/2 k , where Da is the average nucleotide divergence between lineages and 2 k is the divergence rate (Heckel et al. 2005). We utilised the mutation rate of Drosophila (Walton et al. 2000, Yang et al. 2011) to estimate these parameters and their significance based on the coalescence (DnaSP v. 5.10) (Librado & Rozas 2009) and MEGA v. 5.1 (Tamura et al. 2011).

A parsimony network for An. triannulatus indicated that a minimum of six mutational steps separate three potential populations based on geography (Fig. 3). The most common haplotypes were 1 (n = 18), 2 (n = 15), 3 (n = 11) and 6 (n = 9) (Fig. 3), with seven haplotypes (15.6%) shared among localities and 38 (84.4%) unique (Supplementary data). Overall, An. triannulatus was consistent with a category I phylogeographic pattern (Avise 2000) characterised by pronounced genetic gaps between spatially circumscribed haplogroups. A high proportion of singletons (23, 51.1%), in addition to individual populations within An. triannulatus containing star-shaped nodes (surrounding haplotypes 2, 3 and 5), supports demographic expansions, background selection or selective sweeps (Slatkin & Hudson 1991, Fu 1997). Based upon the 600-bp fragment, D between the northern and southern states is 0.0118 [standard deviation (SD) = 0.003]. Therefore, the estimated divergence is 437,500-737,500 years ago, in the Pleistocene. The two PA populations had a D = 0.014 (SD = 0.004), indicating more ancestral divergence (congruent with Bayesian tree), approximately 500,000-900,000 years ago.

An. goeldii COI is consistent with a category III phylogeographic pattern ( Avise 2000) where most haplotypes are closely related, yet some geographic localisation exists. The most common haplotypes were 46 (n = 21), 47 (n = 33) and 48 (n = 13) (Fig. 3), with 14 haplotypes (41.2%) shared among localities, and 20 (58.8%) unique (Supplementary data). Similar to An. triannulatus , there was a relatively high proportion of singletons (14, 41.2%), indicating a demographic expansion, background selection or selective sweep (Slatkin & Hudson 1991, Fu 1997), particularly in PA (singletons arising from haplotype 56). In contrast, haplotypes from RR exhibited longer branches, missing haplotypes and a near equal distribution of shared haplotypes and single mutations, consistent with a signal of an older lineage with balancing selection (Fig. 3). D between northern and southern states is 0.008 (SD = 0.002). Therefore, the estimated divergence is 316,000-516,500 years ago (Pleistocene).

Networks of An. triannulatus nuclear data ( white gene and ITS2) indicate a history of contractions and expansions. Similar to the COI gene, the ITS2 network depicts a category I pattern whereas the white gene reveals a much more homogeneous haplotype distribution. Overall, the white gene follows a category II intraspecific pattern, with pronounced gaps between some branches and some clustering of anciently separated geographic populations (Fig. 4).

Fig. 4

statistical parsimony network of white gene (A) and internal transcribed spacer 2 sequences (B) at 95%, except ITS2 An. goeldii which is at 94% sequence support. black: localities in the state of Roraima; grey: Mojú and Goianésia samples from the state of Pará (PA); n: number of sequences; white: samples from Goianésia, Tucuruí, Jacundá and Marabá (PA).

Networks of An. goeldii nuclear data ( white gene and ITS2) reveal distinct geographical populations. The white gene has a distinct category I pattern separated by five mutational steps, while the ITS2 haplotype clusters (Fig. 4) are less pronounced, but still correspond primarily to RR and PA with a few outliers. Polymorphism analyses of the COI sequences found haplotype and nucleotide diversity was greater among An. triannulatus populations (Supplementary data) in contrast to An. goeldii .

Genetic polymorphism analyses of the COI sequences indicated that haplotype and nucleotide diversity was greater among An. triannulatus populations (Supplementary data) in contrast to An. goeldii . The white gene exhibited lower nucleotide diversities compared to COI , but surprisingly An. triannulatus continued to demonstrate high haplotype diversity (Supplementary data), which combined with relatively low nucleotide diversity may signal a rapid expansion ( Avise 2000, de Jong et al. 2011). High COI variation between states and geographic populations is strongly supported by population differentiation statistics (Table III), although the G values were low, suggesting a history of gene flow (Table III).

TABLE III

Inter and intrapopulation differentiation of Anopheles triannulatus and Anopheles goeldii populations based on statistical parsimony network

	Comparison
	An. triannulatus	An. goeldii
	RR vs. PA	Three populations	RR vs. PA
HS	0.89a	0.86a	0.86a
K S *	1.77a	1.41a	1.60a
Z*	7.43a	7.00a	7.83a
Snn	0.99a	1.00a	0.97a
χ2	129a	258a	125.84a
GST	0.06	0.09	0.045
Kt	10.09	10.09	7.43

a : p < 0.0001; G ST : genetic differentiation; H S : genetic differentiation based on haplotype data; K S *: differentiation based on sequence data; K t : average number of nucleotide differences; n: number of total individuals compared; PA: state of Pará; RR: state of Roraima; S nn : measures how often the nearest neighbours of sequences are found in the same population; χ 2 : genetic differentiation based on allele frequencies; Z *: rank statistic to analyse sequence similarity.

In general, neutrality tests were negative for both species (Table IV), suggesting an excess of rare polymorphisms consistent with either negative/positive selection or an increase in population size, although only one of these tests ( F ) was significant for An. triannulatus, i.e., populations 2 and 3. Mismatch distribution for sudden expansions did not support this finding and, instead, revealed marginally significant expansion patterns for only population 1, corresponding to RR. Sudden and spatial populations for An. goeldii could not be rejected (Fig. 5). Estimated expansions for both species are during the Pleistocene, An. triannulatus 129,693 years ago (59,857-194,877) and An. goeldii 112,387 years ago (17,111-187,449).

TABLE IV

Multiple statistical analyses evaluating the probability for a recent population expansion

Species	n	D T	D a	F a	F S	R 2
Anopheles triannulatus	129	-0.22187	-1.28038	-0.99732	-10.329a	0.08650
Roraima	57	-0.37846	-1.32424	-1.17022	-1.488	0.09200
Pará	72	-0.47543	-0.97094	-0.93438	-5.842	0.08900
Population 1	56	-0.08718	-1.11433	-0.88956	-1.209	0.10200
Population 2	36	-1.06921	-0.55581	-0.85324	-3.582a	0.07770
Population 3	37	-1.02686	-0.40080	-0.71980	-4.809a	0.07800
Anopheles goeldii	133	-0.30013	-1.24177	-1.01960	-5.867	0.08260
Roraima	103	0.06971	-1.6295	-1.1537	-1.677	0.0989
Pará	30	-0.0415	0.15645	0.11023	-3.321	0.1185

a : significance, p < 0.50; D : Fu and Li's D (Fu 1993); D T : Tajima's D (Tajima 1989); F : Fu and Li's F (Fu 1993); F S : Fu's Fs (Fu 1997); n: sample size; R 2 : Ramos-Onsin's and Roza's R 2 (Ramos-Onsins & Rozas 2002).

Fig. 5

representative mismatch distributions with calculated τ and estimated time of sudden expansion. Bars indicate observed values. Line represents the model for all Anopheles triannulatus (A), all Anopheles goeldii (B), An. triannulatus population 2 (C) (Mojú and Goianésia samples) and An. triannulatus population 3 (D) (Goianésia, Tucuruí and Jacundá samples). ya: years ago.

DISCUSSION

Distribution patterns of mosquito species are dependent on the availability of suitable aquatic habitats (Grillet 2000). Overlapping geographic distributions, in addition to the absence of significant difference in species abundance across sampled habitats, suggest that An. goeldii is a suitable comparison species for An. triannulatus . The lack of any significant differentiation between sympatric habitats, to those unique to a species or absent of either, suggests these species may be more broadly adapted (Eisenberg et al. 2000, Chase & Knight 2003). However, co-occurrence of species does not imply the same ability to exploit shared habitats (Diabate et al. 2005).

The combined sequence results support distinct demographic histories for the two species despite similar distribution and larval habitat co-occurrence. Across molecular markers An. goeldii consistently depicts two moderately supported populations corresponding to the two Brazilian states, however, An. triannulatus appears to have undergone historical fragmentation, subsequent secondary contact and homogenisation followed by more recent divergence due to allopatry. Either environmental changes (e.g., habitat availability and climatic oscillations) alone cannot account for population divergence or differences in biology allow these two species to respond differently to climatic oscillations.

The clear separation of populations based on Brazilian states in An. goeldii strongly supports the alternative hypothesis of geographic isolation driving divergence rather than environmental changes. While often contested in the neotropics, the Pleistocene refugia hypothesis (Haffer 1969) postulates that climatic fluctuations resulted in the fragmentation of continuous forest into refuges separated by expanses of grass-dominated savannah or desert. Alternatively, the Miocene marine incursion hypothesis suggests sea level changes during the late Tertiary (23.7-1.8 million years ago) created three isolated regions: the Brazilian Shield, the Guiana Shield and the eastern Andean slopes (Webb 1995, Bates 2001, Hoorn et al ). The last commonly explored mechanism of divergence is the riverine barrier hypothesis, where major rivers act as barriers to gene flow, thus promoting speciation between populations on either side (Wallace 1852, Aleixo 2004, Mirabello & Conn 2008).

Although no boundary was indicated by SAMOVA, this does not discount the possibility that an historical boundary existed, separating north and south populations. Although one study found evidence for waterways over 4 km wide being complete barriers to gene flow (Pedro & Sallum 2009), others have determined that the Amazon is not an effective barrier (Fairley et al. 2002, Mirabello & Conn 2008) and that geographic distance and/or differing demographic histories may be the main forces responsible for partitioning genetic variation (Mirabello & Conn 2008).

The discrepancy regarding the inclusion of Mojú and some of Goianésia among An. triannulatus RR population, but their inclusion with geographically similar PA samples in An. goeldii may be attributed to individual species' responses to dispersal opportunities (Brouat et al. 2004). The two An. triannulatus populations identified in Goianésia could also be the result of a suture zone. Species dispersal activity and differential selection can cause hybrid zones to appear mobile and settle as boundaries between environments or in regions of low population density (Barton & Hewitt 1985, Mallet 1993, 2010). Suture zones that span several species suggest divergence was driven by a common influence of environmental change (Morgan et al. 2010). Further evaluation of species in the Goianésia area could resolve the contribution of environmental change to population diversity.

Overlapping geographic distributions, in addition to the absence of significant difference in species abundance across sampled habitats, implies little habitat differentiation between species and a broad ecological range. Although biological differences may drive divergence uniquely in each of these sympatric species, a common cause for diversity patterns cannot be ruled out. Geographic boundaries/distance, variation in dispersal activity and ecological range may each contribute to the distinct patterns of divergence seen in An. triannulatus and An. goeldii . As species complexes are common in anophelines, additional studies are needed to resolve the extent of divergence. Continued studies using multiple markers at a local and continental scale will resolve the status of speciation and may clarify what role, if any, a given population, lineage or species has on malaria transmission, because population differences influence peak biting times, feeding preferences and vectorial capacity (Lounibos & Conn 2000).

TABLE

Species distribution of anopheline larvae in different larval habitats by locality over three transmission seasons (2009-2011)

Locality	Habitat type	Anopheles triannulatus (n)	Annealing (°C)
Amajarí (1)	Marsh	0	0
	Pool	1	1
Alto Alegre (2)	Ditch	0	0
	Marsh	4	0
	Pond	0	7
	Seepage	0	0
	Stream margin	2	1
Boa Vista (3)	Ditch	0	0
	Marsh	0	0
	Pond	0	4
	Pool	0	14
	Seepage	0	6
	Stream margin	0	0
Mucajaí (4)	Ditch	0	0
	Pool	0	62
Iracema (5)	Ditch	10	45
	Pond	40	11
Petrolina (6)	Ditch	0	0
	Pond	0	0
	Stream margin	0	3
Martins Pereira (7)	Stream margin	1	2
	Swamp	14	4
Ecuador (8)	Pond	0	2
	Seepage	0	5
	Stream margin	0	0
	Swamp	0	0
Mojú (9)	Pond	19	19
	Marsh	66	31
	Stream margin	0	0
Goianésia (10)	Marsh	0	0
	Pond	14	0
	Pool	0	0
	Stream margin	0	0
	Swamp	0	0
Tucuruí (11)	Pond	7	2
	Stream margin	6	0
Jacundá (12)	Marsh	5	0
	Pond	11	0
	Stream margin	0	2
Marabá (13)	Ditch	0	6
	Marsh	0	0
	Pond	2	3
	Stream margin	1	0
Total		203	230

1-13: see Fig. 1 for more information.

TABLE

Description of shared cytochrome oxidase I haplotypes and genetic polymorphism statistics

	Site	n	Haplotypes	π (SD)	1-∑ƒ i 2 (SD)	P	K
An. triannulatus	1	1	1	-	-	-	-
	2	4	1 (2) , 40, 41	0.01 (0.00)	0.83 (0.22)	12	7.67
	5	37	1 (15) , 2 (6) , 7, 8, 9 (2), 35, 36 (2), 37, 39 (4), 43, 45	0.01 (0.00)	0.80 (0.05)	20	5.61
	7	15	2 (9) , 14, 38, 42 (3), 44	0.00 (0.00)	0.63 (0.13)	14	2.27
	9	30	3 (11), 4 (7), 10, 11, 30, 31 (2), 32, 33, 34 (4)	0.00 (0.00)	0.81 (0.05)	12	1.95
	10	14	6 (3) , 12, 13 (3), 15, 16 , 17 (2), 24 , 25 , 28	0.01 (0.00)	0.92 (0.05)	20	6.99
	11	13	5 (2) , 6 , 16 , 18 , 19 , 20 , 21 , 22 , 26 , 27	0.01 (0.00)	0.96 (0.04)	16	4
	12	15	5 (4) , 6 (4) , 23 (2), 24 , 25 (2 ), 29	0.01 (0.00)	0.83 (0.06)	9	2.99
		129		0.02 (0.00)	0.95 (0.01)	57	10.1
An. goeldii	2	4	46 , 47 (3)	0.00 (0.00)	0.50 (0.27)	2	1
	3	15	46 (4) , 47 (7) , 48 , 49 , 77, 78	0.01 (0.00)	0.74 (0.09)	15	3.26
	4	18	46 (3) , 47 (6 ), 48 (2) , 49 (2) , 50 , 63 , 67, 72	0.01 (0.00)	0.86 (0.06)	18	6.78
	5	50	46 (7) , 47 (14) , 48 (10) , 49 (3) , 50 (2) , 51, 52, 63 (3) , 70, 71, 73 (2) , 74 (2), 75	0.01 (0.01)	0.86 (0.03)	24	6.29
	6	3	47 (2) , 72	0.01 (0.01)	0.67 (0.31)	11	7.33
	7	6	46 , 50 , 68, 69 , 73 , 76	0.01 (0.00)	1.00 (0.10)	18	8.47
	8	7	46 (5) , 47 , 53	0.01 (0.00)	0.52 (0.21)	12	3.81
	9	17	53 , 56 , 57, 59, 60 (2) , 61 (3) , 62 (2), 64 , 65, 66	0.00 (0.00)	0.93 (0.04)	10	2.18
	11	2	54 , 63	-	-	-	-
	13	11	54 (3) , 55 (2), 58, 60 , 63 (2) , 64	0.01 (0.00)	0.91 (0.07)	13	5.6
		133		0.01 (0.01)	0.90 (0.002)	44	7.43

K : average number of differences; n: number of sequences; P : polymorphic sites; SD: standard deviation; 1-Σf i 2 : haplotype diversity; Π: nucleotide diversity; 1-13: see Fig. 1 for more information. Bold indicates shared haplotypes and indicates diversity among lineages or species.

TABLE

Description of shared white haplotypes and genetic polymorphism statistics for Anopheles triannulatus and Anopheles goeldii

Species	Site	n	Haplotypes	π (SD)	1-∑ƒ i 2 (SD)	P	K
An. triannulatus	1	1	B	-	-	-	-
	2	3	B , E (2)	-	-	-	-
	5	10	B (3), C (4), D (2), E	0.00 (0.00)	0.51 (0.16)	7	2.47
	7	4	E (2), P, T	-	-	-	-
	9	9	A (8), U	0.00 (0.00)	0.22 (0.17)	1	0.22
	10	9	E , N (2), O, Q (3), R, S	0.01 (0.00)	0.81(0.12)	9	3.94
	11	10	E , F, I , J (2), K, L (2), N (2)	0.01 (0.00)	0.93 (0.06)	11	4.02
	12	8	G (2), H, I , J , L (2), M	0.00 (0.00)	0.93(0.08)	4	1.54
		54		0.01 (0.00)	0.93 (0.02)	24	5.43
An. goeldii	2	4	V , W (2), CC	-	-	-	-
	3	10	V (3), W (6), X	0.00 (0.00)	0.600(0.13)	4	1.07
	4	10	V , W (9)	0.00 (0.00)	0.20 (0.15)	1	0.20
	5	9	W (9)	0.00 (0.00)	0.000(0.00)	0	0
	8	1	W	-	-	-	-
	9	10	Y (9), Z	0.00 (0.00)	0.200(0.15)	2	0.4
	13	8	Y (5), AA (2), BB	0.00 (0.00)	0.61 (0.16)	3	1.04
		52		0.01 (0.00)	0.66 (0.05)	12	2.85

K : average number of differences; n: number of sequences; P : polymorphic sites; SD: standard deviation; 1-Σf i 2 : haplotype diversity; Π: nucleotide diversity; 1-13: see Fig. 1 for more information. Bold indicates diversity among lineages or species.