|Home | About | Journals | Submit | Contact Us | Français|
We assessed the structure and latitudinal selection that might result in sensitivities to critical day-lengths that trigger diapause between Culex pipiens populations distributed along North-South and East-West axes in eastern North America. Strong population structure between Cx. p. pipiens and Cx. p. quinquefasciatus existed. Among Cx. p. pipiens, a 100-km increase in the latitudinal change resulted in an increased square root of FST by 0.002. A 100-km increase in the longitudinal change caused an increased square root of FST by 0.035. A lack of latitudinal influence on the structure between Cx. p. pipiens populations suggests a uniform signal using the 12 microsatellite markers, which might increase the risk of West Nile virus (WNV) transmission toward northern areas because of longer breeding season, extend host-seeking period, and larger population size. Northern Cx. p. pipiens may have undergone additional generations before diapause is triggered, magnifying population size when WNV amplification is peaking.
The recent spread of West Nile virus (WNV) in North America has focused attention on the major mosquito vectors that transmit this infection among birds and humans, such as Culex pipiens pipiens Linnaeus (Northern house mosquito) and Culex pipiens quinquefasciatus Say (Southern house mosquito). These mosquitoes are allopatric in North America. Cx. p. pipiens enters diapause and is typically found north of 39° N latitude, whereas Cx. p. quinquefasciatus does not undergo diapause and is restricted to areas south of 36° N latitude. Both nominal taxa and their intermediates exist between 36° and 39° N latitude.1,2 In these regions of overlap, the Cx. p. pipiens population peaks earlier in the summer season, whereas the Cx. p. quinquefasciatus population does not achieve peak abundance until the end of summer.3
The seasonality of WNV transmission corresponds closely to that of the abundance of anautogenous (i.e., requires a blood meal to reproduce) Cx. p. pipiens. Infection in birds becomes intense in the northeastern United States mainly during late August or September. Presumably, the force of WNV transmission among local birds is driven by the rising mid-summer abundance of these mosquitoes. Near the end of the breeding season, Cx. p. pipiens mosquitoes begin emerging as adults that will enter diapause and that will not seek hosts until spring.4 This diapause is triggered by a combination of short day length and low temperature perceived by the fourth larval instar and pupal stages.3,5
Day length correlates with latitude, and the number of longer days increases with increasing latitude.6 At the latitude of Boston, MA (42.4° N), the abundance of anautogenous Cx. p. pipiens crests in early August, shortly after the local day-length has fallen to < 14.25 h/day.3 This critical diel falls earlier at southern latitudes than in the north. Thus, it may be that Cx. p. pipiens populations and the force of WNV transmission may increase as latitude increases and may wane at lower latitudes because of this disparity in the length of the breeding season.5
Studies on the population structure of Cx. pipiens complex in North America are very important to potentially enable the development of new strategies for controlling the diseases transmitted by these mosquitoes. Such knowledge improves our understanding of how these populations interact with each other geographically, and whether developmental cues, such as diel length, may select for subpopulations that respond differently to changes in day-length. Such changes, if they occur, would tend to be expressed genetically as greater structuring along a latitudinal rather than longitudinal gradient. Previous studies with the use of microsatellite markers7–9 and through enzyme electrophoresis10–14 have not shown the presence or absence of such a pattern.
In this study, we analyze the effects of latitude and longitude on the population structure of Cx. pipiens s.l. mosquitoes based on variation at 12 microsatellite loci15,16 among populations along the North-South (N-S) and East-West (E-W) axes in eastern North America during the breeding season in 2004 and 2005. We explore the possibility that these mosquito populations may be structured more distinctly on a N-S than an E-W axis in eastern North America.9
We collected egg rafts of Cx. pipiens s.l. by placing black polyethylene dishpans partially filled with a bacterial hay infusion overnight in sheltered outdoor locations near homes.17 We conducted this sampling monthly from July to September 2004 and 2005 whenever weather permitted. We chose collection sites (Figure 1) to represent a broad range of latitudes and longitudes across eastern North America, so that we could compare groups along the N-S and E-W axes. Along the N-S axis, collection sites in 2004 included those at Portland, ME (ME; 43.41° N, 70.18° W); Cambridge, MA (MA; 42.2° N, 71.05° W); Lehigh County, PA (PAL; 40.5° N, 75.42° W); Fairfax County, VA (FVA; 38.51° N, 77.19° W); Richmond, VA (RVA; 37.34° N, 77.27° W), Norfolk, VA (NVA; 36.54° N, 76.18° W), and Columbia, SC (SC; 34° N, 81° W). In 2005, we established collection sites at Montreal, Québec, Canada (QUE; 45.30° N, 73.27° W) in lieu of Portland, ME, and at the zone of introgression, i.e., Memphis, TN (MTN; 35.10° N, 90° W) and Nashville, TN (NTN; 36.10° N, 86.5° W). Along the E-W axis, collection sites in 2004 included Urbana-Champaign, IL (IL; 40.07° N, 88.14° W); Columbus, OH (OH; 39.59° N, 83.03° W); Central Pennsylvania in Harrisburg and York (PAC; 40.18° N, 76.49° W); and Lehigh, PA. In 2005, we retained all sites comprising the E-W axis except for the one in Ohio.
We isolated individual egg rafts from each collection in plastic petri dishes and allowed them to hatch. We examined first-instar larvae from each family of mosquitoes under a compound microscope to morphologically differentiate Cx. pipiens s.l. from Cx. restuans18 because they both share similar larval breeding habitats. We randomly selected 10 cohorts of first-instar Cx. pipiens s.l. larvae from each collection site per month, and we shipped them to Harvard School of Public Health Laboratory. We reared the larvae further until they emerged as adults (F0 generation) under standard insectary conditions (28–30°C). We froze male and female adult specimens at −20°C and processed them for microsatellite DNA analysis. We fed representative male siblings per family with 10% sugar solution for 2 days, and later we froze them at −20°C. We examined their copulatory structures for morphologic identification of the members of Cx. pipiens s.l.
The appearance of the egg breaker on the first-instar larvae serves to separate Cx. pipiens s.l. from morphologically similar mosquito species.17,18 We examined the copulatory structures of six siblings of 2-day-old adult male mosquitoes taken from each family under a compound microscope to determine their DV/D ratio.19,20 DV is the distance from the tip of the ventral arm of the phallosome to its intersection with the dorsal arm, and D is the distance between the tips of the dorsal arms of the phallosome, measured at their intersection with the ventral arms. We performed polymerase chain reaction (PCR) as per protocol21 to confirm microscopic identification of Cx. pipiens s.l.
We processed two siblings of adult mosquitoes per family (F0 generation) each month from each collection site, thus we processed a total of 500 indi viduals of Cx. pipiens s.l. each year in 2004 and 2005 for microsatellite DNA analysis; however, we included 12 microsatellite genotypes of one sibling per family only from each population in all data analyses, except for temporal differentiation to resolve sample bias as explained below. We used the adult mosquito DNA extract of Cx. pipiens s.l. (1 μL 5 μg) as a template in the PCR reaction volume of 20 μL, which contained 2 μL of 10× PCR buffer, 0.25 μL of 10 mmol/L dNTP (premixed, 10 mmol/L each), 0.12 μL of labeled forward primer, 0.12 μL of unlabeled reverse primer, 0.1 μL of Taq DNA polymerase (5 U/μL), and 2 μL of 5× Taq enhancer (Eppendorf, Westbury, NY). We amplified microsatellite fragments for desired loci using 12 loci-specific primers. We used nine microsatellite loci that we developed15 and three loci (CxpGT4, CxpGT9, CxpGT46) from Keyghobadi and others.16 The procedure for amplifying these microsatellite loci has been described by Edillo and others.15
We measured population genetic diversity by the number of alleles and heterozygosities at each locus within populations using GENEPOP 3.4.22,23 We tested each locus separately for goodness-of-fit for Hardy-Weinberg equilibrium (HWE) using the probability test and inbreeding coefficient (FIS).22,24 We performed Fisher exact test in GENEPOP 3.4 to detect linkage disequilibrium for pairwise loci in each population. We examined the population structure by using the F-statistics, FST,25,26 calculated by using GENEPOP 3.4.23 For purposes of this study, FST is a better measure of genetic differentiation than RST, in which differentiation is primarily influenced by drift.25 We tested the significance of FST by using a log-likelihood (G) based exact test.27 We determined the temporal differentiation of FST between months in each of the breeding seasons and for each population by taking the average of FST estimates calculated from each of the five data sets randomly derived via a Monte Carlo bootstrapping method to resolve sample size bias. Each of the five data sets was composed of 10 adult mosquitoes for each month in a breeding season and for each population (i.e., N = 30 per data set with 12 microsatellite genotypes per individual) randomly and repeatedly derived via a Monte Carlo bootstrapping method from each array of siblings processed per family (F0 generation) by using custom software code written in J version 6.01 (J Software, Iverson Software Inc., Toronto, ON, Canada).
To determine the effects of the differences in latitude and longitude, controlling for year-to-year difference, on FST estimates between population pairs, we estimated the asymptotic variance-covariance matrix and the weighted least squares algorithm with the following model: √FST = β0 + β1 × Δ latitude + β2 × Δ longitude + β3 × year, in which Δ latitude is the N-S distance and Δ longitude is the E-W distance. We used the square root transformation of FST estimates to normalize them. The non-independence of both the FST estimates between population pairs, and the N-S and E-W distances required the use of this model instead of simple linear regression techniques (D. Graham, unpublished data).
We determined gene flow (Nm) from FST estimates to provide the corrected multilocus estimates of the effective number of migrants per population pairs.28 To determine the effects of the differences in latitude and longitude, controlling for year-to-year differences, on gene flow between population pairs, we used the same model: √Nm = β0 + β1 × Δ latitude + β2 × Δ longitude + β3 × year.
We identified the mosquitoes to subspecies based on the average DV/D ratio of six sibling males taken from every family that was subjected to microsatellite analysis. All sites north of the zone of introgression had Cx. p. pipiens, whereas south of the zone of introgression in SC had Cx. p. quinquefasciatus in both years. We found both Cx. p. pipiens and intermediates between Cx. pipiens and Cx. quinquefasciatus in sites within the zone of introgression, although we collected two intermediate specimens from IL, suggesting a slight expansion of the hybrid zone at this 40.07° latitude.
We measured population genetic diversity by the number of alleles and heterozygosity. The number of alleles per locus varied from 3 to 25 in 2004 (Supplemental Table 1, available online at www.ajtmh.org) and from 4 to 28 in 2005 (Supplemental Table 2, available online at www.ajtmh.org). The average observed heterozygosities for all loci within each of the populations were consistently high; in 2004, they ranged from 0.68 to 0.76, whereas those in 2005 ranged from 0.68 to 0.83.
Within each population, we calculated locus-specific departures from HWE in 2004 (Supplemental Table 1) and in 2005 (Supplemental Table 2), and we applied sequential Bonferroni procedure. In 2004, 33 of 108 tests (30.56%) for all populations (excluding the separate analysis at the zone of introgression) deviated from HWE, whereas in 2005, it was slightly lower, 29 of 108 tests (26.85%). All populations had relatively similar number of loci that departed from HWE from year-to-year, except those from SC and IL. One locus from SC population departed from HWE in 2004, and it increased to seven loci in 2005. Six loci from IL population were away from HWE in 2004, and it decreased to one locus in 2005.
Deviations from HWE generally are attributed to differential selection and non-random mating.29,30 We observed that the CxqCAG5 locus deviated from HWE across seven or eight populations (excluding the separate analysis in the zone of introgression) in 2004 and in 2005, respectively, and showed consistently large positive FIS (Supplemental Tables 1 and 2). Separate analysis of Cx. p. pipiens and the intermediates in the zone of introgression for both years resulted into a consistent deviation from HWE at CxqCAG5 locus, apparently indicating the presence of null alleles.
Alleles of different loci are not randomly associated if they occur together in individuals with a probability higher than would be expected by chance alone. There were 66 independent comparisons for each Cx. p. pipiens population, so we would expect one pair of loci (0.66%) to be significant at the 0.01 level by chance alone. After incorporating a sequential Bonferroni correction, mosquito samples from PAL and IL had all loci in linkage equilibrium in both years. One pair of loci showed linkage disequilibrium from each of the populations (P < 0.0001) in 2004 from ME (CxqGT2 and CxpGT9), MA (CxqATG9 and CxpGT9), PAC (CxqCA9 and CxqCA115), and FVA (CxqCA9 and CxqCA115), and in 2005 from QUE (CxqGT108 and CxpGT46), VA (CxqATG9 and CxpGT4 for pooled samples of Cx. p. pipiens and intermediates; CxqATG9 and CxqCA9 for Cx. p. pipiens samples only), and SC (CxqCTG10 and CxqCAG5) populations. These suggest that independent segregation of these loci in each population was within what should be expected by chance alone.
We found that differences in latitude (N-S distance) and longitude (E-W distance) were significant predictors for increases in the square root of FST estimates between Cx. p. pipiens populations at the 0.05 level. A 100-km increase in the latitudinal change resulted in an increase in the square root of FST by 0.002 (Figure 2A), and a 100-km increase in the longitudinal change resulted in an increase in the square root of FST by 0.035 (Figure 2B). We found greater genetic structuring between Cx. p. pipiens and Cx. p. quinquefasciatus populations in both years (Figure 2C). In 2004, the FST estimate between these populations at the N-S extremes of our samples was 0.108, greatly exceeding 0.012, which represented the E-W extremes (Table 1). Likewise in 2005, the FST estimate between these populations at the N-S extremes was 0.075, exceeding 0.05, which represented the E-W extremes (Table 1). We found a significant difference of FST estimates between different population pairs in 2004 and 2005 (P < 0.01; Figure 2D). In the hybrid zone, we found little genetic differentiation between those neighboring populations of Cx. pipiens from FVA and the intermediates from VA in both years (FST = 0.016 in 2004, FST = 0.019 in 2005; P < 0.01) and those intermediates from VA and TN (FST = 0.047; P < 0.01) in 2005 (Table 1). Genetic differentiation became moderate as geographic distance increased between those populations of Cx. pipiens from FVA and the intermediates from TN (FST = 0.061, P < 0.01) in 2005 (Table 1).
Because the CxqCAG5 locus exhibited the highest heterozygote deficits, and null alleles were likely present that might have caused a biased estimate, we excluded this locus and reanalyzed FST estimates. Patterns of population differentiation did not differ when the CxqCAG5 locus was excluded from the re-analysis of FST estimates The mean of FST estimates in 2004 (Table 1) was 0.028, whereas the mean of the adjusted FST estimates was 0.027. Likewise, the mean FST estimate in 2005 (Table 1) was 0.037, whereas that of the re-analyzed FST estimates was 0.036.
We calculated the multilocus estimates of the effective number of migrants per generation between population pairs indirectly from FST estimates (Table 1). Estimates of gene flow were lower between Cx. p. quinquefasciatus and Cx. p. pipiens populations (mean ± SE: 2.78 ± 0.25 in 2004 and 2.63 ± 0.12 in 2005) than between population pairs of Cx. p. pipiens (mean ± SE: 4.20 ± 0.15 in 2004 and 3.12 ± 0.10 in 2005), indicating the existence of greater genetic structuring between sibling species.
Between Cx. p. pipiens population pairs, we found that the latitudinal and the longitudinal changes between collection sites remained significant predictors for gene flow at the 0.05 level. A 100-km increase in the N-S distance resulted in a decrease in the square root of Nm by 0.179, and a 100-km increase in the E-W distance resulted in a decrease in the square root of Nm by 0.007. The difference in gene flow across years remained significant (P < 0.01).
We determined whether there was temporal differentiation of FST estimates between months of the breeding season. This was performed by taking the average of FST estimates calculated from each of the five data sets derived by Monte Carlo bootstrapping repeated sampling method taken from two siblings per family of 10 in each month for each population (Table 2). Six of the eight Cx. p. pipiens populations sampled in July 2004 were moderately different26 from those sampled in September 2004; this trend was not apparent in 2005. Cx. p. quinquefasciatus population in SC did not generally differ between months of each year's breeding season.
In this study, we found that when Cx. p. quinquefasciatus population from SC was included in the pairwise comparisons with Cx. p. pipiens populations from each of the collection sites, as expected a much greater distinction of FST estimates was observed between latitudinal pairs rather than longitudinal pairs in both 2004 and 2005. The greater gradient observed in FST estimates may result from introgression northward from the hybrid zone and into the mid-Atlantic and New England regions. Flight of Cx. p. quinquefasciatus was documented as far as 12.6 km.31 Dispersal of Cx. p. quinquefasciatus at distances greater than their average flight may be associated with human activity, such as long distance transport by commercial trucks.32,33 Concordant with our FST analyses, estimates of gene flow between these mosquito populations also showed the existence of greater genetic structuring. Barriers to gene flow exist between Cx. p. pipiens and Cx. p. quinquefasciatus because of their distinct latitudinal distribution.3 The inability of Cx. p. quinquefasciatus to hibernate seems to limit their northward expansion into the latitudinal range of Cx. p. pipiens where winters are severe.34
The hypothesis that these Cx. p. pipiens populations may be structured more distinctly on a N-S than an E-W axis in eastern North America was not shown. If latitudinal selection for diel response was solely responsible for these differences, we would have expected distinct population structure between Cx. p. pipiens populations at varying latitudinal distributions. Our results apparently indicated that with the very minimal effect of latitudinal selection for diel response, one would expect northern populations of these mosquito vectors to have fewer restraints on their population growth because the environmental cue (14.25 hours of daylight) would be reached later in the breeding season at higher latitudes, allowing additional generations of host-seeking mosquitoes to be generated before fourth-instar larvae and pupae are triggered to enter an overwintering diapause in which they do not seek hosts, and a metabolic switch from blood feeding to sugar gluttony follows.35 The lack of a clear trend of latitudinal genetic structuring among Cx. p. pipiens populations lends support toward a uniform signal among North American Cx. p. pipiens mosquitoes, although we cannot assume that the 12 microsatellite markers used reflect the overall genome. Such an effect would increase risk of WNV transmission towards northern areas because of the longer breeding season, an extended period for host seeking, and the larger size of the maximum population of mosquitoes that would result from this combination of factors.
The latitudinal and the longitudinal changes between collection sites remained significant predictors for gene flow among Cx. p. pipiens mosquitoes. Barriers to gene flow in Cx. p. pipiens populations may include both by geographical distance and topography, especially in species with low vagility. Mountainous terrain (e.g., Appalachian Mountains, Notre Dame Mountains) and river (e.g., St. Lawrence River in Québec) are particularly important barriers; however, dispersal may occur in mosquitoes during wind storms or through transportation by humans.32,33
Six of the eight Cx. p. pipiens populations sampled in the early and later period of the breeding season in 2004 showed moderate temporal differentiation of FST, which may be associated with their ability to hibernate, but not in 2005 when temperature was warmer. Cx. p. quinquefasciatus population sampled in SC between months of each breeding season did not differ, which might be apparently consistent with its inability to hibernate. It may be that Cx. p. pipiens that were highly fit to hibernate (i.e., ability to respond to critical diel and temperature) were selected for and might have introduced a differential factor in the later part of the breeding season in 2004.
Temperature records taken from the nearest weather stations of the collection sites showed that the mean (±SE) monthly temperature during our mosquito collections from July to September 2004 along the E-W transects was 21.41 ± 0.25°C, whereas along the N-S axis it was 23.05 ± 0.95°C. Mean (±SE) monthly temperature from July until September 2005 was higher than that in 2004; along the E-W sites, it was 23.52 ± 0.23°C, whereas along the N-S axis, it was 25 ± 0.99°C. Most likely, warmer temperatures may have influenced the genetic similarity of Cx. p. pipiens populations sampled between months of the breeding season in 2005. Previous reports 6,36 showed that Cx. p. pipiens in the laboratory after feeding at short photoperiod (12 hours) and warmer temperatures (25°C) were able to terminate diapause, whereas cold temperatures (15°C) increased the incidence of diapause. Moreover, warmer temperature also limits the number of larval development sites that may act as important barrier to dispersal of these mosquitoes. The varying amount of shading caused by vegetation cover and tree canopy in our collection sites may have also influenced the exposure of mosquitoes to latitudinal day length.
In conclusion, Cx. p. pipiens and Cx. p. quinquefasciatus populations have strong population structure; however, Cx. p. pipiens populations are not structured more distinctly on a N-S than an E-W axis in the eastern North America, although latitudinal and longitudinal changes between their collection sites are significant predictors for gene flow. A lack of latitudinal effect on the structure between Cx. p. pipiens populations suggests a uniform signal considering the 12 microsatellite markers used. A latitude change may have exerted a stronger effect on the length of the breeding season. Thus, northern Cx. p. pipiens may have undergone additional generations before diapause is triggered, magnifying population size when WNV amplification is peaking, which may increase the risk of WNV transmission towards northern areas in eastern North America.
The authors thank N. Grefe, S. Hutter, M. Holman, R. Robich, M. Tabibi, N. Whitehurst, and M. Reddy for their help in obtaining and/or rearing mosquito samples.
Financial support: This research was supported by Grant RO1A1 52284 from the National Institute of Allergy and Infectious Diseases and funds provided by the Centers for Disease Control and Prevention under Grant RO1 A1 44064 from the National Institutes of Health to the late A.S. Statistical analyses of latitude and longitude were supported by Grant T32 AI007535 from the National Institute of Allergy and Infectious Disease to J.M. and Grant RO1 EB 006195 from the National Institute of Biomedical Imaging and Bioengineering to M.P. This work is dedicated in memory of Andrew Spielman, a famous scientist and professor who died after completing the molecular work by F.E.E.
Frances Edillo, Harvard School of Public Health, Department of Immunology and Infectious Diseases, Boston, MA 02115.
Anthony Kiszewski, Harvard School of Public Health, Department of Immunology and Infectious Diseases, Boston, MA 02115.
Justin Manjourides, Harvard School of Public Health, Department of Biostatistics, Boston, MA 02115.
Marcello Pagano, Harvard School of Public Health, Department of Biostatistics, Boston, MA 02115.
Michael Hutchinson, Department of Environmental Protection, Harrisburg, PA 17105.
Andrew Kyle, Department of Environmental Protection, Harrisburg, PA 17105.
Jorge Arias, Fairfax Department of Health, Fairfax, VA 22030.
David Gaines, Virginia Department of Health, Office of Epidemiology, Richmond, VA 23218.
Richard Lampman, Illinois Natural History Survey, Division for Biodiversity and Ecological Entomology, Champaign, IL 61820.
Robert Novak, Division of Infectious Diseases, University of Alabama at Birmingham School of Medicine, Birmingham, AL 35294.
Ivo Foppa, Department of Epidemiology, Tulane School of Public Health and Tropical Medicine, New Orleans, LA 70112.
Charles Lubelcyzk, Maine Medical Center Research Institute, Vector-borne Disease Laboratory, South Portland, ME 04106.
Robert Smith, Maine Medical Center Research Institute, Vector-borne Disease Laboratory, South Portland, ME 04106.
Abelardo Moncayo, Tennessee Department of Health Communicable and Environmental Disease Services, Vector-borne Disease Section, Nashville, TN 37216.