Mosquito samples

Culex pipiens taxa discrimination

Culex quinquefasciatus and Cx. pipiens were discriminated using a specific PCR assay based on the acetylcholinesterase-2 gene (ACE2-assay) [18]. Identification of the form molestus and the form pipiens of Cx. pipiens was made on the basis of the ovarian status of females (autogeny) and genetically by both CQ11 [19] and COI [25] assays. The expression of autogeny was studied in the laboratory, for most underground and above ground populations from Russia, insectary colony T7 from France, and in samples from Portugal by the respective collectors (Table 1). Autogeny of Cx. pipiens from the other populations was not studied. Individuals whose autogeny status had not been determined, and that came from collections which included either pipiens or molestus, as well as hybrid forms according to CQ11 analysis, were denominated as Cx. pipiens “hybrid” and those to whom neither autogeny determination nor CQ11 assay had been performed or the results of CQ11 were inconclusive, were denominated as “unknown” in Table 1.

Molecular typing

DNA was extracted from mosquitoes preserved in ethanol using the DIAtom™ DNA Prep Kit (Isogen, Moscow, Russia). Polymerase chain reactions were run in thermocyclers GeneAmpR PCR System 2700 (Applied Biosystems, Foster City, CA, USA) with amplification kits GenePak™ PCR Core (Isogene, Moscow, Russia).

Mitochondrial DNA typing

These sequences were compared with 24 previously studied full-size DNA sequence of the COI gene (1548 bp) of both forms of Cx.pipiens originated from 10 geographically distinct sample sites from Russia (Gene Bank accession numbers FN395171-FN395190) and Cx. quinquefasciatus originated from two sample sites from India (FN395201- FN395204) [40].

Culex mtDNA of 580 specimens from 35 geographical populations was also genotyped using a series of specific PCR-RFLP (restriction fragment length polymorphism) assays based on the DNA variability of COIgene [25, 43]. The 5′ region of the COI gene of 603 bp was amplified using primers CulexCOIF and CulexCOIR [25]. PCR conditions were the following: primary denaturing – 5 min at 94°С; 35 cycles: denaturing at 94°С – 30 s, annealing at 55°С – 30 s, synthesis at 72°С – 40 s; final synthesis at 72°С – 10 min. HaeIII digestion of the COI PCR products allowed the discrimination of D haplotype from A, B, C and E haplotypes [25, 43]. However, HaeIII has no recognition site for GG’CC on the COI sequence of type D and PCR-product of amplification remains unchanged, 603 bp, whereas the COI fragment of other types resulted in two fragments - 397 and 206 bp respectively (Additional file 1A) - after restriction with HaeIII. After digestion with AluI, the COI PCR-products of A, B, C and D types resulted in 8 fragments (189, 171, 99, 67, 45, 15, 12, 5 bp), 5 of which are visible by electrophoresis in a 2 % agarose gel. The COI fragment of type E and E1 is cut into 7 fragments (189, 171, 144, 67, 15, 12, 5 bp) only because the mutation at position 206 (Fig. 1) inactivates the AluI restriction site (AG’CT). Consequently, the 144 bp fragment is diagnostic for haplotypes E and E1 allowing their discrimination (Additional file 1B). Reactions were carried out in a restriction mixture consisting of 5 μl COI PCR product, 0.2 μl (2 units) of the enzyme, 3 μl buffer, 0.3 μl BSA and 21.5 μl ddH2O. Both HaeIII and AluI restriction mixtures were incubated for 2 hours at 37 °C. At least two technical replicates were performed. Results were visualized by electrophoresis in a 2 % agarose gel. As it is not known whether or not the insectary lines descended from one or more female founders, mtDNA polymorphisms were studied for 10–20 individuals per population.

Nuclear DNA typing

Complex species and form identification, using ACE2 (with primers B1246s, ACEpip and ACEquin) and CQ11 assays, respectively were performed as described by authors [18, 19]. After amplification with primers B1246s, ACEpip and ACEquin, PCR product of 610 bp is characteristic for Cx. pipiens and of 274 bp for Cx. quinquefasciatus [18]. The CQ11 PCR-product approximately 200 bp in size is characteristic for Cx. pipiens form pipiens and 250 bp for form molestus. However, CQ11 amplicons from Cx. quinquefasciatus yields a PCR product of 250 bp too; therefore the authors recommend the use of a combination of both tests, CQ11 and ACE2 in areas of sympatry of the two species [19].

ACE2 sequencing

For amplification of ACE2 gene of samples from Kos we used the primers F1457 5′–GAGGAGATGTGGAATCCCAA–3′ and B1246 5′–TGGAGCCTCCTCTTCACGGC–3′ and PCR conditions described earlier [16]. Amplicons of the ACE2 gene were approximately 710 bp. PCR products were excised from a 1 % agarose gel and purified using a QIAquick Gel Extraction kit (Qiagen, USA) according to the manufacturer’s instructions. The PCR products from two mosquito were cloned using the kit pGEM-T Easy Vector Systems (Promega, USA). The clones were screened for the presence of different ACE2 alleles by PCR-RFLP test: the restriction enzyme Sau3AI (Fermentas) cuts the “pipiens” allele into three fragments (330 bp, 213 bp and 167 bp), and the “quinquefasciatus” allele into two fragments (543 bp and 167 bp). Bacterial cells suspension after denaturation in boiling water bath used as a DNA template in PCR with F1457 and B1246 primers and 10 μL of the PCR product after amplification were digested with 2 units of enzyme for 2 hours at 37 °C. PCR products from individual bacterial clones were sequenced from both fragment’s ends using the equipment ABI PRISM 310 and the BigDye Termination kit (Applied Biosystems, USA), according to the manufacturer’s instructions. Sequences were analyzed using Chromas software (http://www.technelysium.com.au) and two different alleles of ACE2 gene were deposited in GenBank under numbers KU163609-10.

Wolbachia infection typing

The wPip infections were genotyped de novo in a subsample of 274 individuals representative of all COI haplotypes and assigned to one group (wPip-I to wPip-V) using PCR-RFLP assays based on two Wolbachia pipientis markers, ank2 and pk1, as previously described [30]. PCR-products about 310 and 510 bp were obtained after ank2 amplification. Specific pk1 PCR-products were approximately 1350 bp in size. The digestion with endonuclease HinfI of the ank2 PCR products provided three alleles: а (313 bp), b (217, 195, 98 bp) and с (293, 217 bp). After digestion of the pk1 PCR products with endonuclease TaqI, three alleles were obtained: a/e (903, 430 bp), c (851, 498 bp) and d (497, 251, 107 bp). The alleles a and e of pk1 gene have the same fragment sizes and therefore needed additional treatment of the pk1 PCR products by restriction endonuclease PstI, resulting in two alleles: a (903, 303, 141 bp) and e (903, 430 bp). After this additional digestion with PstI, alleles a, c and d of the pk1 gene were obtained. According to Atyame and coauthors (2011) [32], different alleles of pk1 and ank2 genes correspond to one of five groups, wPip-I to wPip-V.

Data analysis

Forty-eight COI sequences, 1150 bp, the origin of which is shown in Table 2, were analyzed using the software Chromas (http://www.technelysium.com.au) and aligned and analyzed using MEGA v. 6.0 [44].

Mitochondrial COI haplotype network analysis was performed for 48 sequences using statistical parsimony with the program TCS v.1.21 [45]. The network connection limit was set at 95 %. The resulting networks identify both the relationship between the different haplotypes as well as the number of substitution connecting haplotypes.

The evolutionary history was inferred by using the Maximum Likelihood method based on the Hasegawa-Kishino-Yano model [46]. The tree with the highest log likelihood (−1539.7457) is shown. Initial tree(s) for the heuristic search were obtained by applying the Neighbor-Joining method to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach. The tree was drawn to scale, with branch lengths measured in the number of substitutions per site. Bootstrap coefficients were calculated for a 1000 repeats. All positions containing gaps and missing data were eliminated.

Chi-square for association between COI haplotype and mosquito taxa was tested with GraphPad InStat (www.graphpad.com accessed in May 2014) based on the RLFP analysis of 544 samples, specimens for which the taxonomic status was not determined, i.e. 36 specimens denoted as “unknown” (Table 1), were excluded from the analysis.

Polymorphism of the DNA sequence of the mitochondrial gene COI

Based on the differences in the nucleotide composition of the gene COI, we found 6 mitochondrial haplotypes in mosquitoes of the Culex pipiens complex, denoted as A, B, C, D, E, and a new one found in this work, E1. These haplotypes vary in seven substitutions (Fig. 1). Haplotype D is characterized by two fixed substitutions in positions 119 and 896 when compared with A, B, C or E and E1 types. By an additional single mutation haplotype B (in position 467) and haplotype C (in position 677) differ from haplotype A. Haplotypes E and E1 differ from A, B, C and D in positions 206 and 848, while haplotype E1 differs from E, by a single additional mutation in position 830.

Six sequences from Portuguese mosquitoes were only sequenced in the second half of the gene (one sequence with 486 bp and the other five with ca. 705–784 bp), which however, included the site where haplotypes E and E1 differ. These were 4 molestus with haplotype E1, 1 pipiens and 1 molestus with haplotype E. So, we temptatively included them into analysis, for a total of 54 sequences (Table 2). Among the 54 samples investigated, 6 specimens had haplotype A, 2 - haplotype B, 6 - haplotype C, 22 - haplotype D, 12 - haplotype E and 6 - haplotype E1 (Table 2). Haplotypes A and C were found in samples from Russia and Germany. Haplotype B was found in one population only: Luzki, from Russia. Haplotype D was found in Russia, Germany, northern Italy (Piedmont), Tunisia and Morocco. Haplotype E was found in India, Italy (Viterbo), Israel and Greece (Cyprus). In Portugal, two haplotypes were found: E and E1.

PCR-RFLP assays of the COI gene

The differences in the nucleotide composition of the 5′ region of the COI gene made it possible to choose restriction endonucleases for PCR-RFLP assay [25, 43]. Characteristic 603 bp amplification products were obtained for all Culex spp. mosquitoes being studied. A first assay using HaeIII restriction endonuclease made it possible to identify COI type D. Two hundred and 32 mosquitoes with haplotype D were all originated from underground sampling sites from Russia, indoor sites from Germany and also open habitats from Italy, Tunisia and Morocco.

The second assay using AluI restriction endonuclease made it possible to identify E and E1 haplotypes. Sequences of E and E1 types cannot be distinguished using the PCR-RFLP method. Haplotypes E (E1) were found in 125 of the specimens being examined: mosquitoes from Italy, Portugal, Greece, Israel, Morocco and India.

Nuclear locus analysis

Since in some cases inconsistencies between the taxonomic status of Cx. pipiens and expected type of COI (as with French T7 laboratory line or with Portuguese mosquitoes) have been detected (Table 1), we tested nuclear DNA polymorphism of such markers as ACE2 [18] and microsatellite marker CQ11 [19] to clarify the taxonomic status and to reveal possible cases of hybridization. ACE2 and CQ11 assays could not be performed for all studied individuals, namely from Israel and Viterbo, Italy, due to limitation of DNA availability and budget constraints. In some populations, taxa were defined according to known and previously verified autogeny (namely Russian and Portuguese) and known lab line origin (Indian), of which CQ11 was determined in a subsample of 274 specimens (Table 3).

The ACE2 marker allows for differentiation between Cx. pipiens (without separation into forms) and Cx. quinquefasciatus. After the amplification with the primers B1246s, ACEpip and ACEquin the majority of mosquito samples yielded the PCR product of 610 bp, characteristic for Cx. pipiens (Table 3). The exceptions were specimens collected on the Greek island of Kos (13 of the 24 samples) [13], and specimens from Tanger, Morocco (7 of the 13 samples) in which we found specific PCR products for Cx. pipiens (610 bp) and for Cx. quinquefasciatus (274 bp), i.e. these samples corresponded to hybrids between these taxa (Additional file 2).

The test based on polymorphism on the flanking region of the microsatellite locus CQ11 has been designed to identify both forms of Cx. pipiens (form pipiens and form molestus) and their hybrids. However, the same size PCR-product is obtained for Cx. quinquefasciatus and Cx. pipiens form molestus. Therefore, we took into account the results of both methods, ACE2 and CQ11 as is recommended by authors [19]. In Tanger, Morocco, 6 samples were Cx. pipiens/Cx. quinquefasciatus hybrid and 2 samples were Cx. pipiens form pipiens, according to both assays, 4 samples were Cx. pipiens by ACE2 and hybrid by CQ11 and 1 sample was Cx. pipiens/Cx. quinquefasciatus hybrid by ACE2 and Cx. pipiens form pipiens by CQ11 PCR-results (Additional file 2). The discrepancy in the results of the analysis based on ACE and CQ11 loci were also obtained in samples from Kos, Greece [13]. Such cases indicate recombination processes in the hybrid population and CQ11 hybrids from these collections are likely to be regarded as hybrids between Cx. pipiens and Cx. quinquefasciatus rather than hybrids between pipiens and molestus (11 samples in Tanger and 22 samples in Kos, total 33 in Table 4).

In order to ascertain these PCR results, a larger part of ACE2 gene was cloned and sequenced. Particularly the region in which ACE marker is included, namely part of exon 2, intron 2 and part of exon 3, as described in [16]. Analysis of nucleotide sequences confirmed the occurrence of two alleles of the ACE2 gene in one individual mosquito (alignment shown in Additional file 3). Following Blast analysis, one allele of ACE2 gene of sample Kos1 in our study was completely similar to Cx. pipiens ACE2 gene sequences (Accession No. AY196910), while the second allele of ACE2 gene of the same sample Kos1 was similar to sequence of this gene from Cx. quinquefasciatus, 99 % identity with all published sequences in GenBank (for example, Accession No. AY196911), confirming PCR results suggestive of “pipiens/quinquefasciatus” hybrid.

According to the CQ11 analysis of 274 specimens, hybrids between pipiens and molestus were detected in almost all samples collected in the Mediterranean region, irrespective of the type of cytoplasmic structures (Table 3). Hybrids between pipiens and molestus determined by this method were also found earlier in populations from Portugal [20] and Morocco [23], from which some specimens examined in this work were taken.

Typing of Wolbachia polymorphism

The association between COI and symbiotic intracellular bacterium Wolbachia pipientis was studied in 274 Cx. pipiens sl. individuals representative of all COI haplotypes (Table 3). Infection with Wolbachia was detected in all specimens examined and specific ank2 and pk1 PCR-products were observed. In our collection we found the three known ank2 alleles (a, b and c) and the three known pk1 (a, c and d) alleles (Additional file 4). Using specific PCR-RFLP assays enabled us to genotype and assign the wPip infections of each specimen to a group, from wPip-I to wPip-V. All individuals with mtDNA haplotypes A, B or C appeared to be infected with wPip-II whereas those with mtDNA haplotype D with wPipIV, and those with mtDNA haplotype E and E1 with wPipI (Table 3).

Based on the fact that the haplotypes A, B and C are close, suggesting that haplotypes B and C are derived from haplotype A (Fig. 1) and also on the fact that these haplotypes are transmitted in Cx. pipiens together with the bacterium wPip-II, we have combined them into a group of mitochondrial haplotypes here denominated - haplogroup A. Similarly, haplotype E1 is likely derived from haplotype E and both are transmitted in association with the wPipI infection, so we named them haplogroup E. The only other haplotype we detected was haplotype D, which is linked with wPipIV infection. Thus, we have identified three groups of mitochondrial haplotypes of COI gene that are associated, respectively, with three groups of symbiotic bacteria in Culex pipiens complex mosquitoes, similar with polymorphism of W. pipientis and groups of mitochondrial haplotypes of cytB [35].

A highly statistically significant correlation was observed between COI type and Taxa (Table 4). Haplogroup A (A, B, C) and wPipII was found in 201Cx. pipiens form pipiens specimens, 11Cx. pipiens form molestus and 10 hybrid based on CQ11 assay from southern France and Morocco, and 1 sample from France, (St-Nazaire de Pezan) with negative CQ11 results, the so called 0 allel, denoted as “unknown” in Table 1. Haplotype D and wPipIV was found in 179 Cx.pipiens form molestus from northern European counties and in 20 Cx.pipiens form pipiens specimens, 21 hybrid and in 12 specimens with “unknown” status from Italy, Piedmont and Tunisia. Haplogroup E (E, E1) and wPipI was found in 43 Cx.quinquefasciatus samples, 33 Cx. pipiens/Cx. quinquefasciatus hybrids, 9 Cx. pipiens form pipiens and 14 Cx.pipiens form molestus, 3 pipiens /molestus hybrids and 23 Cx. pipiens mosquitoes with “unknown status” (1 from Greece, Cyprus with negative result after CQ11 assay, and 22 Cx. pipiens mosquitoes from Italy and Israel whose status was not determined by CQ11 assay) (Table 4). Within the 544 specimens with known taxonomic status (580 minus 36 “unknown” status (Table 1), a strong association was observed between COI haplotype or group and taxa, being group A more frequent in Cx. pipiens form pipiens (91 %), type D in Cx. pipiens form molestus (81 %), and group E in Cx.quinquefasciatus (42 %) and in its hybrids with Cx. pipiens (32 %) (Chi-square = 732.71, d.f. = 8, P < 0.0001).

Phylogenetic analysis

Molecular Phylogenetic analysis of 48 sequences with 1150 bp of COI gene from Cx. pipiens taxa, was carried out using the Maximum Likelihood method based on the Hasegawa-Kishino-Yano model (Fig. 2). The average total A + T content was 70.4 % and both variable or Parsimony informative sites were 0.006 % (7/1150). The tree topology supports the data from our RFLP analysis. However, it does not confirm a phylogenetic relation between taxa groups, as the bootstrap values uniting the different groups are quite low (<75), with the exception of haplotype D and haplotypes E and E1 with bootstrap values of 87.