Phylogenetic analysis and temporal diversification of mosquitoes (Diptera: Culicidae) based on nuclear genes and morphology

1Mosquitoes (Diptera: Culicidae) are a monophyletic group of true flies, recognizable by their elongate adult mouthparts through which the females of most species feed on vertebrate blood. Mosquitoes occur throughout temperate and tropical regions, and well beyond the Arctic Circle, but are most diverse in tropical forest environments. A bewildering amount of morphological diversity parallels their spectacular radiation into virtually every conceivable collection of water, ranging from a few droplets trapped by plant parts to large bodies of fresh and brackish surface water, making mosquitoes "as ubiquitous as water". The relationship between human health and those species that are medically important (<200 of 3,524 currently recognized; http://mosquito-taxonomic-inventory.info/) has driven most mosquito research. Within this small subset of disease vector species, morphological similarities between close relatives (e.g., cryptic or sibling species complexes) continue to pose practical and academic challenges to disease control, conventional taxonomy, and phylogenetic inference. Ironically, it is not the morphological similarity but rather the morphological diversity of mosquitoes that has confounded efforts to delimit many supraspecific groups and reconstruct their evolutionary history.

2The Culicidae had an ancient origin, probably in the Jurassic, consistent with the fossil record of their sister group Chaoboridae. Unfortunately, the sparse mosquito fossil record sheds no light on evolutionary relationships in the family. Traditional classification of Culicidae based on the phenetic framework of Edwards resulted in arbitrary groupings reflecting intuitive interpretation of morphological similarities. More modern classifications have incorporated important revisions of select genera and tribes based on explicit methodology, but in the absence of a comprehensive application of quantitative methods across the family, the result still does not entirely reflect evolutionary history [reviewed in]. Current classification divides the family into two subfamilies (Anophelinae and Culicinae), 11 tribes, and a minimum of 44 genera (Table 1).

3Generic-level relationships across all Culicidae have rarely been studied. The first attempt was based on comparative bionomics and morphology using intuitive methods typical of that time. Surprisingly, further attempts using modern cladistic methods were not made for nearly 50 years. The most comprehensive of these phylogenetic re-analyses employed 73 morphological characters to examine the relationships of the 38 genera then recognized. In general, almost none of the hypotheses raised by the 1951 phylogeny were supported in the 1998 reconstruction, with few exceptions including the monophyletic and basal position of the subfamily Anophelinae, and the monophyly of the tribes Sabethini and Culicini. However, most characters were homoplastic - some extensively - and many relationships were inadequately resolved. Although the Harbach and Kitching study challenged traditional generic groupings and reinforced the need for reappraisal, it also suggested that robust recovery of generic-level relationships of Culicidae may be difficult with morphological characters alone.

4Only four higher-level phylogenies of Culicidae based on gene sequences have been published, each of which were very taxon-limited in scope. All were able to show that, in agreement with the morphological phylogenies, Anopheles was sister to other sampled genera. The study by Miller et al., based on four mosquito species, was inconclusive and the more comprehensive study of Shepard et al., based on 18S rDNA sequences of 39 species representing nine genera, was unable to resolve deeper relationships. Of the four studies, only Besansky and Fahey employed a single-copy nuclear protein coding gene to assess relationships. Their study, based on the white gene, included 13 species representing nine genera. When third codon positions were excluded, Anophelinae was recovered as a basal lineage and Sabethini, Culicini and Aedini were recovered as monophyletic, suggesting the potential of protein-coding sequences for reconstructing generic-level relationships within Culicidae.

5The importance of sampling multiple genes when attempting to reconstruct species phylogenies is well recognized. Mitochondrial DNA and nuclear ribosomal DNA are convenient targets, due to conserved primer binding sequences and ease of amplification based on their typically high copy number. However, both can be problematic for resolving deep phylogenetic relationships. Mitochondrial DNA may exhibit a high mutation rate, skewed base composition, and even symbiont-induced biases. Beyond base compositional bias, ribosomal DNA also can be exceedingly difficult to align. Given the increasing availability of completely sequenced mosquito genomes, protein-coding nuclear genes represent a viable alternative as well as a rich and largely untapped resource.

6In the study reported here, we explored the phylogenetic utility of six nuclear protein-coding genes: arginine kinase, CAD, catalase, enolase, hunchback, and white. As noted above, only white was used previously in mosquitoes. All except catalase have been used in other insect groups: CAD in bees, empidoid flies, and lacewings, among others; arginine kinase in hymenopterans; enolase in beetles; and hunchback in Hawaiian drosophilids. We sequenced these six genes from 26 mosquito species representing 25 genera, and two chaoborid outgroup species. In addition, 80 morphological characters were scored from these mosquito and outgroup species. Our goals were twofold: (1) to estimate a generic-level phylogeny of Culicidae based on molecular and morphological evidence, and (2) to use fossils and sequence data to infer divergence times for major culicid lineages.

7Twenty-six species of mosquitoes representing 25 genera were used as ingroup taxa (Table 1). Two chaoborid midges, Eucorethra underwoodi and Chaoborus astictopus, were used as outgroup taxa based on the sister-group relationship between Chaoboridae and Culicidae. Specimens were preserved in 70-100% ethanol at -20°C. Sampling of additional ingroup genera for this study was precluded by unavailability of specimens adequately preserved for molecular analysis; DNA from pinned museum specimens or specimens stored at room temperature for prolonged periods was found to be excessively degraded.

8Sequences were obtained from six nuclear protein-coding genes (arginine kinase, CAD, catalase, enolase, hunchback, and white) from VectorBase http://www.vectorbase.org for the three mosquito species with completely sequenced genomes (Anopheles gambiae, Culex quinquefasciatus, and Aedes aegypti), from GenBank for available white sequences (GenBank accession numbers U73829, AF318199, AF318200, U73834, U73827, AF318206, U73835, U73837, AY055811, AF318193, AF318209), or by PCR amplification and direct sequencing. Genes were amplified and sequenced using degenerate primers previously reported in the literature, or designed based on amino acid alignments of respective genes from the three sequenced mosquito species. Primer sequences and their sources are provided in Table 2. Primers were located in regions that would allow amplification within a single large exon, or across exons separated by small introns, to facilitate amplification from genomic DNA templates. Primers were hemi-nested, whereby the first round of amplification based on the outermost pair was followed by alternative second round PCR from which two internal, overlapping fragments were amplified (Figure 1).

9Genomic DNA was extracted from whole mosquitoes using the DNeasy kit (Qiagen, Valencia, CA). PCR reactions (50 μl) contained 20 mM Tris-HCl (pH 8.4), 50 mM KCl, 1.5 mM MgCl2, 200 μM each dNTP, 2.5 U Taq polymerase, 2.4% DMSO, 0.5% BSA, 25 pmol each primer, and 1 μl template DNA (~1/50th-1/500th of the amount extracted from a single mosquito). Amplifications consisted of 1 cycle at 93°C for 1 min; 35 cycles of 94°C for 20 sec, 45°C for 20 sec, and 72°C for 2 min; and a final extension cycle of 72°C for 5 min. After inspection of an aliquot by electrophoresis through a 1.5% agarose gel, products of the expected size were excised from the gel and purified using GeneClean Kit (MP Biomedicals, Irvine, CA), or sequenced directly if only one band was observed. Excess primers and dNTPs were removed by adding 2 U of exonuclease 1, 1 U of shrimp alkaline phosphatase, and 1.8 μl of H2O to 8 μl of PCR product, incubating at 37°C for 15 min, and inactivating at 80°C for 15 min before sequencing. Sequencing was carried out using an Applied Biosystems 3730 × l DNA Analyzer and Big Dye Terminator v3.1 chemistry. Electropherograms were inspected and trimmed using SeqMan II (DNASTAR, Madison, WI). Sequences were deposited in GenBank under accession numbers: arginine kinase, GQ906806-GQ906829; CAD, GQ906830-GQ906854; catalase, GQ906855-GQ906874; enolase, GQ906879-GQ906901; hunchback GQ906902-GQ906926; white, GQ906927-GQ906937.

10Morphological structures were examined in the adult, pupal and larval (fourth-instar) stages. Heads were removed from adult mosquitoes for comparative studies of structures not readily visible in intact specimens. Heads were cleared in 5% sodium hydroxide solution, stained in acid fuchsin and mounted frontodorsal side uppermost in Euparal on microscope slides. Pinned adults were examined under simulated natural light. Dissected genitalia, larvae, and larval and pupal exuviae were studied with differential interference contrast optics. The morphological terminology follows Harbach and Knight and Harbach and Kitching.

11The 28 species were coded for 80 characters (see additional file 1: FileS1) derived from fourth-instar larvae (24), pupae (12), adults (36), female genitalia (1), and male genitalia (7). The data matrix is shown in Additional file 2: TableS2. Characters were coded from direct observations except in a few cases where structures were missing from available specimens. Some missing data for the immature stages of Chaoborus astictopus were coded from literature sources. Characters that could not be scored due to missing data or absence of homologous structures (e.g., character 20 in three anophelines and two chaoborid species that lack a siphon) were denoted by a "?". All multistate characters were treated as unordered. Most of the structures and states of characters are illustrated in Harbach & Knight and Harbach & Kitching.

12Nucleotide sequence alignments were guided by the corresponding amino acid alignments, using utilities within the program suite EMBOSS; http://embossgui.sourceforge.net/demo/. Inferred amino acid sequences were aligned using the program "emma," which provides an interface to ClustalW. The Gonnet matrix was chosen and the resulting alignment was followed by limited manual adjustment. "Tranalign" was then used to align the nucleotide sequence based on the previously aligned amino acid sequence. Introns were removed from aligned gene regions before analysis.

13Basic sequence information (pairwise sequence divergence, base composition, statistical tests of homogeneity of base composition, number of variable and parsimony informative characters) was obtained using PAUP* v4.0b10. In addition, plots of transitions and transversions versus divergence at each codon position were based on observed (uncorrected) p-distances from PAUP*.

14Maximum parsimony (MP) analyses were implemented in PAUP* on a phylogenetic data set that included concatenated genes with/without morphological characters. Third codon positions for each gene were removed and gaps were treated as missing data. Heuristic searches consisted of 1000 random sequence additions with tree bisection-reconnection (TBR) branch swapping. Bootstrap support values were based on 500 replicates, each with 10 random additions and TBR branch swapping.

15Maximum likelihood (ML) analyses were performed on molecular data sets only, which included both individual/concatenated genes, with/without third codon positions. The ML heuristic searches were performed in PAUP*, using the model of nucleotide substitution and parameter values selected via Modeltest. Values for the substitution matrix, base composition, gamma distribution of among-site rate variation (G) and the proportion of invariant sites (I) are available from the authors on request. Bootstrap resampling was conducted using 1000 replicate neighbour-joining (NJ) trees based on the ML substitution matrix. The Shimodaira-Hasegawa test [; data not shown] was used to test for incongruence between phylogenies suggested by individual genes, or successive combinations of congruent genes.

16Bayesian (BI) phylogenetic tree searches were performed in MrBayes 3.1.2 on concatenated gene and gene + morphology data sets using aligned nucleotides, both including and excluding third codon positions, and using concatenated aligned amino acids. For concatenated genes (nucleotides) and genes (nucleotides) + morphology, a mixed model approach was used with model parameters specified per gene partition according to Modeltest and a Markov K + G model for morphology, with branch lengths unlinked and estimated for each partition. Bayesian tree searches using aligned amino acids were carried out in MrBayes 3.1.2 using the WAG model of amino acid evolution (WAG+I+G). Each Bayesian search was carried out for 10,000,000 generations (sampling every 1000) using four chains (default heating parameters) and a 30% burn-in value. The included Bayesian sets of trees were sampled after likelihood scores reached convergence and the mean split difference values were below 0.02.

17Estimates of divergence times for mosquito lineages were calculated using the parametric Bayesian-relaxed clock approach implemented in the programs ESTBRANCHES and MULTIDIVTIME and using the combined gene data set (nucleotides) including third codon positions. Branch lengths and evolutionary rate priors were estimated from the data using the BASEML program in the PAML software package and ESTBRANCHES. Tree topology, minimum and maximum root node age, and fossil-based minimum age constraints are set as user-defined analysis priors. For the tree topology we used the tree recovered from the BI search of combined amino acids (see Phylogenetic Analyses, below). The root max-min age prior between Culicidae and outgroups was set as 230-187 Ma corresponding to the hypothesized age of the Diptera and a fossil assignable to the Chaoboridae [187 Ma; ref], and three lineages were constrained according to fossil-based minimum ages (Toxorhynchites mexicanus, 16 Ma; Culex winchesteri, 34 Ma; Anopheles dominicanus, 34 Ma; http://mosquito-taxonomic-inventory.info/category/fossil-culicidae/fossil-culicidae). We followed the analytical procedure described in Rutschmann et al. and in the MULTIDIVTIME readme files, and ran the Markov chain for 1.1 × 106 cycles with samples collected every 100 cycles and discarded the first 100,000 cycles as burn-in. We performed the MULTIDIVTIME analysis multiple times from different initial conditions to confirm convergence of the Markov chain on highly similar resulting time estimates and posterior intervals.

18Across six genes, the molecular data matrix consisted of 5352 aligned characters, of which 2839 were variable and 2259 were parsimony informative (Table 3). Not surprisingly, most of the variation was found in the third codon (nt3) position. Analysis of base composition of the combined genes for each major taxonomic grouping revealed significant departures from homogeneity at the nt3 position, owing to three groups: Anophelinae, Culicini (Cx. quinquefasciatus) and Ficalbiini (Mimomyia luzonensis) (Table 4). Moreover, plots of transitions and transversions against uncorrected pairwise nucleotide divergences at each codon position suggested saturation of transitions at the nt3 position (Figure 2). These results prompted us to perform ML and BI phylogenetic analyses both with and without the nt3 partition, and to exclude this partition in MP analyses (see Phylogenetic analyses).

19Uncorrected pairwise sequence divergence within the mosquitoes sampled ranged widely, from 10-27% (summarized by major taxonomic groupings in Table 5). However, the average distance between tribes and subfamilies was at the upper end of the range (20% and 23%, respectively), approaching that between mosquitoes and their sister group, the chaoborid midges (26%).

20Evidence was found for incongruence among some genes or gene combinations via the Shimodaira-Hasegawa test; data not shown. However, examination of phylogenies resulting from individual genes revealed that topological incongruence was generally limited to certain poorly supported nodes. It has been suggested that different data sets may have a common phylogenetic signal recoverable only upon combined analysis, and, under the hypothesis that combining data from multiple genes may potentially overcome misleading signal in individual genes, we conducted further analyses to compare results from a concatenated data set with results obtained from ML and BI.

21Relationships inferred by MP and ML, including BI, are summarized in Figure 3 and Table 6. All three algorithms, as applied to various data partitions (± nt3; ± morphological characters; nucleotides or amino acids), gave overwhelming support for the monophyly of Culicidae (node O), the monophyly and basal position of the subfamily Anophelinae (node A; gray box), and the monophyly of the tribe Sabethini (node I; gray box). Less conclusive support by ML, but reasonable support by MP and BI, was observed for the monophyly of the tribe Aedini (node C). These results confirm the conclusions of Harbach regarding what was already known about the phylogeny of mosquitoes.

22Support varied for relationships within these well-supported clades. The subfamily Anophelinae was represented in this study by three species from two genera: Bironella (Bironella gracilis) and Anopheles (Anopheles atroparvus, subgenus Anopheles; An. gambiae, subgenus Cellia). Not all analyses supported the monophyly of the genus Anopheles; an alternative relationship of Bi. gracilis + An. atroparvus was also recovered. There is precedence for the paraphyly of Anopheles relative to Bironella in previous morphological and molecular studies. Reliable inference of relationships between these groups may be problematic due to conflicting signals or contemporaneous radiations, but the suggestion of Sallum et al. to redefine Bironella as an informal group within Anopheles seems premature.

23Within tribe Sabethini, Malaya occupied the most basal position among the taxa sampled, although this placement was not recovered in a subset of BI analyses. The genus Maorigoeldia, containing only a single species exclusive to New Zealand, was sister to Tripteroides (Oriental, Australasian and Palaearctic species), in all cases with 100% bootstrap support or posterior probability of 1.0. These relationships were not recovered in cladistic analyses of morphological data that included representatives of a larger number of genera. Belkin regarded Maorigoeldia to be sister to all other sabethine species. This was supported in the studies of Harbach and Kitching, Harbach and Peyton, and Harbach et al., but not in the study of Judd, which placed Maorigoeldia as the sister group to the New World genera of Sabethini. Whereas Malaya was recovered as the sister of genus Topomyia in the first two of these four studies, it was paired with Limatus in the most derived clade of Sabethini when Harbach et al. included the new genus Kimia in the data set of Harbach and Peyton. Also, Tripteroides (Old World) was recovered as sister to Trichoprosopon (New World), which is supported by shared morphological characters that are unique to these two genera.

24Decisive support for monophyly of the New World genera (Limatus, Sabethes, Shannoniana, Trichoprosopon, and Wyeomyia) was found, in agreement with previous studies. Among these genera, a close relationship between Limatus, Sabethes and Wyeomyia was strongly supported by all analyses, but other nodes were unstable. Although not apparent in Figure 3, Shannoniana and Trichoprosopon showed a sister relationship in all but the BI amino acid analysis. However, as indicated in the previous paragraph, the results of cladistic analyses based on morphological data and a larger sample of sabethine taxa casts doubt on these relationships.

25The remarkably large tribe Aedini (1255 species, http://mosquito-taxonomic-inventory.info/taxonomy/term/6065) has been the subject of recent efforts to infer higher-level relationships based on morphological characters of all life stages. Although this has resulted in major changes to classification, phylogenetic resolution has been limited. In the present study, as in the cladistic analyses of extensive morphological data by Reinert et al., Psorophora was recovered as sister to all other Aedini. Sister-group relationships strongly supported in most cases were Aedes (Stegomyia) + Eretmapodites and Haemagogus + Ochlerotatus. Other relationships within Aedini were less clear. Moreover, no consensus could be reached regarding affinities of any other genera within Culicidae as a whole, outside of Aedini and Sabethini.

26A chronogram for Culicidae is given in Figure 4, and corresponding divergence time estimates are provided in Table 7. Based on the taxa sampled and three fossil constraints, earliest divergence within mosquitoes - between the lineages leading to Culicinae and Anophelinae - dates to ~226 Ma. This estimate is in reasonable agreement with Krzywinski et al., who determined the split between Anopheles and Aedes (Stegomyia) to have occurred ~145-200 Ma based on mitochondrial DNA sequences. Although 226 Ma is substantially older than the 118 Ma divergence between Chaoboridae and Culicidae estimated by Bertone et al., the 95% credibility ranges overlap between studies. Because the latter study was aimed at deeper divergences within lower Diptera, Bertone et al. only included two mosquitoes and one chaoborid, possibly accounting for the discrepancy. Moreover, Bertone et al. estimated divergence times from a single gene (28S rDNA) and used only a few fossil calibration points (none close to Culicidae), which may also have contributed to differences in age estimates. We favor the older divergence estimates, as they are consistent with other evidence suggesting that mosquitoes likely originated in the Jurassic. As early as 1923, Edwards surmized that the "origin and phylogenetic history of the Culicidae must go back to well into the Mesozoic Era."

27Unexpectedly, the split between Anopheles and Bironella was a remarkably shallow 51 Ma. This estimate contradicts Krzywinski et al., who estimated a substantially greater (90-106 Ma) divergence between the Anopheles subgenera, Anopheles and Cellia, represented in this study by An. atroparvus and An. gambiae, respectively. To what extent this incongruence can be explained by limited taxon sampling or biases in the molecules used to infer divergence dates (e.g., strong base composition bias in Anophelinae; Table 4), is presently unknown. However, given the universal agreement that Anopheles occupies an early-branching position among the Culicidae, it seems likely that our crown-group estimates do not accurately reflect the age of this group.

28The only well-supported clades of the subfamily Culicinae in the phylogenetic analyses are the tribes Aedini and Sabethini, which apparently arose at similar times (roughly 112 and 115 Ma, respectively) and diversified more recently. Eight genera of the subfamily Culicinae (Aedeomyia, Coquillettidia, Culiseta, Culex, Mimomyia, Orthopodomyia, Toxorhynchites, Uranotaenia), whose relationships were not strongly or consistently recovered, represent the deeper internal branches of the tree. The nodes connecting these branches are not only ancient (exceeding 127 Ma), but also relatively close together in time, occurring within a ~30 million year interval between 127-158 Ma. If these estimates are corroborated in the future by denser taxon sampling and more fossil-based age constraints, they will support the notion of rapid radiations early in the diversification of mosquitoes, potentially explaining the difficulty in attaining a stable phylogeny for these lineages. An early Cretaceous timing for these rapid radiations is consistent with the appearance of angiosperms, a group of plants whose nectar is exploited as an energy source by mosquitoes, and whose water-filled parts are the sole habitats occupied by the immature stages of many groups of mosquitoes, notably members of the tribe Sabethini.

29This study represents one of the few attempts to reconstruct generic-level relationships within Culicidae as a whole, and the only attempt to combine morphological data and molecular characters from multiple genes. Among molecular phylogenetic studies of the family, it more than doubled the number of taxa sampled to date. Yet results were mixed. The ability to recover previously known clades (Anophelinae, Sabethini, and Aedini) was encouraging. However, the deeper relationships among genera could not be resolved unambiguously, potentially due to ancient and rapid radiation, as hypothesized for other insect groups. There has been much debate regarding whether better resolution and support of relationships is achieved through broader taxonomic sampling or sequencing of more loci. The current explosion of sequencing whole genomes from organisms, including mosquitoes, promises many potentially informative genes beyond those included here. On the other hand, the recent study by Wiegmann et al. successfully resolved even deeper divergences and a longstanding controversy in the phylogeny of holometabolous insects, using only six single-copy nuclear genes comprising a similar number of base pairs to that compiled for the present study. Although more molecular, as well as morphological, characters may well prove useful, there is little doubt that broader taxonomic sampling is now the key roadblock. Considering that the mosquito diversity housed in museums is almost invariably preserved in a fashion that has impeded conventional molecular data collection, this roadblock may be substantial. There is an urgent need for fresh museum collections, particularly from under-sampled yet high-biodiversity regions worldwide, and their cryo- or ethanol preservation with vouchers. The present study was limited by what was available in existing collections. Broader taxon sampling is crucial not merely because it may help break up long branches. To the extent that the current generic system of classification includes paraphyletic and polyphyletic groups containing numerous species, it is clear that inclusion of only one or few generic exemplars can be misleading, and that more representative sampling is needed.

30BI: Bayesian; Ma: million years ago; ML: maximum likelihood; MP: maximum parsimony; nt: nucleotide; PCR: polymerase chain reaction.

31NJB and REH conceived and designed the experiments. KRR designed PCR primers and collected molecular data. REH collected morphological data. KRR, SC, MAB and BMW performed experiments and analyzed data. REH contributed specimens. KRR and NJB wrote the paper with contributions by all authors.