Anopheles (Nyssorhynchus) marajoara Galvão & Damasceno, a member of the Albitarsis Complex, was recently recognized as the primary vector of malaria parasites in northeastern Amazonia, Brazil (Conn et al. 2002). Other species in the complex are An. albitarsis Lynch-Arribálzaga, An. deaneorum Rosa-Freitas, and an unnamed species “B” (Kreutzer et al. 1976; Linthicum 1988; Narang et al. 1993; Rosa-Freitas and Deane 1989; Wilkerson et al. 1995a, b). From studies carried out in Rondônia State, Brazil, there is also evidence to support the importance of An. deaneorum as an important malaria vector (Klein et al. 1991a, b). Wilkerson et al. (1995a, b) separated the four largely isomorphic species by using species-specific random amplified polymorphic DNA (RAPD-PCR) markers. Their analysis relied on an empirical assumption of multiple correlated “fixed” markers to hypothesize reproductive isolation, similar to the use of correlated morphological characters for the same ends, but they did not address confounding factors inherent in RAPD markers such as possible nonhomology of comigrating bands or linkage of markers. Lehr et al. (2005), based on complete mitochondrial DNA cytochrome oxidase I (mtDNA COI) sequence, question the validity of this approach, suggest a fifth species, and present a COI gene tree for all five. Using specimens from the Wilkerson et al. (1995a, b) RAPD-PCR studies, Merritt et al. (2005) analyzed a portion of the white gene that contains the white fourth intron. They found that the intron was present in An. marajoara but not in the other three species (also noted by Krzywinski et al. (2001) for An. albitarsis). Phylogenetic analysis of coding sequence in the area of the fourth intron, correlated with intron loss hypotheses, resulted in strong support for a single intron loss event in this species complex and gave a topology different from that found both by Lehr et al. (2005) and from that presented here.

Anopheles subgenus Nyssorhynchus includes 33 species (Harbach 2004), including the two most important vectors in the New World tropics: An. darlingi Root (Linthicum 1988) and An. albimanus Wiedemann (Faran 1980). The subgenus is divided into three sections based on morphological characters: the Argyritarsis Section, which includes the Albitarsis Complex and An. darlingi; the Albimanus Section, which includes An. albimanus; and the little known Myzorhynchella Section (Peyton et al. 1992). Relationships among the 33 included species are not well resolved, but it is known that the Argyritarsis and Albimanus Sections are paraphyletic relative to each other because one putative clade contains both An. darlingi and An. albimanus, suggesting a possible evolutionary link to vector capacity (Conn 1998, Sallum et al. 2000). Cryptic species are common in An. (Nyssorhynchus) and in Anopheles in general, and most groups that are closely studied yield new taxa (Rosa-Freitas et al. 1998), with An. darlingi being an apparent exception (Manguin et al. 1999). Considering the medical importance of this complex and our general lack of knowledge regarding the relationships within subgenus Nyssorhynchus, we undertook this study to corroborate results produced by RAPDs and to investigate how the Albitarsis Complex species are related to each other. We report here a molecular phylogenetic analysis of the four species, initially separated by RAPDs, by using two ribosomal DNA (rDNA) sequences, internal transcribed spacer two (ITS2) and the D2 expansion of the 28S subunit (D2), and partial sequence from two mitochondrial genes, NADH dehydrogenase 4 (ND4) and COI.

Table 1 List of species and source of specimens used in this study

Table 2 Sequences of COI, ND4, D2, and ITS2 primers used in this study

PCR reactions were carried out in a total volume of 50μl by using standard protocols (Palumbi 1996). PCR temperature profiles to obtain the above-mentioned sequence were initial denaturation at 95°C for 3 min followed by 40 cycles of denaturation at 94°C for 40 s; annealing at 56°C for 40 s (ITS2 and COI) or 52°C (D2 and ND4); and extension at 72°C for 1 min and final extension at 72°C for 10 min.

For sequencing, PCR products were purified using polyethylene glycol (PEG) precipitation (20% PEG 8000 and 2.5 M NaCl). Sequencing reactions were carried out directly on both strands of DNA by using ABI Big Dye chemistry (Applied Biosystems, Foster City, CA), and the sequences were generated with an ABI 377 automated sequencer. The sequences were analyzed and questionable base calls resolved using Sequencher 3.0 (Gene Codes Corp., Ann Arbor, MI). Sequences were initially aligned using ClustalX, version 1.8 (Thompson et al. 1997) and then compared and visually aligned using Se-Al version 2.0a9 (Sequence Alignment Editor, A. Rambaut, University of Oxford) or MacClade (Maddison and Maddison 2000). GenBank accession numbers for all sequences are in Table 1.

For maximum likelihood (ML), we used PAUP 4.0b10 (Swofford 2004). Starting models were chosen with ModelTest (Posada and Crandall 1998) by using the Akaike Information Criterion (AIC). The resulting tree was saved and also used to test site-specific models. Some maximum likelihood determinations and data manipulations were done using p4 (Foster 2004). Maximum likelihood searches started with a neighbor-joining tree, on which we then optimized parameters and fixed the values for those parameters for branch swapping on that tree. Several ML search rounds were carried out until the parameters were fully optimized.

The program p4 (Foster 2004) was used to bootstrap the data, which allowed bootstrapping of partitioned data under a site-specific model. A consensus of the trees from 200 bootstrap replicates was made in PAUP. Branch lengths of that consensus tree were optimized using the search model and the original data, and bootstrap support values from the PAUP tree bipartitions table were placed on the tree using p4.

For Bayesian analysis, we used the program MrBayes (Huelsenbeck and Ronquist 2001). Where possible, the same model used for ML analysis was used in Bayesian analysis, however when that model was not implemented in MrBayes then the next more complex model available in MrBayes was used. When a sitespecific model was used, we used site rates, rate matrices, compositions, and among-site rate variation specific to each partition. Markov chain Monte Carlo (MCMC) runs were 500,000 generations long, sampling every 250 generations, for a total of 2,001 samples. Of these, the first 1,001 sample were discarded as burn-in, which is well past the point where the likelihood plot reached a plateau.

To estimate likely placements of the root for the four taxa, we used p4 to compile the sampled MCMC trees with outgroup attachment information preserved to show the distribution of root positions.

In total, 1,846 sites were included in the analysis (Table 3). Because the outgroup was very divergent, we suspected that it might cause us to choose inappropriate models or be responsible for long branch effects. To test this, the analyses were conducted with and without the outgroup. To better fit models to the data, we reasoned that the available genes could be separated into two groupings based on apparent relative substitution rates, “fast” and “slow”, with the mtDNA position 3 partition (fast) in one group and the remaining, mtDNA position 1 and rDNA (ITS2, D2), in the other group (slow). The data were not subdivided further because of lack of variation. The data groupings as described above as well as the number of parsimony informative characters are given in Table 3.

Table 3 Constant, variable, and parsimony informative sites in the ingroup only and with the outgroup included

Table 4 ML boostrap support and posterior probability (PP) for relationships within the An. albitarsis complex, when the outgroup was both included and excluded

Fig. 1 Results of combined data. The data were placed in two partitions as described in the text, consisting of mtDNA position 3 from COI and ND4 in partition 1, and mtDNA position 1 and ribosomal sequences D2 and (more ...)

Brief mention will first be made here of the separate analyses of fast (mtDNA position 3) and slow (mtDNA position 1 plus rDNA) partitions (Table 4). Using the mtDNA position 3 partition (fast) only, the topologies of the ML and Bayesian analyses were the same, and showed resolution of An. albitarsis and An. albitarsis B as separate groups, and good resolution of the (An. deaneorum, An. marajoara) group from the others (designated “A,”“B,”“D,” and “C,” respectively in Figs. 1 and 2). However, sequences of An. marajoara and An. deaneorum were not recovered into two separate nonexclusive clades (de Queiroz 1998) because An. marajoara (C10) and An. deaneorum (D17) clustered together in a poorly supported clade that was closer to An. marajoara than to An. deaneorum. Using the slow partition by itself (mtDNA position 1 and rDNA), there was more ambiguity, and the topology for the ML analysis differed from the Bayesian analysis. Support for separate groups was generally poor, with the highest support for the branch separating the entire An. marajoara sequence-group from the rest of the tree.

Fig. 2 Outgroup attachment distribution. From postburn-in samples from two MCMC runs using the combined data including the outgroup, a consensus tree was made such that the position of the attachment point of the (more ...)

Parsimony analysis of the combined rDNA and mtDNA data sets generated six most passinionious trees (MPTs) (not shown). The strict consensus tree generated from those six MPTs recovered three well-supported groups: An. albitarsis, An. albitarsis B, and (An. marajoara + An. deaneorum). The latter also were recovered as separate groups but with less support: An. deaneorum with 72% and An. marajoara with 91% bootstrap support.

ML and Bayesian analyses were carried out using the mtDNA position 3 in 1 partition, and mtDNA position 1 plus rDNA data in another partition. For the unpartitioned data, ModelTest suggested the TVM + I model; however, the TVM + SS, a site-specific model based on the two partitions, gave a better likelihood, showing an increase of 87.6 log units, so the TVM + SS model was used in all ML analyses. The ML topology, including the nonparametric bootstrap support values, is shown in Fig. 1. For Bayesian analysis, the data were partitioned in the same way. The general time-reversible (GTR) rate matrix and composition parameters were unlinked between the two partitions. The applicability of this model was confirmed by ML in p4, which showed an increase of 137 log units by allowing free rates and compositions in the two partitions, at a cost of eight parameters, compared with having the same rate matrix and composition in both partitions. As suggested by ModelTest, gamma-distributed among-site rate variation was applied to the mitochondrial position 3 partition, and a pInvar model was applied to the mitochondrial plus ribosomal partitions. The Bayesian tree has the same topology as the ML bootstrap tree. Bayesian posterior probabilities are shown in Fig. 1. Three groups (An. albitarsis, An. albitarsis B, and An. marajoara) were recovered with high support, and a fourth group, An. deaneorum, was recovered with somewhat less support. The support for the split between (An. albitarsis, An. albitarsus B) and (An. marajoara, An. deaneorum) was also high.

For the combined data, ModelTest suggested the GTR + I + G model. However, the GTR + SS model, where the data were partitioned into mitochondrial position 3 partition, and mitochondrial position 1 plus ribosomal partitions, had a better likelihood, with an increase of 86 log units, and so the GTR + SS model was used for ML analysis in PAUP. Using ML two major groups were recovered [(An. deaneorum, An. marajoara), (An. albitarsis, An. albitarsis B) outgroup]. Support for placement of the outgroup is low (51%), but the groups consisting of An. albitarsis and An. albitarsis B are moderately and strongly supported (92 and 99%, respectively) (Table 4). For the Bayesian analysis, a site-specific model was used, but a test was made in p4 using maximum likelihood to determine whether individual rate matrices and compositions in the two partitions were better than using a single overall rate matrix and composition. The increase in the log likelihood owing to a separate rate matrix and composition was 205, which is highly significant, and so this strategy was used in the Bayesian analysis. The settings were as described above for combined data without the outgroup. That includes using gamma model for the mitochondrial position 3 partition, and a pInvar model for the ribosomal partition. The support for relationships among the sequences of each ingroup taxon was generally lower when the outgroup was included in both ML and Bayesian analyses than when the outgroup was excluded (Table 4).

To estimate the distribution of likely root positions, the combined postburn-in MCMC samples from two runs were reanalyzed to obtain a consensus tree based on retained root information (Fig. 2). The largest number of input trees (667/2,000) had the outgroup attached on the branch separating the An. albitarsis B group from the rest of the tree. However, many other input trees had the outgroup attached on the branches leading to the (An. marajoara, An. deaneorum) clade, and many other input trees had the outgroup attached on the branch leading to the An. albitarsis clade. Also, many trees had the outgroup attached on the branches leading to sequences C10 (An. marajoara) and A4 (An. albitarsis), which we interpret to be caused by long branch effects and disregard. Few trees had the root attached along the branches leading to either the An. marajoara or An. deaneorum clades separately, showing little evidence for placement of the root on this part of the tree. However, note that the ML bootstrap tree has this root placement (not shown). Pending further observations we conclude that the rooting is as shown or on the branches leading to the (An. marajoara, An. deaneorum) clade or to the An. albitarsis clade. Stated another way, we find that the largest number of sampled trees had the outgroup attach between An. albitarsis B and [An. albitarsis (An. marajoara, An. deaneorum)], but with many trees having the outgroup attach at the base of the An. albitarsis clade, making [outgroup, (An. albitarsis, (An. albitarsis B, (An. marajoara, An. deaneorum)))] or on the branch separating An. albitarsis and An. albitarsis B from An. marajoara and An. deaneorum, making [outgroup ((An. albitarsis, An. albitarsis B), (An. marajoara, An. deaneorum))]. A topology with An. albitarsis B basal to the others was not recovered in any of the other parsimony, ML, or Bayesian consensus trees.

In a search for the best evolutionary hypothesis for the Albitarsis Complex, we used partial sequences of two mitochondrial genes (COI and ND4), and two ribosomal DNA fragments (ITS2 and D2 expansion region of the 28S subunit) and compared maximum parsimony, maximum likelihood, and Bayesian analyses with several combinations of data partitions. Individual genes failed to give well-resolved trees, possibly because of the low number of variable sites. Also, to optimize our results we analyzed different data partition combinations, finally settling on two partitions: mtDNA position 3 alone because it was more variable and presumably faster evolving, and mtDNA position 1 plus rDNA because they were less variable and presumably slower evolving. In addition, because the outgroup was highly divergent, we tested for long branch effects by carrying out all analyses with and without the outgroup.

The strongest support for the evolutionary relationships among the four species tested was retrieved when all four genes were combined and partitioned as described above. Analyses excluding the outgroup, presumably more independent of long branch effects, offer our best hypothesis for the ingroup topology (Fig. 1). In summary, four major evolutive lines were recovered. Two groups, An. albitarsis/An. albitarsis B and An. marajoara/An. deaneorum, were usually recovered, but not in all analyses. The latter group includes two species that are important vectors of human Plasmodium in localities in the Amazonas region of Brazil (Klein et al. 1991a, b; Conn et al. 2002). This suggests a possible phylogenetic link associated with the ability to transmit human malaria parasites. Similar conclusions have been made for An. albimanus and An. darlingi (Conn 1998).

When the mtDNA position 3 partition was analyzed alone, An. deaneorum specimen D17 and An. marajoara specimen C10 clustered together sister to the remaining An. marajoara individuals. This result suggests the possibility of incomplete lineage sorting, introgression (Donnelly et al. 2004), or even the existence of an additional taxon.

Our analysis of the distribution of possible roots used sampled trees from MCMC runs, including the outgroup. It showed that An. albitarsis and An. albitarsis B are in an ambiguous position with relation to the root but that (An. marajoara, An. deaneorum) generally formed a clade. Therefore, we think that the root is either as shown in Fig. 2 or as evidenced by the large number of MCMC sampled trees rooting there, either at the base of the An. albitarsis B clade, or at the base of the (An. marajoara, An. deaneorum) clade. A more definitive answer to this question requires more data.

Lehr et al. (2005) using the entire COI gene sequence data of 29 individuals of An. albitarsis, An. albitarsis B, An. deaneorum, and An. marajoara recovered results similar to those generated in the current study with the mtDNA position 3 data partition, in that there was nonexclusivity of the An. deaneorum clade with respect to An. marajoara. They showed the exclusivity of the sequences of An. albitarsis, An. albitarsis B, and a sister group relationship of these two taxa. They also found four individuals of An. marajoara that fell outside the remaining sequences of the Albitarsis Complex (in Bayesian topology), which they suggest represents a fifth species. The individuals that were used to generate these sequences were collected in Roraima State, Boa Vista, Brazil, and Venezuela. Interestingly, Lehr et al. (2005) also recovered a nonexclusive clade consisting of individuals of An. marajoara and An. deaneorum. A similar grouping was found when we analyzed the mtDNA data partition for the current study (An. marajoara specimen C10 clustered with An. deaneorum specimen D17). Lack of exclusivity of sequences of An. marajoara and An. deaneorum are similar to our results and are also suggestive of ancestral introgression or perhaps a recent speciation event that could not be detected by partial sequences of the mitochondrial genes COI and ND4.

The current study was based on conclusions about species boundaries reached using fixed RAPD markers (Wilkerson et al. 1995a, b). Our results corroborate the RAPD evidence that indicates four putative species: An. albitarsis, An. albitarsis B, An. marajoara, and An. deaneorum. No additional taxa were detected. This is not surprising given that the same genetic material was used in both studies. The existence of a fifth species as reported by Lehr et al. (2005) in Boa Vista, Brazil, was not directly tested by us using sequence as described in this study. However, we assumed the Boa Vista specimens to be An. marajoara based on comparison to one or two diagnostic RAPD markers (data not shown). The possibility of a fifth species is supported by independent data sets: Kreutzer et al. (1976) (chromosomes), Rosa-Freitas et al. (1990) (isozymes), and possibly Narang et al. (1993) (allozymes, mtDNA restriction fragment length polymorphisms). However, in support of a hypothesis for the existence of An. marajoara as a single widespread species is an extensive data set of rDNA ITS2 sequence from throughout its range, and taxon-specific PCR primers based on that sequence (Li and Wilkerson 2005). If in fact there is a fifth member of this complex, the RAPD results should be revisited because assumptions about the wide distribution of An. marajoara (Venezuela to southern São Paulo State) were based on the existence of seven RAPD markers that were found in nearly all individuals tested from all parts of its putative range (Wilkerson et al. 1995a, b).

The topology of the gene tree reported by Merritt et al. (2005), who used coding sequence of a portion of the white gene containing its fourth intron, varied significantly from that found by us and Lehr et al. (2005). They found good statistical support for a single loss of the fourth intron in the species complex (present in An. marajoara, absent in the other species), but weak evidence that An. marajoara is basal relative to [An. albitarsis B, (An. albitarsis, An. deaneorum)]. The alternative topology placed An. marajoara sister to An. An. albitarsis B. There was high support for the sister relationship of An. albitarsis and An. deaneorum. This is in contrast to our results that give high support for a close relationship between An. albitarsis and An. albitarsis B and between An. marajoara and An. deaneorum.

The above-mentioned conflicting results will certainly require additional data to resolve. A recent report by Besansky et al. (2003) addressed the issue of conflicting data sets in the resolution of species boundaries and phylogenetic relationships in the An. gambiae complex that may be germane in solving the issue of An. marajoara, and resolving the phylogenetic relationships of the Albitarsis Complex. They found evidence supporting both introgression and reproductive isolation as well as different tree topologies, depending on which sequence was sampled. They concluded that adoption of a “total evidence” approach for phylogenetic analysis of closely related species runs a risk of recovering a highly supported wrong answer and suggested that at the level of closely related species, it would be better to do a careful locus-by-locus assessment of sequence divergence rather than just adopt a total evidence approach. They were able to use various genes on all the chromosomes, inside and out of inversions, as well as mitochondrial genes, for their conclusions. The approach of Besansky et al. (2003) provides a model for future research into the phylogenetic relationships of the Albitarsis Complex.

We greatly appreciate the support of those who helped collect and rear many of the specimens used in this study, including Terry Klein and Eric Milstrey who were then members of the United States Army Medical Research Unit in Rio de Janeiro. But we especially thank Jose Bento Lima without whose extra efforts many of samples would not have been available to us. We also thank Lee Weigt at the Smithsonian Institution, Laboratory of Analytical Biology, for use of laboratory and sequencing facilities and Jan Conn and Meg Lehr for thoughtful review of a draft of this manuscript. This research was performed under a Memorandum of Understanding between the Walter Reed Army Institute of Research and the Smithsonian Institution, with institutional support provided by both organizations. This work was supported by National Institutes of Helath (AI) R0154139 to J. C.