Evidence for a Population Expansion in the West Nile Virus Vector Culex tarsalis

Abstract

Population genetic structure of the West Nile Virus vector Culex tarsalis was investigated in 5 states in the western United States using 5 microsatellite loci and a fragment of the mitochondrial reduced form of nicotinamide adenine dinucleotide dehydrogenase 4 (ND4) gene. ND4 sequence analysis revealed a lack of isolation by distance, panmixia across all populations, an excess of rare haplotypes, and a star-like phylogeny. Microsatellites revealed moderate genetic differentiation and isolation by distance, with the largest genetic distance occurring between populations in southern California and New Mexico (F_ST = 0.146). Clustering analysis and analysis of molecular variance on microsatellite data indicated the presence of 3 broad population clusters. Mismatch distributions and site-frequency spectra derived from mitochondrial ND4 sequences displayed pattern's characteristic of population expansion. Fu and Li's D* and F*, Fu's F_S, and Tajima's D statistics performed on ND4 sequences all revealed significant, negative deviations from mutation-drift equilibrium. Microsatellite-based multilocus heterozygosity tests showed evidence of range expansion in the majority of populations. Our results suggest that C. tarsalis underwent a range expansion across the western United States within the last 375,000–560,000 years, which may have been associated with Pleistocene glaciation events that occurred in the midwestern and western United States between 350,000 and 1 MYA.

Introduction

Culex tarsalis is implicated as a major mosquito vector of West Nile Virus (WNV) in the western United States (Goddard et al. 2002, Reisen et al. 2004, Bell et al. 2005, Turrell et al. 2005, DiMenna et al. 2006). The geographic range of C. tarsalis spans from the west coast to the Mississippi River and from northern Mexico to Canada (Darsie and Ward 1981). Larvae can develop in various freshwater habitats (Reisen and Reeves 1987) and adult populations commonly reach very high density, becoming important biting nuisance pests. Historically, C. tarsalis has played a prominent role in the maintenance and transmission of arboviruses including St. Louis Encephalitis Virus and Western Equine Encephalitis Virus. Laboratory studies have revealed that C. tarsalis is among the most efficient vectors of WNV known, with high rates of both oral and vertical transmission (Goddard et al. 2002, 2003). Culex tarsalis not only feeds primarily on birds (Reisen and Reeves 1987) and is hypothesized to play an important role in WNV amplification, but also feeds opportunistically on mammals (Reisen and Reeves 1987) and is abundant in peridomestic habitats, possibly serving as a bridge vector to humans and domestic animals. Given its broad range, large population sizes, and high capacity for WNV transmission, C. tarsalis is likely in part responsible for recent epidemics in western United States including outbreaks in Nebraska, Colorado, Arizona, New Mexico, and California (Hayes et al. 2005).

Studies of field-caught individuals indicate that there is significant spatial and temporal heterogeneity in the ability of C. tarsalis to transmit WNV. Northern California populations exhibited higher transovarial and oral transmission than populations in southern California (Goddard et al. 2002, 2003). WNV infection rates also appear to vary at different times of year (Reisen et al. 2006). The extent to which these observed temporal and spatial differences in WNV transmission are rooted in genetic variation and/or environmental influences is currently unknown. It is possible that part of the variation that C. tarsalis exhibits in viral transmissibility is due to genetic differences among populations. However, except for a single allozyme study (Gimnig et al. 1999), little is known about the genetic structure of C. tarsalis populations. This study was undertaken to begin characterizing population structure of C. tarsalis in the western United States.

We have characterized the genetic structure of C. tarsalis populations across a broad geographical range using both mitochondrial and nuclear markers. We conducted a genetic analysis of 12 C. tarsalis populations from California, Washington, Colorado, New Mexico, and Nebraska using 5 microsatellite markers and mitochondrial sequence data. Microsatellites revealed moderate genetic structure and isolation by distance, whereas mitochondrial sequences indicated panmixia across the study region. Both types of markers exhibit signatures indicative of a range expansion. Our findings have important implications for future studies on C. tarsalis population dynamics, the potential role of genetic structure in explaining variation in viral transmission, and the effect of past evolutionary events on current estimates of gene flow.

Materials and Methods

Mosquito Collections and Extractions of DNA

Adult mosquitoes were collected in California, Washington, New Mexico, Colorado, and Nebraska between June and August 2005 using Centers for Disease Control light traps (fig. 1). Specimens were either transported from the collection localities to the Johns Hopkins Bloomberg School of Public Health on dry ice or preserved in 100% ethanol and stored at −80° C until processed for DNA extraction. Specimens were identified as C. tarsalis according to Darsie and Ward (1981). DNA was extracted from individual mosquitoes using the DNeasy Tissue Kit (Qiagen, Valencia, CA) according to Appendix G of the manufacturer's protocol for insects or by salt extraction/ethanol precipitation as previously described (Black and DuTeau 1997). Extracted DNA was suspended in nuclease-free water and DNA concentration was quantified using a NanoDrop spectrophotometer (NanoDrop Technologies, Wilmington, DE), adjusted to 25–50 ng/μL, and stored at −20° C until used for Polymerase Chain Reaction (PCR).

ND4 Mitochondrial DNA Amplification and Sequencing

A 389 bp fragment of the reduced form of nicotinamide adenine dinucleotide dehydrogenase gene subunit 4 was amplified using primers ND4+ and ND4− as previously described (Gorrochotegui-Escalante et al. 2000). Amplifications were carried out using DNA from 25 to 32 individuals from each California population and 16 individuals from 1 population each in Washington, Colorado, and Nebraska, for a total of 8 populations (N = 170). PCR conditions were as previously described (Gorrochotegui-Escalante et al. 2000). Amplicons were separated by 1% agarose gel electrophoresis, purified from the gel using Qiaquick columns (Qiagen), and sequenced in both directions at the Johns Hopkins Medical Institutes Core Sequencing Facility (Baltimore, MD). Because primer sequences contain degenerate sites, these portions were excluded from sequence analysis.

Microsatellite Amplification and Fragment Size Determination

Five trinucleotide repeat microsatellite loci were selected from 12 previously identified markers for analysis (CUTC6, CUTC12, CUTD107, CUTD113, and CUTD120) (Rasgon, Venkatesan et al. 2006). Primer sequences and PCR conditions were as previously described (Rasgon et al. 2006). Amplifications were carried out using DNA from 24 to 38 individuals from each of 12 populations (N = 356). Forward primers were 5′ labeled with either HEX or 6FAM for marker multiplexing. Labeled amplicons were resolved on an ABI Prism Genetic Analyzer 3100 Avant (Applied Biosystems, Foster City, CA). Allele sizes were automatically estimated with an internal ROX-500 size standard (Applied Biosystems) using GeneScan v. 3.1 and Genotyper software (Applied Biosystems).

Diversity and Population Structure Analyses

ND4 haplotype polymorphism and nucleotide diversity, estimated by π (the average number of nucleotide differences per site between 2 sequences, calculated according to Nei (1987), were determined using DnaSP 4.0 (Rozas et al. 2003). Microsatellite allele frequencies were analyzed for deviation from Hardy–Weinberg equilibrium and linkage disequilibrium using Arlequin 2.0 (Schneider et al. 2000) and GenePop software (Raymond and Rousset 1995). Genetic diversity values were not normally distributed, so diversity among populations was compared using the Wilcoxon matched-pairs sign-rank test for microsatellite loci and a Chi-square goodness of fit test for ND4 sequences.

Arlequin 2.0 was used to compute F_ST values (Wright 1951, Weir and Cockerham 1984) and significance determined by permutation tests (1,000 permutations). An unweighted pair group method with arithmatic mean (UPGMA) phylogenetic tree inferred from microsatellite-derived multilocus F_ST values was constructed using MEGA v. 3.1 (Kumar et al. 2004). Branch lengths were determined using the program PHYLIP (Felsenstein 1989). A statistical parsimony (Templeton et al. 1992) haplotype network of unique ND4 sequences was constructed using the program TCS 1.21 (Clement et al. 2000) with a 95% cutoff.

Structure 2.0 (Pritchard et al. 2000) was used to assign individuals from all populations to a predetermined number of clusters (K) based on multilocus microsatellite data. For each run, a burn-in period of 30,000 steps was followed by 1,000,000 iterations under the admixture model and the assumption of correlated allele frequencies among populations (gamma = 0.3468, calculated in 5 runs at K = 1). For each K of 1–6, 10 runs were performed. Log-likelihood scores were averaged across runs and compared to determine the posterior probability of each K. Arlequin 2.0 was used to partition genetic variation among regions of restricted gene flow (based on microsatellite-derived clustering results and F_ST values) using analysis of molecular variance (AMOVA). Regressions of linearized F_ST values on ln-transformed geographic distances among populations (Rousset 1997) were performed for microsatellite and ND4-derived data using STATA 8.0. Mantel probabilities were calculated using IBD 1.52 (Bohonak 2002).

Methods to Infer Population Expansion

Several test statistics have been developed to determine whether microsatellite and mitochondrial haplotype frequencies deviate from evolutionary neutrality. These statistics are based on a comparison of observed variability (measured by heterozygosity, haplotype or nucleotide diversity, and number of pairwise differences or segregating sites) to what would be expected under mutation-drift equilibrium or population growth. H and K tests (Depaulis and Veuille 1999) were performed on ND4 sequence data using the program Allelix (provided by S. Mousset). These tests compare observed haplotype diversity and number with an expected distribution of values under equilibrium conditions. Strobeck's S statistic (Strobeck 1987), the probability of having an equal number or fewer haplotypes than observed based on the gene frequency distribution derived from the inferred mutation rate θ, was computed for each population using the program DnaSP v. 4.0. Mismatch distributions (Rogers and Harpending 1992, Rogers 1995) and site-frequency spectra (Donnelly et al. 2001) produce characteristic patterns believed to reflect population expansion. Distributions and frequency spectra were plotted across all samples and for each population using DnaSP v. 4.0 along with the corresponding raggedness statistic (Harpending 1994) to quantify the smoothness of mismatch distributions.

Fu and Li's D* and F * (Fu and Li 1993), and Fu's F_s statistics test for deviation of sequence variation from evolutionary neutrality. Fu's F_s is based on the probability of the observed number of haplotypes or greater occurring under conditions of neutrality, whereas Fu and Li's D* and F* compare estimates of theta based on mutations in internal (older) and external (younger) branches of a genealogy. Negative values of these statistics are most often attributed to positive selective sweeps (such as genetic hitchhiking), recent population growth, or background selection. These events may all result in similar sequence patterns such as an excess of rare haplotypes, but different tests exhibit varying levels of sensitivity in detecting potential causes. Used in combination, these tests can provide evidence for or against particular evolutionary mechanisms. Fu and Li's D* and F* and Fu's F_s tests and confidence intervals (CI) calculated from coalescent simulations were conducted using DnaSP v. 4.0.

A multilocus test of expected homozygosity (and conversely heterozygosity) has been developed for microsatellite evolution under a variety of mutation models (Cornuet and Luikart 1996). Expected heterozygosity under Hardy–Weinberg allele frequencies (H_e) is compared with expected heterozygosity calculated from the expected number of alleles and sample size (H_eq). Under conditions of rapid population growth, high mutation rates and low loss of variation from drift would cause both heterozygosity and allele number to increase, though allele number would increase at a more rapid rate. As a result, H_eq would be greater than H_e (Cornuet and Luikart 1996), providing a signature of population expansion. Microsatellites generally follow the stepwise mutation model (SMM) or have a small percentage of multistep mutations (2-phase model or TPM) (Weber and Wong 1993, Di Rienzo et al. 1994, Primmer et al. 1998). Tests under strict SMM and TPM with 5%, 10%, and 15% multistep mutations were performed using the program Bottleneck (Cornuet and Luikart 1996).

Results

ND4 Haplotype and Nucleotide Polymorphism

Sixty-four ND4 haplotypes were detected in 170 C. tarsalis individuals across 8 populations. Forty-nine of 349 nucleotide positions were variable and only 5 substitutions were nonsynonymous. The number of haplotypes per population ranged from 9 to 17, with a mean of 12.75 (table 1). Once haplotype richness was weighted by sample size, no significant differences in haplotype number were observed. Haplotype diversity (Hd) ranged from 0.78 to 0.96 among populations and nucleotide diversity (π) ranged from 0.0050 to 0.0079. Haplotype diversity (Hd) and nucleotide diversity (π) did not vary significantly among populations.

H and K tests indicated that although overall ND4 haplotype diversity did not deviate from the expected distribution (observed Hd = 0.88; 95% CI of expected Hd = 0.81–0.94; P = 0.37), the number of observed haplotypes, at 64, is much greater than what would be expected under the neutral model (95% CI of expected number of haplotypes = 16–35; P < 0.0001).

Hardy–Weinberg Equilibrium, Linkage Disequilibrium, and Microsatellite Polymorphism

Significant deviation from Hardy–Weinberg equilibrium was found in 8 of 60 possible tests (Bonferroni adjusted α = 0.0008; table 2). All deviations were due to heterozygosity deficiency. None of the pairs of loci exhibited significant linkage disequilibrium (Bonferroni corrected α = 0.0004) in single pairwise tests. Fisher's global test for each pair of loci across all 12 populations also did not reveal significant linkage disequilibrium.

The 5 microsatellite loci had an average of 5.3 alleles per locus per population and ranged between 3 and 10 alleles in each population. Average heterozygosity per locus ranged between 0.12 and 0.71. No significant differences in allele number or heterozygosity were detected among populations.

Phylogenetic Analyses

An ND4 haplotype network was constructed using statistical parsimony with a 95% connection limit (parsimony cutoff = 7 mutational steps) (fig. 2). Approximately 1/3 of the individuals in each population shared haplotype X, which is inferred to be the ancestral haplotype (fig. 2). A star-like pattern of closely related haplotypes appeared to radiate from haplotype X. Of the remaining 63 haplotypes, 47 were unique, occurring in only one individual from the entire sample. Geographic clustering of haplotypes was not observed.

Population pairwise F_ST values based on microsatellite allele frequencies were used to generate a UPGMA population phylogeny, in which a constant rate of evolution across populations is assumed and genetic distances between populations are relative. Based on an F_ST cutoff of 0.05, considered to represent moderate genetic isolation (Wright 1978), 3 population clusters are present. These include 1) 4 California populations (CTM, CTG, CTY, and CTK) along with Colorado and Washington (CTOa, CTOb, CTWa, and CTWb), 2) New Mexico (CTNM) and Nebraska (CTNa and CTNb), and 3) population CTC in California (fig. 3). Population pairwise F_ST values based on ND4 haplotype frequencies ranged from 0.00 to 0.008. Because pairwise ND4 F_ST values were very low and nonsignificant, a phylogenetic tree constructed from these data is not presented here.

Population Genetic Structure

The posterior probability for K = 3 clusters in the microsatellite-based multilocus clustering analysis was 0.95. Posterior probabilities for K = 1, 2, and 4–6 were nearly 0. At K = 3, population CTC (Coachella, CA) appeared to be distinct, as did the cluster of CTNM, CTNa, and CTNb (fig. 3), corroborating results from the F_ST-based UPGMA tree. The remaining populations appear to share similar proportions of ancestry from each of the other 2 groups as well as from a 3rd group. In all cases, high levels of admixture were present within groups and individuals.

Microsatellite-based phylogenetic and cluster analyses were used to designate 3 genetically distinct groups of populations (fig. 3). An AMOVA conducted by partitioning variation among and within groups and populations revealed that ND4 haplotype variation was attributed solely to within-population differences (100%) and not to differences among populations or the 3 geographic regions delineated by microsatellite-derived genetic distances (F_CT = −0.002, P = 0.45). In contrast, the microsatellite-based AMOVA indicated that 93% of observed variation in microsatellites was attributed to within-population differences. Approximately 6% of the variation in microsatellite allele frequencies was associated with differences among the 3 regions (F_CT = 0.055, P < 0.001), whereas the remaining variation (<2%) was attributed to among-population differences (F_ST = 0.073, P < 0.001).

Gene Flow among Populations

The Mantel regression for ND4 sequence data (fig. 4) did not yield a significant relationship between genetic and geographic distance (Mantel probability P = 0.823). The Mantel regression of linearized microsatellite F_ST values on the natural logarithm of pairwise distances among collection sites (fig. 4) revealed isolation by distance. The relationship was positive (slope = 0.054; R² = 0.084) and statistically significant (Mantel probability P < 0.004).

Pairwise F_ST values for ND4 haploytpe data ranged from 0.000 to 0.027, yielding estimated N_m's of 17.73 to infinity. Pairwise ND4 F_ST values were not significantly different among any of the populations. Mean pairwise F_ST values for microsatellite allele frequencies ranged from 0.000 to 0.146 resulting in N_m (estimated number of migrants between populations per generation) values between 1.5 and infinity among populations (table 3). Markov chain pairwise exact tests (Goudet et al. 1996) revealed significant differences between population CTC (Coachella, CA) and all other populations including neighboring population CTM (Los Angeles, CA). Both Nebraska populations (CTNa and CTNb) were genetically distinct from populations in all other states except for CTNM in Otero City, New Mexico. Patterns in pairwise F_ST values were consistent with population structure based on AMOVA and multilocus clustering analyses.

Evolutionary Neutrality and Population Expansion

The mismatch distribution of pairwise nucleotide differences among ND4 sequences among all populations exhibited a smooth unimodal distribution characteristic of a large population expansion (Slatkin and Hudson 1991, Rogers and Harpending 1992) (fig. 5). Similar patterns were observed in the individual distributions for each population except for CTK, which exhibited a more ragged distribution (data not shown). The raggedness statistic r (Harpending 1994) was lower than expected under conditions of constant growth (P (r < 0.0336) = 0.088). Study-wide site-frequency spectra reveal an excess of singleton mutations when compared with expected frequencies under neutrality and stable population size (fig. 5). Population-specific spectra produced similar patterns (data not shown).

Fu's F_s was negative (F_s = −92.81; P < 0.0001) which occurs when an excess of rare haplotypes is present and suggests that either population expansion or genetic hitchhiking has taken place (Fu 1997). Fu and Li's D* and F* were negative as well (D* = −4.68; P < 0.02; F* = −4.38; P < 0.02), indicating an excess of recently derived haplotypes and suggesting that either population expansion or background selection has occurred (Fu and Li 1993). Strobeck's S, the probability of obtaining equal or fewer haplotypes based on gene frequency and mutation rate was very high in all populations, ranging from 0.995 to 1.000. These results are consistent with deviation from neutrality due to either selection or population expansion.

We tested for signatures of population expansion in microsatellites using the heterozygosity test assuming 2 alternative models of microsatellite evolution. Both models gave results with qualitatively similar trends. Under the SMM, assuming that only single-repeat mutations occur, expected heterozygosity based on the number of alleles and sample size (H_eq) was significantly higher than expected heterozygosity based on allele frequencies (H_e) in 10 of the 12 populations (table 4) Under the TPM, assuming that most mutations are single repeats but multistep mutations occur at a frequency of 5%, 10%, or 15%, the number of populations with significant excesses in H_eq fell to 9, 8 and 6 of 12 populations, respectively. All populations with significant H_eq > H_e exhibited this excess for at least 4 of the 5 loci.

The date of the population growth event in mutational time, or τ, was calculated to be 1.285 μ generations, where μ is the sequence's mutation rate per generation. Calculation of τ was based on the average number of pairwise differences among ND4 sequences (DnaSP v. 4.0). We used this value to estimate the approximate timing of the population expansion by the equation τ = 2μt (Slatkin and Hudson 1991), where μ was conservatively estimated as 5.7 × 10⁻⁸ substitutions/site/year (Tamura 1992) and t is time in generations. Solving for t, t = τ/2μ = 1.285/(2 × 5.7 × 10⁻⁸) = 11,271,930 generations. Assuming an average of 20–30 generations per year, we estimate that the population expansion occurred approximately 375,000–560,000 years ago. This date range should be regarded as preliminary due to the small size of the analyzed sequence and rough estimation of the gene's mutation rate.

Discussion

Sequence analysis of the ND4 mitochondrial DNA fragment from C. tarsalis populations across 4 states revealed a lack of isolation by distance, genetic structuring, and clustering of haplotypes by population or geographic region. Taken alone, ND4 results suggest that either: 1) current gene flow is very high across a large geographic range, 2) C. tarsalis experienced a selective mitochondrial sweep, or 3) C. tarsalis underwent an expansion that genetically homogenized populations after which there has been insufficient time for signatures of mitochondrial population structure to be detected. Previous allozyme data (Gimnig et al. 1999) and our microsatellite results suggest that populations are not panmictic across the collection range, refuting the 1st hypothesis. Evidence for selective mitochondrial sweeps in conjunction with invasion of maternally inherited Wolbachia endosymbionts have been observed in the congeneric Culex pipiens species complex (Rasgon, Cornel, Scott 2006), but because C. tarsalis is not infected with Wolbachia (Rasgon and Scott 2004) hypothesis 2 is not a likely explanation of our results.

The most likely explanation for our observations is that C. tarsalis underwent a population expansion some time in the relatively recent past. Statistical signatures of population expansion are detectable in both ND4 and microsatellite data. ND4 statistical parsimony analysis shows a star-shaped phylogeny (Slatkin and Hudson 1991) and a large number of rare haplotypes separated by single-nucleotide differences. H and K tests, Strobeck's S statistics, and site-frequency spectra were consistent with the presence of a rapidly growing population, whereas mismatch distributions indicated that a single or few haplotypes spread and mutated to generate several closely related haplotypes, as is often observed during population expansion (Rogers and Harpending 1992). Fu's F_s and Fu and Li's D* and F* all indicated a significant negative deviation from evolutionary neutrality. Though it is impossible to distinguish between selection and expansion as the cause of negative F_s, D*, and F* values, these statistics in combination with our microsatellite data are indicative of a deviation from equilibrium which was likely caused by population expansion.

Several populations demonstrated an excess in expected microsatellite heterozygosity based on microsatellite allele number (H_eq) when compared with expected heterozygosity based on allele frequencies (H_e), which occurs under conditions of rapid population growth (Cornuet and Luikart 1996). Heterozygote excess based on allele number was most pronounced under the SMM, though populations exhibited moderate to widespread H_eq excess and no significant deficiencies across all models tested, providing corroborative evidence of population expansion. An alternate explanation for H_eq > H_e is the influx of rare alleles from genetically distinct populations (Cornuet and Luikart 1996). This phenomenon is unlikely because multiple populations over a large geographic range demonstrated the same pattern.

Our results represent an example of how current estimates of gene flow, particularly those based on mitochondrial sequences, may be exaggerated due to past events such as range expansion. Evidence of population expansion has recently been found in many disease vectors including Triatoma infestans, Amblyomma americanum, Anopheles funestus, Anopheles darlingi, Simulium tani, Anopheles gambiae, and A. dirus (Walton et al. 2000, Donnelly et al. 2001, Pramual et al. 2005, Bargues et al. 2006, Michel et al. 2006, Mirabello and Conn 2006, Mixson et al. 2006), suggesting that it may not be an uncommon phenomenon in arthropods. Because both extensive gene flow and range expansion can produce similar patterns in mitochondrial sequence data, lack of isolation by distance and significant partitioning of variance among groups or populations may not necessarily be indicative of current panmixia. Instead, they may reflect historical events and should be interpreted carefully.

ND4 sequences suggest that the expansion of C. tarsalis in the western United States occurred between 375,000 and 560,000 years ago. Other studies of Anopheles population dynamics have suggested that range expansions in Africa and Asia may be associated with changes in agricultural practices (Walton et al. 2000, Donnelly et al. 2001). Our estimate predates human civilization and the development of agriculture in North America by well over 300,000 years, suggesting that the expansion was likely not influenced by anthropogenic forces. Instead, it is more likely to have been associated with climatic changes such as the Nebraskan or Kansan stages of Pleistocene glaciation (Dreeszen 1970) that occurred in the midwestern and western United States between 350,000 and 1 MYA.

Microsatellites in our study reflected both signatures of population expansion (nonequilibrium conditions) and isolation by distance (consistent with evolutionary equilibrium). Microsatellites have been shown to have higher rates of mutation than mitochondrial genes, ranging in insects and mammals from 10⁻³ to 10⁻⁶ mutations per locus per generation (Dietrich et al. 1992; Edwards et al. 1992; Weber and Wong 1993; Schug et al. 1997, 1998), whereas mitochondrial gene mutation rates have been estimated at 9.4 × 10⁻⁸ substitutions per 4-fold degenerate site in mammals (Brown et al. 1982) and at 5.7 × 10⁻⁸ substitutions in Drosophila species (Tamura 1992). The discrepancy between microsatellite and mitochondrial-derived genetic structure may be attributed to the fact that after range expansion rapidly mutating microsatellites have begun to differentiate among populations, whereas the more slowly mutating mitochondrial gene does not reflect the presence of genetically isolated populations.

Microsatellite-based population structure did not appear to be associated with major geographic features such as the Cascade, Sierra Nevada, or Rocky Mountain ranges. The Tehachapi Mountain range has previously been associated with genetic structure among southern California C. tarsalis populations (Gimnig et al. 1999). Although we observed that populations south of the Tehachapis (CTC and CTM) were genetically distinct from each other, the divergence did not seem to be associated with this mountain range. Interestingly, although the 2 Nebraska populations appeared to be panmictic and cluster with New Mexico, they were genetically distinct from both Colorado populations, which fall on either side of the Rocky Mountains. These results suggest that there may be a barrier to gene flow between the 2 neighboring states that is not associated with the Rockies. The arid Great Plains, which ranges across eastern Colorado and New Mexico to Iowa, may limit migration between populations in Colorado and Nebraska.

Generally, patterns of gene flow inferred from microsatellite data exhibit isolation by distance and geographical clustering but are not sufficiently explained by topographical barriers or geographic distance alone. This phenomenon may result from an insufficient number of loci, movement of mosquitoes based on agricultural practices or other environmental features beyond distance and elevation, or residual effects of range expansion. The inclusion of more loci and populations from as yet unsampled regions will likely address the first 2 possibilities. Regardless, it is still likely that microsatellites reflect only a partial return to evolutionary equilibrium after the expansion event because they exhibit signatures of both equilibrium and nonequilibrium conditions. The fact that microsatellites may underestimate current levels of genetic isolation among populations in the western United States should be taken into account when gathering and interpreting further information on the genetic structure of C. tarsalis.

Culex tarsalis is known to exhibit significant spatial and temporal heterogeneity in a variety of phenotypic traits including autogeny, feeding behavior, and vertical and horizontal arboviral transmission. Thus far, genetic and environmental determinants of such phenotypic variation, including WNV transmission, are largely unknown. Although our findings of microsatellite-based genetic structure may not reflect the full extent of reproductive isolation across the western United States, they are indicative of significant genetic differentiation among mosquito populations and thus allow for the possibility of a relationship between genetic and spatial variation in WNV vector competence. Genetic determinants of vector–pathogen interactions have been identified in several other systems by quantitative trait locus (QTL) analysis, including transovarial transmission of LaCrosse virus in Ochlerotatus triseriatus, midgut infection and escape of dengue virus in Aedes aegypti, and Plasmodium resistance in A. gambiae (Bosio et al. 2000, Graham et al. 2003, Zheng et al. 2003, Gomez-Machorro et al. 2004). Recently, Plasmodium resistance was mapped in naturally occurring populations of A. gambiae in Kenya (Menge et al. 2006) and Mali (Riehle et al. 2006). QTL analyses such as these may lay the foundation for developing biomarkers of vector susceptibility and novel targets for genetic manipulation and disease control. Future studies on genetic and phenotypic spatial variation of C. tarsalis will provide the foundation for understanding genetic determinants of viral susceptibility and disease transmission in this important arboviral vector.

Supplementary Material

The sequence data for this paper appear in GenBank under accession numbers EF125799–EF125862.

We thank W. K. Reisen for kindly providing colony and field specimens. We thank H. Lothrop, D. Lemenager, A. Moser, M. Weissman, M. DiMenna, W. L. Kramer, and M. Stromer for providing field collections. This work was funded by National Institute of Environmental Health Sciences Training Grant T32 ES07141, National Institutes of Health Grant R01AI067371 and the Johns Hopkins Malaria Research Institute. Research conducted at The W. Harry Feinstone Department of Molecular Microbiology and Immunology, Bloomberg School of Public Health, Johns Hopkins University, Baltimore, MD.

References

Author notes