- Open Access
High genetic diversity in the Culex pipiens complex from a West Nile Virus epidemic area in Southern Europe
Parasites & Vectors volume 9, Article number: 150 (2016)
The Culex pipiens complex includes the most widespread mosquito species in the world. Cx. pipiens is the primary vector of the West Nile Virus (WNV) in Europe and North America. Cases of WNV have been recorded in Italy since 1998. In particular, wet areas along the Po River are considered some of the most WNV affected areas in Italy. Here, we analyzed the genetic structure of ten Cx. pipiens populations collected in the last part of the Po River including the Delta area.
We assessed the genetic variability of two mitochondrial markers, cytochrome oxidase 1 (COI) and 2 (COII), for a total of 1200 bp, and one nuclear marker, a fragment of acetylcholinesterase-2 (ace-2), 502 bp long. The effect of the landscape features was evaluated comparing haplotype and nucleotide diversity with the landscape composition.
The analysis showed a high genetic diversity in both COI and COII gene fragments mainly shared by the populations in the Delta area. The COI-COII network showed that the set of haplotypes found was grouped into three main supported lineages with the higher genetic variability gathered in two of the three lineages. By contrast, ace-2 fragment did not show the same differentiation, displaying alleles grouped in a single clade. Finally, a positive correlation between mitochondrial diversity and natural wetland areas was found.
The high mitochondrial genetic diversity found in Cx. pipiens populations from the Po River Delta contrasts with the low variability of inland populations. The different patterns of genetic diversity found comparing mitochondrial and nuclear markers could be explained by factors such as differences in effective population size between markers, sex biased dispersal or lower fitness of dispersing females. Moreover, the correlation between genetic diversity and wetland areas is consistent with ecosystem stability and lack of insecticide pressure characteristic of this habitat. The mtDNA polymorphism found in the Po River Delta is even more interesting due to possible linkages between the mitochondrial lineages and different biting behaviors of the mosquitoes influencing their vector ability of arboviral infections.
The Culex pipiens complex includes two of the most widespread mosquito species in the world, Cx. pipiens L. (the northern house mosquito) in temperate climates and Cx. quinquefasciatus Say (the southern house mosquito) in tropical and subtropical climates. Although still under debate, two Palaearctic forms or biotypes are usually recognized within Cx. pipiens, the pipiens and the molestus forms [1, 2]. The two forms are morphologically very similar, nonetheless the two forms exhibit important behavioral and physiological differences, such as host choice and breeding sites [3, 4].
The mosquitoes of the Cx. pipiens complex are considered to be the primary vectors of the West Nile Virus (WNV) in Europe and North America [5, 6]. The cycle of infection involves birds and mosquitoes. In particular, migratory birds are instrumental in the introduction of the virus to temperate areas of Eurasia during spring migrations . The mosquitoes belonging to the Cx. pipiens complex feed on both avian and mammalian hosts. This behavior assigns to mosquitoes of the complex the role of bridge-vector for the transmission of the WNV to birds and mammals, including humans [5, 8]. In Italy, after the first report in 1998, WNV became endemic in the North-East after 2008 [9, 10] and since then, in this area, it has been constantly detected in humans, in animals and in the vector mosquitoes . In particular, the areas considered the most affected by WNV outbreak are the wet areas surrounding the Po River, in the Veneto and Emilia Romagna regions [12, 13]. Moreover, the recent finding of simultaneous circulation of two WNV lineages in this area confirms North-Eastern Italy as a high-risk area for WNV emergence .
This scenario has drawn attention to mosquitoes of the Cx. pipiens complex inhabiting the last part of the Po River. In this area, urban zones are mostly alternated with large cultivated areas, rice paddies and natural wetlands. Here, during summer, Cx. pipiens achieves high population densities, creating much annoyance. Despite its epidemiological relevance, population genetic data for the Cx. pipiens complex in this WNV outbreak area are still lacking. We consider it particularly important to understand the population structure and the distribution patterns of Cx. pipiens in a WNV epidemic area. Genetic variability in local Cx. pipiens populations could underlie biological differences in biting and flight behavior, and in spatial and temporal distribution, with effects on the epidemiological patterns of the WNV infection.
Up to date, some surveys about genetic variability for Cx. pipiens complex have been carried out in several European countries using allozymes and mitochondrial markers [15–17]. These studies showed a low genetic variability in most of the populations analyzed compared to the closely related species Cx. torrentium Martini, morphologically very similar to Cx. pipiens and occurring in sympatry with this species throughout Europe and some parts of Asia . Other genetic studies focused on the sympatric occurrence of the molestus and pipiens forms, underlining the presence of asymmetric introgression between the two forms and its implications for the WNV transmission [18, 19].
In this study we inferred the genetic structure of the Cx. pipiens complex along the last part of the Po River using both mitochondrial and nuclear markers. In detail, we studied the genetic diversity of cytochrome c oxidase subunit 1 (COI) and the cytochrome c oxidase subunit 2 (COII) in the populations considered. The cytochrome c oxidase fragment has been frequently used to study mosquito species complexes as well as to compare phylogeographic patterns within closely related taxa e.g. [20–22] due to the high number of copies per cell of the mitochondrial genome, its maternal inheritance and its rare recombination rate. Moreover, in order to investigate the diversity of Culex mosquitoes and to identify possible cryptic taxa, we analyzed the acetylcholinesterase-2 gene (ace-2), a protein-coding gene often used for discrimination within the Culex complex species [21, 23, 24]. Finally, we assessed the possible effect of landscape composition on the genetic diversity. Landscape composition is known to influence mosquito density and diversity [25, 26] but, to our knowledge, studies associating genetic variability and landscape composition are still lacking for Cx. pipiens. This information could help to better identify landscape features that may affect the population structure of mosquitoes, mainly in a high-risk area for arbovirus outbreaks.
Mosquito collection was carried out in the Province of Rovigo (Veneto Region-North Eastern Italy), an area of about 1,500 km2. The Rovigo Province consists of a plain located between the lower portion of the rivers Adige and Po, the two longest Italian rivers. The study region is a wide lowland area with large zones below sea level and a network of embanked rivers and channels.
The eastern part, facing the Adriatic Sea, terminates in the Po River Delta, which covers an area of about 380 km2. The Po Delta area is characterized by intensive agriculture (rice paddies, maize, wheat) due to its flat topography and the abundance of surface waters for irrigation. Moreover, lagoons and wooded wetland patches, with high biodiversity, are frequent.
The Po Delta is an area located at the crossroads of the most important bird migration routes connecting Europe, the Mediterranean Basin and Africa; moreover, several equine stables are present [27–30]. This area, with more than 300 species of birds, is considered the most important ornithological area in Italy and one of the most relevant in Europe. For these reasons, since 1997 the Po Delta wetlands have been protected by the institution of the Regional Po Delta Park.
Mosquito specimens were collected twice a month during summer and autumn in 2012. A total of ten geo-referenced sampling sites were chosen for insect collection: nine sites located in the Rovigo Province and one, used as an out-group, located in the Venice Province (Table 1).
Mosquitoes were collected as larvae or pupae from drainage or natural ditches using the dipper method. At each visit, larvae and pupae were sampled with at least five dips for each collection site, using half a liter of water, at randomly located points taken at about 30 m intervals along the ditches to avoid sampling bias. After collection, the samples were transferred to the laboratory, reared to the adult stage, identified and then stored in 80 % ethanol before being analyzed.
No ethical approval is required for the experimental methods used in this study.
The morphological distinction of females of Cx. pipiens from other mosquitoes belonging to the genus Culex is almost impossible . Thus, only male specimens belonging to the Cx. pipiens complex were used to carry out molecular analysis. Being generally accepted that the main diagnostic character of adult mosquitoes is the structure of male genitalia , a phallosome analysis was carried out before the genetic analysis in order to identify the samples morphologically. Phallosome was extracted following the protocol described by  with some modifications. Mounting of phallosome on microscope slides and further morphological examination of the specimens were made using the identification keys of the Italian Culicidae adults . Observations and measurements were taken with the help of ocular micrometer at 40X of the microscope with phase contrast.
Legs from each mosquito were removed for molecular analysis, while the rest of the body was preserved for future morphological studies.
DNA was extracted according to a previously described salting-out protocol . A region of the mitochondrial DNA corresponding to a fragment of the cytochrome c oxidase subunit 1 (COI), tRNA-Leu, and the cytochrome c oxidase subunit 2 (COII) was amplified using the universal primer pairs LCO-1490/HCO-2198  and C1-J-2195/TK-N-3796 . Amplifications were performed in 20 μl reactions (1x PCR Go Taq Flexi buffer - Promega, 2.5 mM MgCl2, 0.1 mM dNTPs, 0.5 μM for each primer, 0.5 U of Taq polymerase - Promega, 2 μl DNA template). Thermal cycling conditions for the fragment including the 5’ region of the COI were 5 min at 96 °C followed by 4 cycles of 96 °C for 1 min, 47 °C for 1 min, and 72 °C for 1 min, and other 35 cycles of 96 °C for 1 min, 50 °C for 1 min, and 72 °C for 1 min, with a final extension of 72 °C for 5 min; for the fragment including the 5’ region of the COII, PCR conditions were 5 min at 96 °C followed by 35 cycles with a denaturation step of 96 °C for 1 min, an annealing step of 50 °C for 1 min, and an extension step of 72 °C for 1 min, with a final extension of 72 °C for 5 min.
A fragment of the acetylcholinesterase 2 gene (ace-2) was amplified using the primer pairs F1457/B1246 and Acepip/B1246 . The cycling program of the PCR consisted of a first step at 95 °C for 5 min followed by 35 cycles with a denaturation step of 95 °C for 30 s, an annealing step of 60 °C for 30 s, and an extension step of 72 °C for 1 min, with a final extension of 72 °C for 5 min.
PCR products were checked through electrophoresis on 1.0 % agarose gels stained with SYBR® (Invitrogen) and purified using exonuclease and antarctic phosphatase (GE Healthcare) before sequencing.
Sequencing was performed at the BMR Genomics Service (Padova, Italy) on automated DNA sequencers employing for the mitochondrial markers, the primers LCO-1490 and TK-N-3796 and, for the ace-2 fragment, the primers Acepip and B1246.
Sequences were edited and aligned using MEGA 6.06 . A GenBank BLAST analysis of the sequences obtained was run through the NCBI website (www.ncbi.nlm.nih.gov) to assess the identity of the sequences. For the COI fragment, the integrated bioinformatics platform Barcode of Life Data (BOLD) System database (www.barcodinglife.org) was also used to identify sequences.
A partition homogeneity test was performed for the COI and COII fragments using PAUP* 4.0b10 . Haplotype and nucleotide diversity as well as the pairwise genetic distances between populations using a Kimura 2-parameters model were calculated with Arlequin 3.5 . The presence of population differentiation was also tested and accomplished by conducting exact tests of population differentiation with 100,000 steps in Markov chain, with 10,000 dememorization steps.
Phylogenetic relationships among sequences of the COI-COII data set were estimated using the approximate maximum-likelihood (ML) analysis. A GTR + I + G model implemented in PHYML 2.4.4 software  was applied, with neighbor-joining starting tree and 100 bootstrap replications. A haplotype parsimony network of the COI-COII dataset with a probability cut-off at 95 % was reconstructed using the software TCS 1.21 .
The demographic history was inferred within each lineage through the Tajima’s D test  and mismatch distributions of the pairwise genetic differences  using Arlequin 3.5. Populations at demographic equilibrium or decreasing in their size should provide significant positive D values with a multimodal distribution of pairwise differences, whereas populations that have undergone through a sudden demographic expansion usually show significant negative D values with a unimodal distribution [43, 44]. The sudden expansion model was tested through the analysis of the sum of square deviations (SSD) and the raggedness index (r) representing the modality of the distribution obtaining the corresponding P values with a parametric bootstrap approach (10,000 replicates).
The ace-2 haplotypes were reconstructed using a Bayesian statistical method implemented in PHASE 2.1 [45, 46]. Two alleles for every mosquito specimen were reported by this software. Individuals with haplotypes that could not be phased with a probability > 90 % were removed from the data set. In addition, we tested for recombination in the remaining dataset using the methods implemented in the RDP4 software  with default settings. Moreover, the heterozygote samples bearing haplotypes never found in homozygosis were cloned using the pGEM-T Easy vector (Promega, Madison, WI) following the manufacturer’s protocols with at least 5 clones per reaction, in order to confirm the haplotypes found by PHASE 2.1. Gene and nucleotide diversity were obtained using Arlequin 3.5.
The phylogenetic relationships among the different ace-2 haplotypes were calculated using two methods: the ML and the Bayesian inference methods. The ML tree was obtained following the same methodology reported above for the mitochondrial markers. For the BI analysis of the same data set Mr Bayes 3.1.2  was used. Two independent iterations were run for 5,000,000 generations and sampled every 100 generations. The 50 % majority rule consensus tree and Bayesian posterior probability of support were obtained discarding the first 25 % of sampled generations.
Links between genetic diversity and landscape features
In order to assess possible correlations between landscape features and genetic diversity for both nuclear and mitochondrial markers, two variables were considered: landscape composition and air distances from the sampling sites to the closest main natural wetlands.
The correlation between genetic diversity and landscape composition was evaluated comparing haplotype and nucleotide diversity with different land cover categories, quantified within a 1,000 m radius, around the sampling site coordinates using detailed land-use maps (Quadro conoscitivo – Regione Veneto -www.regione.veneto.it/web/ambiente-e-territorio/quadro-conoscitivo) in QGIS software (QGIS Development Team, 2014. QGIS Geographic Information System. Open Source Geospatial Foundation Project, http://qgis.osgeo.org).
For each sampling site, we quantified four land cover categories: 1) urban areas, including residential, commercial and industrial areas; 2) rural areas, comprising arable lands with annual crops; 3) tree crops, including poplar plantations and orchards, and 4) wetlands, comprising surface waters and natural wetland areas. The response variables were haplotype and nucleotide diversity, tested separately across the ten sampling sites. The four land cover categories were tested in order to detect possible collinearity between them. The land cover categories that were not correlated were included in the models as explanatory variables.
The correlation between the distances from the sampling sites to the closest main natural wetlands and genetic variability (haplotype and nucleotide diversity) was tested through a regression analysis. We considered as main natural wetlands those with an area at least 10 km2. These areas were located along the coast and were included in the Delta area or in the Venice Lagoon. All statistical analyses were carried out using R software .
Genetic diversity within COI and COII markers
A total of 158 of the 162 adult mosquitoes analyzed, morphologically identified as members of the Cx. pipiens complex and representing ten populations (Table 1), were successfully sequenced. A fragment of 571 bp of COI and an additional fragment of 629 bp of the 5’ part of COII gene was obtained for all the individuals. Alignment of sequences showed a total of 19 and 21 polymorphic sites for COI and COII respectively.
All sequences were translated with Transeq (EMBOSS: http://www.ebi.ac.uk/Tools/st/emboss_transeq/) to exclude the presence of nuclear mitochondrial pseudogenes. All of the substitutions, except one, resulted to be synonymous thus not affecting the amino acid sequence. In haplotype 14 of COII fragment, the replacement of a G by an A in pos 22 resulted in the amino acid change G8S. Sequences were deposited in the NCBI database with the GenBank accession numbers reported in Additional file 1: Table S1.
A comparison with GenBank and BoldSystem databases showed for the COI a similarity > 99 % with members of the Cx. pipiens complex. In particular, 125 sequences showed a similarity of 100 % with Cx. p. molestus [17, 50], two sequences corresponded to Cx. p. pipiens [50, 51] while the remaining 31 sequences showed only a 99 % similarity with both sequences of Cx. pipiens (undefined form)  and sequences of Cx. quinquefasciatus [17, 52–54]. On the other hand, COII sequences showed homology only with Cx. quinquefasciatus (99–100 %) (i.e. ) as for this gene no sequences of Cx. pipiens were present in GenBank at the time of writing (July 2015).
Partition homogeneity test confirmed that COI and COII fragments bear a homogeneous signal (P = 0.28), allowing data to be pooled for further analyses. Diversity indexes ranged between 0 and 0.92 for the haplotype diversity (H) and between 0 and 0.77 % for the nucleotide diversity (π). Populations S07 and S08 showed the significantly highest H values (H = 0.92 and 0.6 respectively) while S02 and S05 showed the lowest (H = 0). Regarding nucleotide diversity, population S07 displayed the highest value together with populations S06, S08 and S09 while populations S02, S05 and S10 showed the lowest variability. The distribution of haplotypes among the ten populations and other summary statistics are shown in Table 2.
Pairwise genetic distances (F st ) between analyzed populations showed significant F st values only for population S07 when compared to all the other populations except for populations S06 and S08 (Table 3). On the other hand, exact test of population differentiation was significant (P < 0.05) for all the pairwise comparisons, indicating a non-random distribution of haplotypes among populations.
Relationships between mitochondrial haplotypes
A statistical parsimony network obtained joining COI and COII sequences included 28 haplotypes. Three major lineages were identified and reported in the haplotype network (Fig. 1a) through bootstrap supports calculated on an ML tree (Fig. 1b and Additional file 2: Figure S1). Lineage 1 consisted in four haplotypes including three unique haplotypes (B1, B2, and C6) and the most frequent haplotype (A1) present in approximately 79 % (n = 125) of the mosquitoes analyzed. This haplotype was found in all the populations analyzed and in particular it was exclusive in populations S02 and S05 (Table 2; Fig. 1c). Lineage 2 grouped together eight rare haplotypes mainly present in the populations close to the Delta area (S03, S06, S07, S08, and S09) (Fig. 1c). This clade was separated from lineage 1 by nine mutational steps with missing intermediate haplotypes. Moreover, within the clade, seven haplotypes were separated from each other by several missing intermediate haplotypes (i.e. 11 mutational steps separated the haplotype Q15 from M18). Lineage 3 included 16 rare haplotypes with 12 of them connected in a star-shape manner with an inner haplotype (G7) shared by only two samples. Only three of the haplotypes were represented by two sequences each, while the remaining 13 were represented by a single sequence. The haplotypes of this clade were present at higher frequency in populations of the Delta area (S06, S07, S08 and S09). Interestingly, more than 50 % of the samples of population S07 bore haplotypes of this clade (Fig. 1c).
Tajima’s D test was applied in each of the three lineages identified in order to check for past demographic events. The null hypothesis of neutrality was rejected only in lineage 1 (D = -2.07; P < 0.01) suggesting a past population expansion.
The analysis of mismatch distributions for the three lineages revealed differences among their demographic histories. Mismatch distribution of lineage 1 was unimodal and characterized by small observed mean (0.17 ± 0.01 S.E.M.). In contrast, the mismatch distributions in lineage 2 and lineage 3 were broadly multimodal, with higher mismatch observed means of 5.82 ± 0.39 S.E.M and 4.47 ± 0.21 S.E.M, respectively as expected under a model of relative constant population size (Fig. 2). Moreover, lineage 1 and lineage 2 had SSD and raggedness index values that did not reject a sudden expansion model. Conversely, lineage 3, although not significant for the SSD, showed a significant raggedness index (r = 0.17, P < 0.01) suggesting that this lineage may have had a relatively stable population size over time.
Genetic diversity within ace-2 marker region
A subset of 74 individuals of the 158 mosquitoes previously sequenced for the mitochondrial DNA was used for the nuclear analysis. In particular, all the samples showing rare mitochondrial haplotypes and a representative part of the samples with the most frequent haplotype (A1) in the combined data set (COI-COII) were chosen for these analyses. Sequencing of the ace-2 gene yielded a 502 bp fragment with the exon 3 spanning between positions 313–502. Through the analysis carried out with the software PHASE, 12 samples were removed since their alleles could not be correctly assigned (P < 90 %). Thus a total of 16 unique alleles of ace-2 gene were identified in the remaining 62 samples. Cloning confirmed the haplotypes found by PHASE whereas the RDP4 software did not detect any statistically significant recombinant among the ace-2 haplotypes. Variable sites were mainly found in intron 2 (15 polymorphic sites out of 312 bp) while only one nucleotide substitution out of 190 bp was found in exon 3, that did not change the inferred amino acid sequence. The sequences were deposited in the NCBI database with the GenBank accession numbers reported in Additional file 1: Table S1. A BLAST analysis of these sequences showed a similarity between 99–100 % with Cx. pipiens [5, 24].
Gene diversity related to ace-2 was high in all the populations analyzed, with values ranging between 0.61 (population S01) and 0.92 (populations S07). Nucleotide diversity (π) showed values between 0.33 and 0.42 for all the populations except for populations S06 and S07, which showed the highest values (π = 0.52 % and π = 0.64 %, respectively) (Table 2).
Phylogenetic analysis conducted on the ace-2 alleles found in this study (Additional file 3: Table S2), showed that all the alleles were grouped in one highly supported cluster (Bp/Pp = 100/100) including the ace-2 haplotypes retrieved from GenBank corresponding to Cx. p. pipiens [FJ948081, JF430595, JF501651] and Cx. p. molestus [AB294405] (Fig. 3). No internal structure was present within this cluster with the exception of a subclade grouping two haplotypes supported by both maximum likelihood and Bayesian inference (Bp/Pp = 81/97).
Links between genetic diversity and landscape features
The landscape composition quantified within a 1,000 m radius around on the sampling site coordinates showed a high percentage of urban and rural areas in all the sites, with an average percentage of 19 and 66 % respectively. The other two categories (tree crops and wetlands) presented a low coverage in almost all the sites, with the exception of S07 and S08 where a high percentage of the latter category was observed (42 and 33 %, respectively) (Table 1).
Since an initial analysis showed a significant inverse correlation between rural and urban land categories only the latter was included in the models. The analyses for the mitochondrial diversity showed that only wetland areas category was significantly and positively correlated with Cx. pipiens complex genetic diversity, for both the nucleotide (P < 0.01) and haplotype (P = 0.0 < 1) diversity. Conversely, no significant correlations were found when considering the nuclear marker diversity.
Finally, we found a significant inverse effect of the distance from the sampling sites to the main natural wetland areas on the mitochondrial haplotype (R2 = 0.83; P < 0.001) and nucleotide (R2 = 0.77; P < 0.001) diversity of each population, that decrease from the main natural wetlands close to the coast to the inner areas in the mainland (Fig. 4a-c). On the other hand, nuclear marker diversity did not show any significant correlation with the distance from the main natural wetland areas (Fig. 4b-d).
Cx. pipiens populations in the WNV outbreak area along the last part of the Po River showed a high genetic diversity in both mitochondrial and nuclear markers. In particular, for the COI and COII gene fragments, the diversity was mainly concentrated in four populations occurring in the lowlands in the Delta area (S06, S07, S08, S09), with one population (S07) showing also a significant genetic distance from the other populations. A different pattern of variability was found in the ace-2 marker, with all the populations showing homogenous genetic variability values with population S07 displaying the highest nucleotide diversity value (Table 3). This outcome is even more interesting considering that in this study we examined the population genetic structure across a small geographical area.
The high variability found in the Delta area, with the presence of several rare haplotypes, contrasts with the low genetic diversity found in other areas in previous studies. A lowest genetic variability for COI was found for example in most of the populations analyzed in Germany [16, 56] and Russia . Similarly, a low intraspecific variability has been observed for the COII in populations collected around the world and reared in laboratory . This difference in genetic variability could be due to the fact that most of these surveys on the Cx. pipiens complex have been mainly focused on urban and suburban environments and do not include natural habitats. In general, surveillance efforts for potential disease-carrying mosquitoes have necessarily focused on high human population density areas with a greater risk of disease transmission to human populations .
The haplotypes found in this study are structured in three distinct lineages showing differences in their composition (Fig. 1a). Interestingly, the high genetic variability was condensed into two lineages (lineage 2 and 3) mainly linked to the populations of the Delta area. These two lineages were characterized by the presence of several rare haplotypes without any high frequency haplotype. This pattern is consistent with what is expected for populations that have not experienced past demographic expansions as highlighted by both the neutrality test and the multimodality of the mismatch distribution curves of these two lineages . Moreover, clues of past bottlenecks may be inferred for lineage 2 due to the presence of many peaks at larger values of pairwise nucleotide differences, usually generated after drastic contractions of population size [43, 60]. Due to the fine-scale sampling of this study it is premature to speculate on the causes that led to these past demographic events in the two lineages and further phylogeographic analyses would be needed extending the research on a wider range. In contrast with the high variability of these two lineages a reduced mitochondrial variability was found in lineage 1. This lineage was characterized by the presence of one haplotype shared by high number of individuals and only other three rare haplotypes. This may be related to a recent population expansion as suggested by both the neutrality test and the unimodal mismatch distribution. Interestingly, the most common haplotype found in our sampling and belonging to the lineage 1, is reported to be associated with Cx. pipiens form molestus, mainly linked to human-made environments [50, 61]. The presence of a high frequency haplotype could be also the result of selective sweeps due to symbiotic bacteria. Possible influence of Wolbachia producing selective sweeps has been suggested by several authors in order to justify the reduced haplotype diversity in Cx. pipiens [17, 22, 58] and Cx. quinquefasciatus [62, 63].
The occurrence of three mitochondrial lineages may underlie the existence of cryptic species within Cx. pipiens as it has already been suggested for the strong genetic structure found in Cx. torrentium . Noteworthy, some of our haplotypes of lineages 2 and 3 resulted similar to haplotypes retrieved from GenBank of the southern house mosquito, Cx. quinquefasciatus. However, the morphological analysis excluded that the specimens with these haplotypes could belong to Cx. quinquefasciatus as they presented the divergent dorsal arms of the male genitalia characteristic of Cx. pipiens in contrast to the parallel dorsal arms of Cx. quinquefasciatus . Moreover, sequences obtained for the ace-2, a marker successfully used to distinguish the Cx. pipiens and Cx. quinquefasciatus , confirmed the morphological results as they showed a strong differentiation from the Cx. quinquefasciatus sequences deposited in GenBank (Fig. 3).
The high variability found for ace-2 in the populations of Cx. pipiens analyzed is consistent with the diversity found in other studies for this nuclear marker [64, 65] and is largely due to the higher mutation rate of the intron region [66–68]. On the other hand, in contrast to the well differentiated mitochondrial haplotype structure, all the ace-2 sequences grouped in a supported cluster without a strong internal subdivision.
Several studies report discordant patterns between mtDNA and nuclear markers in animals . These incongruences could be explained by diverse factors such as differences in effective population sizes between the two markers, sex-biased dispersal and selection. Discordances between mitochondrial and nuclear markers are expected as mitochondrial genome is haploid and usually uniparentally inherited leading to a fourfold smaller effective population size [70, 71]. Mitochondrial DNA can thus complete the process of lineage sorting faster than nuclear DNA, resulting in genetic differences between populations, whereas nuclear markers do not reflect yet this differentiation due to their longer coalescence times. The discordance between the two markers used could be also the result of a sex-biased dispersal with males spreading on a wider range than females, leading to nuclear gene flow between populations without a corresponding mitochondrial gene flow. A lower dispersal potential of females has been observed in other mosquitoes, such as in Cx. tarsalis , Ochlerotatus sp.  and Anopheles sp.  hypothesizing that mosquitoes may have a strong home range memory and tend to return to their natal breeding sites for ovipositing and blood-feeding [75, 76]. Another factor that is rarely considered but could be potentially relevant, is selection. Environmental characteristics, such as spatial variation in habitat quality, could favor one mitochondrial variant over another. These variants could be associated with different biological or behavioral features of individuals [77, 78]. In our case, we may assume that mosquito females dispersing outside their own habitat may show a lower fitness than local females, due to different behaviors and ecological requirements (e.g. host choice, mating behavior and oviposition site choice).
Understanding landscape effects on genetic structure could provide insight into biological processes such as metapopulation dynamics and species distribution . In our study, a significant effect of landscape composition on genetic diversity was found, with positive correlation between mitochondrial variability and percentage of wetland areas. An even more significant effect was found considering the distance from the main wetland areas to sampling sites, with the mitochondrial genetic variability decreasing from the main wetland areas along the coast to the inner areas (Fig. 4a-c). Large patches of natural wetlands could maintain higher genetic diversity serving as a source of variability for smaller natural patches, suggesting a model of population expansion between the mainland-island and the source-sink model . Natural environments could be more suitable to harbour populations with higher genetic variability due to the presence of more undisturbed and stable breeding sites sustaining permanent mosquito populations. In addition, in this natural and scarcely urbanized context, where the insecticide treatments are low or absent , mosquito populations may undergo to lower selective pressure thus being less prone to bottlenecks. On the other hand, the low genetic differentiation found in inner and more anthropized sites could also be explained by a chemical control of mosquito populations in these environments. A low genetic differentiation induced by insecticide-resistance has been suggested within the genus Culex [16, 55]. In anthropized areas, the mosquitoes may have experienced population bottlenecks caused by programs of vector-borne disease control (i.e. larviciding, fogging and elimination of breeding sites) .
The Delta Po River is on the route of migratory birds that play an important role in the introduction of viruses, such as WNV, into Europe and the Mediterranean Basin . Here, a high genetic variability of the Cx. pipiens complex, one of the primary WNV vectors, was found using both mitochondrial and nuclear markers. In particular, the mitochondrial diversity was positively correlated to the natural wetland areas close to the Delta Po River. The presence of different haplotypes in the Cx. pipiens complex populations of this area could underlie differences in their biology as it has been already observed for the forms pipiens and molestus [18, 19] that show different feeding behaviors. In this scenario, the cryptic diversity found in the Po River Delta is more interesting as we cannot completely exclude a linkage between the three mitochondrial lineages found and the biology of the mosquitoes and hence their vector ability of arboviral infections, complicating in this way the arbovirus transmission network. Since the link between WNV and deltaic areas seems clear , we think that it should be important to further extend genetic studies of the Cx. pipiens complex in delta areas and other natural habitats. In addition, these studies could be coupled with vector competence studies, because of the potential risk of the different forms as vectors of arbovirus.
Basic Local Alignment Search Tool
Barcode of Life Data Systems
cytochrome oxidase subunit 1
cytochrome oxidase subunit 2
- F st :
Glycine to Serine in position 8
- GTR + I + G:
General Time Reversible + Invariant + Gamma
- H :
- P :
polymerase chain reaction
- R2 :
coefficient of determination
standard error of the mean
sums of squared deviations
West Nile virus
- π :
Vinogradova EB. Culex pipiens pipiens Mosquitoes: Taxonomy, Distribution, Ecology, Physiology, Genetics, Applied importance and Control. Sofia: Pensoft; 2000.
Harbach RE. Culex pipiens: species versus species complex – taxonomic history and perspective. J Am Mosq Control Assoc. 2012;8:10–23.
Vinogradova EB. Ecophysiological and morphological variations in mosquitoes of the Culex pipiens complex (Diptera: Culicidae). Acta Soc Zool Bohem. 2003;67:41–50.
Kent RJ, Harrington LC, Norris DE. Genetic differences between Culex pipiens f. molestus and Culex pipiens pipiens (Diptera: Culicidae) in New York. J Med Entomol. 2007;44:50–9.
Fonseca DM, Keyghobadi N, Malcolm CA, Mehmet C, Schaffner F, Mogi M, et al. Emerging vectors in the Culex pipiens complex. Science. 2004;303:1535–8.
Calistri P, Giovannini A, Hubalek Z, Ionescu A, Monaco F, Savini G, et al. Epidemiology of West Nile in Europe and in the Mediterranean Basin. Open Virol J. 2010;4:29–37.
Hubálek Z, Halouzka J. West Nile fever-a reemerging mosquito-borne viral disease in Europe. Emerg Infect Dis. 1999;5:643–50.
Hamer GL, Kitron UD, Brawn JD, Loss SR, Ruiz MO, Goldberg TL, et al. Culex pipiens (Diptera: Culicidae): a bridge vector of West Nile virus to humans. J Med Entomol. 2008;45:125–8.
Autorino GL, Battisti A, Deubel V, Ferrari G, Forletta R, Giovannini A, et al. West Nile virus epidemic in horses, Tuscany region, Italy. Emerg Infec Dis. 2002;8:1372–8.
Gobbi F, Barzon L, Capelli G, Angheben A, Pacenti M, Napoletano G, Veneto Summer Fever Study Group. Surveillance for West Nile, dengue, and Chikungunya virus infections, Veneto Region, Italy, 2010. Emerg Infect Dis. 2012;18:671–3.
Mulatti P, Bonfanti L, Capelli G, Capello K, Lorenzetto M, Terregino C, et al. West Nile Virus in North-Eastern Italy, 2011: Entomological and Equine IgM-Based Surveillance to Detect Active Virus Circulation. Zoonoses Public Hlth. 2013;60:375–82.
Barzon L, Squarzon L, Cattai M, Franchin E, Pagni S, Cusinato R, et al. West Nile virus infection in Veneto Region, Italy, 2008-2009. Euro Surveill. 2009;14:31.
Barzon L, Pacenti M, Franchin E, Lavezzo E, Masi G, Squarzon L, et al. Whole genome sequencing and phylogenetic analysis of West Nile virus lineage 1 and lineage 2 from human cases of infection, Italy, August 2013. Euro Surveill. 2013;18:38.
Savini G, Capelli G, Monaco F, Polci A, Russo F, Di Gennaro A, et al. Evidence of West Nile virus lineage 2 circulation in Northern Italy. Vet Microbiol. 2012;158:267–73.
Becker N, Jöst A, Weitzel T. The Culex pipiens Complex in Europe. J Am Mosq Control Assoc. 2012;28 Suppl 4:53–67.
Werblow A, Klimpel S, Bolius S, Dorresteijn AW, Sauer J, Melaun C. Population structure and distribution patterns of the sibling mosquito species Culex pipiens and Culex torrentium (Diptera: Culicidae) reveal different evolutionary paths. Plos ONE. 2014;9:e102158.
Shaikevich EV, Zakharov IA. Polymorphism of mitochondrial COI and nuclear ribosomal ITS2 in the Culex pipiens complex and in Culex torrentium (Diptera: Culicidae). Comp Cytogenet. 2010;4:161–74.
Gomes B, Kioulos E, Papa A, Almeida APG, Vontas J, Pinto J. Distribution and hybridization of Culex pipiens forms in Greece during the West Nile virus outbreak of 2010. Infect Genet Evol. 2013;16:218–25.
Gomes B, Sousa CA, Novo MT, Freitas FB, Alves R, Côrte-Real AR, et al. Asymmetric introgression between sympatric molestus and pipiens forms of Culex pipiens (Diptera: Culicidae) in the Comporta region, Portugal. BMC Evol Biol. 2009;9:262.
Pedro P, Sallum MAM. Spatial expansion and population structure of the Neotropical malaria vector, Anopheles darlingi. Biol J Linn Soc Lond. 2009;97:854–66.
Hemmerter S, Šlapeta J, Beebe NW. Resolving genetic diversity in Australasian Culex mosquitoes: incongruence between the mitochondrial cytochrome c oxidase I and nuclear acetylcholine esterase 2. Mol Phylogenet Evol. 2009;50:317–25.
Rasgon JL, Scott TW. Wolbachia and cytoplasmic incompatibility in the California Culex pipiens mosquito species complex: parameter estimates and infection dynamics in natural populations. Genetics. 2003;165:2029–38.
Smith JL, Fonseca DM. Rapid Assays for Identification of members of the Culex (Culex) pipiens complex, their hybrids, and other sibling species (Diptera Culicidae). Am J Trop Med Hyg. 2004;70:339–45.
Kasai S, Komagata O, Tomita T, Sawabe K, Tsuda Y, Kurahashi H, et al. PCR-based identification of Culex pipiens complex collected in Japan. Jap J Infect Dis. 2008;61:184–91.
Bisanzio D, Giacobini M, Bertolotti L, Mosca A, Balbo L, Kitron U, et al. Spatio-temporal patterns of distribution of West Nile virus vectors in eastern Piedmont Region, Italy. Parasites Vectors. 2011;4:230.
Pecoraro HL, Day HL, Reineke R, Stevens N, Withey JC, Marzluff JM, et al. Climatic and landscape correlates for potential West Nile virus mosquito vectors in the Seattle region. J Vector Ecol. 2007;32:22–8.
Calistri P, Giovannini A, Savini G, Monaco F, Bonfanti L, Ceolin C, et al. West Nile Virus Transmission in 2008 in North-Eastern Italy. Zoonoses Public Hlth. 2010;57:211–9.
Atkinson PW, Clark JA, Delany S, Diagana CH, du Feu C, Fiedler W, et al. Urgent preliminary assessment of ornithological data relevant to the spread of Avian Influenza in Europe. In: Delany S, Veen J, Clark J, editors. Report to the European Commission. Wageningen, The Netherlands: Wetlands International; 2006.
Rezza G. Chikungunya and West Nile virus outbreaks: what is happening in north-eastern Italy? Eur J Public Health. 2009;19:236–7.
Lelli R, Calistri P, Bruno R, Monaco F, Savini G, Di Sabatino D, et al. West Nile Transmission in Resident Birds in Italy. Transbound Emerg Dis. 2012;59:421–8.
Weitzel T, Collado A, Jöst A, Pietsch K, Storch V, Becker N. Genetic differentiation of populations within the Culex pipiens complex and phylogeny of related species. J Am Mosq Control Assoc. 2009;25:6–17.
Barr AR. The distribution of Culex p. pipiens and C. p. quinquefasciatus in North America. Am J Trop Med Hyg. 1957;6:153–65.
Severini F, Toma L, Di Luca M, Romi R. Le zanzare italiane: generalità e identificazione degli adulti (Diptera, Culicidae). Fragmenta Entomologica. 2009;41:213–372.
Patwary MU, Kenchington EL, Bird CJ, Zouros E. The use of random amplified polymorphic DNA markers in genetic studies of the sea scallop Placopecten magellanicus (Gmelin, 1791). J Shellfish Res. 1994;13:547–53.
Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3:294–9.
Simon C, Buckley TR, Frati F, Stewart JB, Beckenbach AT. Incorporating molecular evolution into phylogenetic analysis, and a new compilation of conserved polymerase chain reaction primers for animal mitochondrial DNA. Annu Rev Ecol Evol Syst. 2006;37:545–79.
Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: Molecular Evolutionary Genetics Analysis Version 6.0. Mol Biol Evol. 2013;30:2725–9.
Swofford DL. PAUP*: Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4. Sunderland, MA: Sinauer Associates; 2002.
Excoffier L, Lischer HEL. Arlequin suite ver 3.5: A new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10:564–7.
Guindon S, Gascuel O. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003;52:696–704.
Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9:1657–9.
Tajima F. The effect of change in population size on DNA polymorphism. Genetics. 1989;123:597–601.
Rogers AR, Harpending H. Population growth makes waves in the distribution of pairwise genetic differences. Mol Biol Evol. 1992;9:552–69.
Slatkin M, Hudson RR. Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations. Genetics. 1991;129(2):555–62.
Stephens M, Smith NJ, Donnelly P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001;68:978–89.
Stephens M, Scheet P. Accounting for decay of linkage disequilibrium in haplotype inference and missing-data imputation. Am J Hum Genet. 2005;76:449–62.
Martin DP, Murrell B, Golden M, Khoosal A, Muhire B. RDP4: Detection and analysis of recombination patterns in virus genomes. Virus Evol. 2015;1(1):vev003.
Ronquist F, Huelsenbeck JP. MRBAYES 3: bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19:1572–4.
R Core Team. R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria; 2013.
Shaikevich EV. PCR-RFLP of the COI gene reliably differentiates C. pipiens, C. pipiens f. molestus and C. torrentium of the Pipiens Complex. European Mosquito Bull. 2007;23:25–30.
Danabalan R, Ponsonby DJ, Linton YM. A critical assessment of available molecular identification tools for determining the status of Culex pipiens s.l. in the United Kingdom. J Am Mosq Control Assoc. 2012;28:68–74.
Azari-Hamidian S, Linton YM, Abai MR, Ladonni H, Oshaghi MA, Hanafi-Bojd AA, et al. Mosquito (Diptera: Culicidae) fauna of the Iranian islands in the Persian Gulf. J Nat Hist. 2010;44:913–25.
Cook S, Moureau G, Harbach RE, Mukwaya L, Goodger K, Ssenfuka F. Isolation of a novel species of flavivirus and a new strain of Culex flavivirus (Flaviviridae) from a natural mosquito population in Uganda. J Gen Virol. 2009;90:2669–78.
Behura SK, Lobo NF, Haas B, Lovin DD, Shumway MF, Puiu D, et al. Complete sequences of mitochondria genomes of Aedes aegypti and Culex quinquefasciatus and comparative analysis of mitochondrial DNA fragments inserted in the nuclear genomes. Insect Biochem Mol Biol. 2011;41(10):770–7.
Low VL, Lim PE, Chen CD, Lim YA, Tan TK, Norma-Rashid Y. Mitochondrial DNA analyses reveal low genetic diversity in Culex quinquefasciatus from residential areas in Malaysia. Med Vet Entomol. 2013;28:157–68.
Werblow A, Bolius S, Dorresteijn AW, Melaun C, Klimpel S. Diversity of Culex torrentium Martini, 1925 - a potential vector of arboviruses and filaria in Europe. Parasitol Res. 2013;112:2495–501.
Fedorova MV, Shaikevich EV. Morphological and molecular-genetic distinctions between adult mosquitoes Culex torrentium Martini and C. pipiens Linnaeus (Diptera, Culicidae) from Moscow Province. Entomol Rev. 2007;87:127–35.
Guillemaud T, Pasteur N, Rousset F. Contrasting levels of variability between cytoplasmic genomes and incompatibility types in the mosquito Culex pipiens. Proc R Soc Lond. 1997;264:245–51.
Kirkman LK, Whitehead EA, Golladay SW, Smith LL, Opsahl SP. A research framework for identifying potential linkages between isolated wetlands and disease ecology. Ecol Res. 2011;26:875–83.
Johnson JA, Dunn PO, Bouzat JL. Effects of recent population bottlenecks on reconstructing the demographic history of prairie‐chickens. Mol Ecol. 2007;16:2203–22.
Vinogradova EB, Shaikevich EV. Morphometric, physiological and molecular characteristics of underground populations of the urban mosquito Culex pipiens Linnaeus f. molestus Forskål (Diptera: Culicidae) from several areas of Russia. Eur Mosq Bull. 2007;22:17–24.
Morais SA, Almeida FD, Suesdek L, Marrelli MT. Low genetic diversity in Wolbachia-infected Culex quinquefasciatus (Diptera: Culicidae) from Brazil and Argentina. Rev Inst Med Trop Sao Paulo. 2012;54:325–9.
Behbahani A. Wolbachia infection and mitochondrial DNA comparisons among Culex mosquitoes in South West Iran. Pak J Biol Sci. 2012;15:54–7.
Fonseca DM, Smith JL, Kim HC, Mogi M. Population genetics of the mosquito Culex pipiens pallens reveals sex-linked asymmetric introgression by Culex quinquefasciatus. Infect Genet Evol. 2009;9:1197–203.
Lee Y, Seifert SN, Nieman CC, McAbee RD, Goodell P, Fryxell RT, et al. High degree of single nucleotide polymorphisms in California Culex pipiens (Diptera: Culicidae) sensu lato. J Med Entomol. 2012;49:299–306.
Brower AV, Desalle R. Practical and theoretical considerations for choice of a DNA sequence region in insect molecular systematics, with a short review of published studies using nuclear gene regions. Ann Entomol Soc Am. 1994;87:702–16.
Friedlander TP, Regier JC, Mitter C. Nuclear gene sequences for higher level phylogenetic analysis: 14 promising candidates. Syst Biol. 1992;41:483–90.
Friedlander TP, Regier JC, Mitter C. Phylogenetic information content of five nuclear gene sequences in animals: initial assessment of character sets from concordance and divergence studies. Syst Biol. 1994;43:511–25.
Toews DP, Brelsford A. The biogeography of mitochondrial and nuclear discordance in animals. Mol Ecol. 2012;21:3907–30.
Hudson RR, Turelli M. Stochasticity overrules the “three‐times rule”: genetic drift, genetic draft, and coalescence times for nuclear loci versus mitochondrial DNA. Evolution. 2003;57:182–90.
Zink RM, Barrowclough GF. Mitochondrial DNA under siege in avian phylogeography. Mol Ecol. 2008;17:2107–21.
Reisen WK, Lothrop HD, Lothrop B. Factors influencing the outcome of mark-release-recapture studies with Culex tarsalis (Diptera: Culicidae). J Med Entomol. 2003;40:820–9.
Renshaw M, Service MW, Birley MH. Host finding, feeding patterns and evidence for a memorized home range of the mosquito Aedes cantans. Med Vet Entomol. 1994;8:187–93.
McCall PJ, Eaton G. Olfactory memory in the mosquito Culex quinquefasciatus. Med Vet Entomol. 2001;15:197–203.
Charlwood JD, Graves PM, Marshall TDC. Evidence for a ‘memorized’ home range in Anopheles farauti females from Papua New Guinea. Med Vet Entomol. 1988;2:101–8.
McCall PJ, Kelly DW. Learning and memory in disease vectors. Trends Parasitol. 2002;18:429–33.
Irwin DE. Phylogeographic breaks without geographic barriers to gene flow. Evolution. 2002;56:2383–94.
Irwin DE. Local adaptation along smooth ecological gradients causes phylogeographic breaks and phenotypic clustering. Amer Nat. 2012;180:35–49.
Storfer A, Murphy MA, Evans JS, Goldberg CS, Robinson S, Spear SF, et al. Putting the ‘landscape’ in landscape genetics. Heredity. 2007;98:128–42.
Harrison S. Local extinction in a metapopulation context: an empirical evaluation. Metapopulation dynamics: empirical and theoretical investigations. 1991. p. 73–88.
Toma L, Menegon M, Romi R, De Matthaeis E, Montanari M, Severini C. Status of insecticide resistance in Culex pipiens field populations from north-eastern areas of Italy before the withdrawal of OP compounds. Pest Manag Sci. 2011;67:100–6.
Hillis DM, Bull JJ. An empirical test of bootstrapping as a method for assessing confidence in phylogenetic analysis. Syst Biol. 1993;42:182–92.
Huelsenbeck JP, Rannala B. Frequentist properties of bayesian posterior probabilities of phylogenetic trees under simple and complex substitution models. Syst Biol. 2004;53:904–13.
The authors thank the Veneto Regional Po Delta Park for facilities and logistic support, during the field work. We thank Gioia Capelli and Edoardo Petrucco Toffolo for constructive and helpful comments. Special thanks to the Prof. Vincenzo Girolami for critical revision of the manuscript and helpful suggestions. This study was partially supported by the ex-60 % projects of the University of Padua (60A08-7389/13) and by the regional project “Pilot project for the control of mosquitoes in the Rovigo Province” funded by Veneto Region and by Parco Delta del Po. The Department of Agronomy, Food, Natural resources, Animals and Environment (DAFNAE) of University of Padova provides all other support for this study.
The authors declare that they have no competing interests.
AS and LM planned the research, MS IMS AS and LM performed the research, MS IMS GC and GS analyzed the data, MS IMS and LM wrote the paper and revised the manuscript. All authors read and approved the final manuscript.
GenBank accession numbers for the haplotypes found. (DOCX 18 kb)
Phylogenetic tree of Culex pipiens complex based on mitochondrial COI and COII markers. (PDF 286 kb)
Summary of Culex pipiens complex samples analyzed with reference to individual haplotypes. (DOCX 50 kb)
About this article
Cite this article
Simonato, M., Martinez-Sañudo, I., Cavaletto, G. et al. High genetic diversity in the Culex pipiens complex from a West Nile Virus epidemic area in Southern Europe. Parasites Vectors 9, 150 (2016) doi:10.1186/s13071-016-1429-1
- Culex pipiens complex
- Mitochondrial DNA
- Cytochrome oxidase
- Nuclear DNA
- Genetic diversity
- West Nile Virus