Introduction
Accurate taxon identification is a crucial element of vector surveillance and control. Namely, different vector species exhibit varying potential to transmit pathogens, a property which depends on the vector's ecology and physiology (Talbalaghi & Shaikevich, Reference Talbalaghi and Shaikevich2011), and it is therefore important to precisely determine the vector composition of a certain area. This problem is particularly prominent in species complexes as they represent groups of taxa that are morphologically difficult or impossible to distinguish and which often showcase sharp differences in physiological and behavioral properties, the traits which influence their vectorial capacity and epidemiologic significance (Bahnck & Fonseca, Reference Bahnck and Fonseca2006). To overcome this setback, a plethora of molecular approaches for taxon identification was developed, the most popular one being DNA barcoding (Hebert et al., Reference Hebert, Cywinska, Ball and de Waard2003).
One of the vector species complexes whose taxonomy represents an outstanding problem is the Culex (Cx.) pipiens complex. It is a cosmopolitan mosquito assemblage that includes vectors of various viral, protozoan, and metazoan pathogens of people and animals. Some of the associated diseases include: West Nile encephalitis, Japanese encephalitis, avian malaria and filariasis (Shaikevich, Reference Shaikevich2017), and the spectra of affiliated pathogens keeps expanding with the emergence of novel studies (see El-Kholy et al., Reference El-Kholy, El-Husseiny, Meshrif, Abou El-Azm and Salem2017 for hepatitis C virus implications). Despite physiological and behavioral differences, the members of the Cx. pipiens complex exhibit a high degree of morphological similarity, making it difficult to precisely discriminate between them (Becker et al., Reference Becker, Jöst and Weitzel2012). The morphological similarity of the taxa included in the assemblage is particularly pronounced in its nominal species, Cx. pipiens L. Namely, it includes two subspecies, Cx. p. pipiens and Cx. p. pallens Coquillett (www.mosquito-taxonomic-inventory.info), with the former further comprised of two morphologically almost identical, yet physiologically and behaviorally distinct ecotypes: molestus Forskål (autogenous, stenogamous, homodynamic, and mammophilic) and pipiens (anautogenous, eurygamous, heterodynamic, and ornithophilic) (Becker et al., Reference Becker, Jöst and Weitzel2012). Their physiological and behavioral differences are interpreted as adaptations to different habitats: while molestus mosquitoes usually inhabit underground, urban areas, pipiens mosquitoes are usually found in aboveground, rural habitats (Chevillon et al., Reference Chevillon, Eritja, Pasteur and Raymond1995).
Due to their epidemiological distinctiveness (Di Luca et al., Reference Di Luca, Toma, Boccolini, Severini, La Rosa, Minelli, Bongiorno, Montarsi, Arnoldi, Capelli, Rizzoli and Romi2016), it is necessary to delineate the two Cx. p. pipiens ecotypes. However, apart from generally unfruitful attempts using morphological characters (see Kent et al., Reference Kent, Harrington and Norris2007 for summary, but also Vinogradova et al., Reference Vinogradova, Reznik and Kuprijanova1996; Krtinić et al., Reference Krtinić, Ludoški and Milankov2012, Reference Krtinić, Ludoški and Milankov2015), their delimitation is problematic since they were defined using phenotypic differences, which are prone to environmental modulation (Chevillon et al., Reference Chevillon, Eritja, Pasteur and Raymond1995). To overcome these difficulties, the issue of the molestus/pipiens differentiation has been tackled by a multitude of molecular approaches. Some of the markers include: the 3’ end of mitochondrial cytochrome c oxidase subunit I gene (Shaikevich et al., Reference Shaikevich, Vinogradova, Platonov, Karan and Zakharov2005; Kent et al., Reference Kent, Harrington and Norris2007), microsatellite loci (Fonseca et al., Reference Fonseca, Keyghobadi, Malcolm, Mehmet, Schaffner, Mogi, Fleischer and Wilkerson2004; Kent et al., Reference Kent, Harrington and Norris2007; Gomes et al., Reference Gomes, Sousa, Novo, Freitas, Alves, Côrte-Real, Salgueiro, Donnelly, Almeida and Pinto2009, Reference Gomes, Sousa, Vicente, Pinho, Calderón, Arez, Almeida, Donnelly and Pinto2013; Kothera et al., Reference Kothera, Nelms, Reisen and Savage2013), the flanking region of CQ11 microsatellite locus (CQ11FL) (Bahnck & Fonseca, Reference Bahnck and Fonseca2006; Gomes et al., Reference Gomes, Sousa, Novo, Freitas, Alves, Côrte-Real, Salgueiro, Donnelly, Almeida and Pinto2009, Reference Gomes, Sousa, Vicente, Pinho, Calderón, Arez, Almeida, Donnelly and Pinto2013; Di Luca et al., Reference Di Luca, Toma, Boccolini, Severini, La Rosa, Minelli, Bongiorno, Montarsi, Arnoldi, Capelli, Rizzoli and Romi2016), intergenic spacer (Shaikevich et al., Reference Shaikevich, Zagoskin and Mukha2013), 12S rRNA, HS60, ND4 (Kent et al., Reference Kent, Harrington and Norris2007), and most recently, period and timeless, circadian rhythm genes (Shaikevich et al., Reference Shaikevich, Karan and Fydorova2016), as well as a genome-wide AFLP approach (Gomes et al., Reference Gomes, Wilding, Weetman, Sousa, Novo, Savage, Almeida, Pinto and Donnelly2015). The issue of molestus/pipiens identification, however, remains a controversial topic since different DNA regions are subject to varying levels of gene flow as a consequence of their mode of inheritance (Petit & Excoffier, Reference Petit and Excoffier2009), which is reflected in their diagnostic value.
Moreover, a regionally varying migration/selection balance was proposed in Cx. p. pipiens ecotypes (Chevillon et al., Reference Chevillon, Rivet, Raymond, Rousset, Smouse and Pasteur1998). Namely, it has been suggested that the degree of genetic exchange between molestus and pipiens populations changes in a latitude-dependent manner, being the greatest in the Mediterranean region and gradually declining toward northern areas (Chevillon et al., Reference Chevillon, Rivet, Raymond, Rousset, Smouse and Pasteur1998). Such a strengthening of the genetic connectedness of below- and aboveground Cx. p. pipiens populations was accompanied by a decreased selection against maladaptive alleles in aboveground populations. Hence, putatively diagnostic alleles may be variably distributed among geographic molestus and pipiens populations, due to the varying balance between homogenizing gene flow and divergent selection (see Chevillon et al., Reference Chevillon, Rivet, Raymond, Rousset, Smouse and Pasteur1998 for Aat-1 locus example).
The success of Cx. p. pipiens ecotype delimitation therefore depends not only on the choice of a genetic locus, but also on the geographic region where the mosquitoes were sampled due to the different degree of interecotypic hybridization. Therefore, the diagnostic value of molecular assays which were only locally tested should be evaluated on a larger spatial scale. The establishment of regional DNA barcode libraries of mosquito species drove the accumulation of mitochondrial and nuclear sequences of molestus and pipiens ecotypes which were not necessarily discussed in regards to their delimitation. Consequently, it became possible to transregionally review molecular assays proposed for Cx. p. pipiens ecotype designation. Therefore, the first goal of our study was to assess the applicability of assays based on the sequence polymorphism in the 5’ end of mitochondrial cytochrome c oxidase subunit I (5’ COI mtDNA) and internal transcribed spacer 2 (ITS2 rDNA) loci, which were so far locally tested [e.g., in Russia (Shaikevich, Reference Shaikevich2007) and the UK (Danabalan et al., Reference Danabalan, Ponsonby and Linton2012) for 5’ COI mtDNA and the USA (Crabtree et al., Reference Crabtree, Savage and Miller1995) and Russia (Vinogradova et al., Reference Vinogradova, Shaikevich and Ivanitsky2007) for ITS2 rDNA] and suggested as diagnostic and not diagnostic, respectively. Next, we investigated whether the patterns of genetic structuring on a global scale using these two markers would correspond to molestus/pipiens sample grouping. The use of both 5’ COI mtDNA and ITS2 rDNA allowed us to study the pattern of genetic differentiation using markers which differed in inheritance mode (uniparental vs. biparental inheritance) and mutational rate (mitochondrial vs. nuclear gene). Finally, concerning molestus/pipiens designation, we wanted to identify Cx. p. pipiens mosquitoes sampled in Serbia using these two assays in the case of their applicability, and compare the results with CQ11FL identification, suggested as the most reliable locus for this purpose (Di Luca et al., Reference Di Luca, Toma, Boccolini, Severini, La Rosa, Minelli, Bongiorno, Montarsi, Arnoldi, Capelli, Rizzoli and Romi2016).
A separate issue important for vector and disease control is the scale at which management strategies should be implemented. Namely, in the case of gene flow between urban and rural areas, pathogens and insecticide resistance alleles could be exchanged between them. Therefore, it is also necessary to investigate the genetic connectedness of urban and rural mosquito samples (see Krtinić et al., Reference Krtinić, Francuski and Milankov2014). Thus, we also genotyped Cx. p. pipiens individuals originating from aboveground urban and rural samples in the Vojvodina Province (Serbia) at COI mtDNA and ITS2 rDNA loci to test whether gene flow kept them connected.
Materials and methods
The transregional review of 5’ COI mtDNA and ITS2 rDNA assays
To review the diagnostic value of 5’ COI mtDNA and ITS2 rDNA assays for Cx. p. pipiens ecotype identification, we analyzed sequences previously published in GenBank (table 1). Out of the total number of 5’ COI mtDNA and ITS2 rDNA Cx. p. pipiens sequences present in GenBank, we selected a subsample where the mosquitoes were identified up to the ecotype level using one or several ecotype-defining criteria (table 1). Additionally, sequences obtained de novo in this study from mosquitoes sampled on the territory of the Autonomous Province of Vojvodina were included (to determine with which ecotypes they cluster together), making the total of analyzed sequences 159 (26 Serbian and 133 from GenBank) and 47 (21 Serbian and 26 from GenBank) for 5’ COI mtDNA and ITS2 rDNA marker, respectively. In the case of 5' COI mtDNA, the primer combination we used (discussed below) amplified a fragment which is larger than the standard, approximately 650 bp-long barcode fragment, amplified by the LCO-1490/HCO-2198 (Folmer et al., Reference Folmer, Black, Hoeh, Lutz and Vrijenhoek1994) primer pair. Furthermore, as different primer combinations were used to sequence 5’ COI mtDNA of Cx. p. pipiens (see Batovska et al., Reference Batovska, Blacket, Brown and Lynch2016), after the alignment of GenBank accessions with our sequences, an overlapping 5’ COI mtDNA fragment of around 500 bp was used for the review of the assay's diagnostic power. For ITS2 rDNA, we used the entire amplified region of around 325 bp as it overlapped with the GenBank sequences. The number of haplotypes and variable positions was determined in DAMBE6 (Xia, Reference Xia2017), while haplotype and allele networks were constructed in Network 5.0.0.0. (Fluxus Technology Ltd, Clare, Suffolk, the UK) using the median joining approach. The diagnostic power of DNA barcode assays was inferred from the existence of a barcode gap. Namely, when the range of interspecific diversity (in our case, divergence between the ecotypes) is greater than intraspecific diversity (here, distance within the ecotypes), the barcode gap is established and the barcoding approach is deemed successful (Puillandre et al., Reference Puillandre, Lambert, Brouillet and Achaz2011). The existence of the barcode gap between molestus and pipiens sequences was investigated by retrieving p distances between them using MEGA7 software (Kumar et al., Reference Kumar, Stecher and Tamura2016).
Ecotype membership, region of origin, source article, ecotype-defining criteria, GenBank accession codes, sample size, and COI mtDNA haplotype/ITS2 rDNA allele representation of Cx. p. pipiens mosquitoes whose DNA sequences were retrieved from the GenBank and analyzed in the study.
1 N – number of individuals per region.
Criteria: A – autogeny, B – habitat type, C – distributional data, D – siphonal index of larvae, E – COI mtDNA assay, F – CQ11FL assay, G – unknown.
Serbian sample genotyping and the analyses of genetic differentiation
The larvae of Cx. p. pipiens were collected from seven areas of the Autonomous Province of Vojvodina (leg. and det. Krtinić, B.) during the period from 2009 to 2010. The sampling design and map with localities of the studied Serbian samples was provided in detail by Krtinić with coauthors (Table S1, 2016). Data on 28 individuals were used for 5’ COI mtDNA and ITS2 rDNA analyses (table 2). The mosquitoes were divided into two groups, defined by the aboveground habitat type where the larvae were sampled (sensu Krtinić et al., Reference Krtinić, Francuski and Milankov2014): SII sample, draining ditches and ponds in the city (urban habitat), and SIII sample, ponds and swamps outside urban area (rural habitat). The belowground, urban sample (SI) from the previous studies (Krtinić et al., Reference Krtinić, Ludoški and Milankov2012, Reference Krtinić, Francuski and Milankov2014, Reference Krtinić, Ludoški and Milankov2015) was not included since only the sample from Novi Sad exists and we were interested in the assays’ behavior at the regional (Vojvodina Province) scale. Following the metamorphosis of larvae into pupae and adults, the collected mosquitoes were determined as Cx. p. pipiens based on morphological characters (Gutsevich et al., Reference Gutsevich, Monchadskii and Shtakel'berg1974; Božičić, Reference Božičić1985).
Population origin, sample membership, sample size, code, and COI mtDNA haplotype/ITS2 rDNA allele representation of the Serbian Culex p. pipiens individuals genetically analyzed in the study.
The maceration of tissue and the extraction of total DNA was carried out using NucleoSpin® Tissue DNA extraction kit (MACHEREY-NAGEL, Düren, Germany) following the manufacturer's protocol, after which the DNA was amplified by PCR using an illustra PuReTaq Ready-To-Go PCR Beads kit (GE Healthcare Life Sciences, Buckinghamshire, the UK). To amplify 5′ COI mtDNA, we used the forward LCO-1490 (5'-GGTCAACAAATCATAAAGATATTGG-3′, Folmer et al., Reference Folmer, Black, Hoeh, Lutz and Vrijenhoek1994) and reverse Inger (5′-AAAAATGTTGAGGGAAAAATGTTA-3′, UEA8 in Lunt et al., Reference Lunt, Zhang, Szymura and Hewitt1996) primer pair, while the amplification of ITS2 rDNA was done using the forward ITS2A (5′-TGTGAACTGCAGGACACAT-3′) and reverse ITS2B (5′-TATGCTTAAATTCAGGGGGT-3′) primer pair (Beebe & Saul, Reference Beebe and Saul1995), with PCR conditions as described in Milankov et al. (Reference Milankov, Ludoški, Ståhls, Stamenković and Vujić2009). To check the success of reactions, amplification products were separated on a 2% agarose gel. PCR products were then purified using ExoSAP-IT™ PCR Product Cleanup Reagent (Thermo Fisher Scientific, Vilnius, Lithuania), and bidirectionally sequenced on ABI3730XL by Macrogen (The Netherlands).
Chromatograms obtained by 5′ COI mtDNA and ITS2 rDNA sequencing were checked in Chromas 2.6 (Tehnelysium Pty Ltd, South Brisbane, Australia) for erroneously called bases, while sequence alignment was performed in BioEdit 7.2.5 (Hall, Reference Hall1999). The retrieved 5′ COI mtDNA and ITS2 rDNA sequences were uploaded to GenBank (MK607060-085 and MK584719-739). The number of haplotypes and variable positions was determined in DAMBE6 (Xia, Reference Xia2017), while haplotype and allele networks were constructed in Network 5.0.0.0. (Fluxus Technology Ltd, Clare, Suffolk, the UK) using the median joining approach. The existence of stop codons within the 5′ COI sequences was investigated by using an online Translate tool (https://web.expasy.org/translate/).
Mosquitoes belonging to the Cx. p. pipiens urban and rural samples from Serbia were also identified to the ecotype level using the CQ11FL assay (Bahnck & Fonseca, Reference Bahnck and Fonseca2006). The ecotype-specific alleles were amplified using primers CQ11F2 (5′-GATCCTAGCAAGCGAGAAC-3′), pipCQ11 (5′-CATGTTGAGCTTCGGTGAA-3′), and molCQ11 (5′-CCCTCCAGTAAGGTATCAAC-3′), while PCR conditions were as defined by Bahnck & Fonseca (Reference Bahnck and Fonseca2006). The reaction products were separated on a 2% agarose gel and visualized under UV light. The band characteristic for the pipiens ecotype is around 200 bp-long and (around 70 bp) shorter than the one characteristic for molestus, while hybrids are characterized by the presence of both alleles (Bahnck & Fonseca, Reference Bahnck and Fonseca2006).
The analyses of genetic differentiation were performed transregionally (global molestus vs. pipiens ecotype sample) and intraregionally (Serbian urban vs. rural sample). We applied the Bayesian model-based clustering algorithm implemented in BAPS 6.0 (Corander & Tang, Reference Corander and Tang2007), jointly assigning 5′ COI mtDNA and ITS2 rDNA GenBank sequences to artificial populations corresponding to the regions of origin. Clustering among populations was performed using Clustering with linked loci and codon linkage model. We ran BAPS for values of K ranging from one to ten, performing five replicates for each K. The results of the population mixture analysis were used as an input file for population admixture analysis using Admixture based on mixture clustering method. The number of iterations was set to 100, the number of reference individuals per population to 200, and the number of iterations per reference individual to 20. The degree of genetic differentiation between the samples in Serbia was also determined based on 5′ COI mtDNA and ITS2 rDNA in BAPS 6.0 using identical settings. Apart from this individual-based analysis, genetic structuring in Serbia was also investigated using population-based analysis of molecular variance (AMOVA) (Excoffier et al., Reference Excoffier, Smouse and Quattro1992) performed in Arlequin 3.11 software (Excoffier et al., Reference Excoffier, Laval and Schneider2005) using both markers. The mosquitoes sampled at different localities were pooled together as members of SII (urban) and SIII (rural) samples. Ten thousand permutations were used to determine the significance of variance components. Genetic differentiation, measured as Wright's F ST (Weir & Cockerham, Reference Weir and Cockerham1984), was also estimated among samples (urban and rural ecotypes of each locality were considered separately) using Arlequin 3.11, and the significance between each comparison pair was evaluated through 1000 permutation procedures.
Results
Identification and delimitation of Culex p. pipiens ecotypes at the transregional scale
The transregional review of the assays’ diagnostic utility encompassed 159 5′ COI mtDNA and 47 ITS2 rDNA sequences. The sequences yielded five 5′ COI mtDNA haplotypes (HI–HV) and 22 ITS2 rDNA alleles (A–V; fig. 1), with four and 45 variable positions in total, respectively. No stop codons were registered in the 5′ COI mtDNA sequences and HII–HV differed from HI in a single substitution each. Furthermore, using the CQ11FL assay, all the SII and SIII Cx. p. pipiens individuals sampled in Serbia were identified as pipiens ecotype (Supplementary Fig. S1).
Spanning network of COI mtDNA sequence haplotypes and ITS2 rDNA alleles derived from the transregional Cx. p. pipiens ecotype analysis. Each circle represents one haplotype/allele, with black and white coloration indicating DNA sequences found in molestus and pipiens ecotypes, respectively, while grey coloration shows sequences found in the Vojvodina Province. The length of the branches connecting circles is proportional to the number of base differences between the respective sequences.
The lineage sorting between Cx. p. pipiens ecotypes was incomplete since molestus and pipiens shared both 5′ COI mtDNA haplotypes (two haplotypes: HI and HII) and ITS2 rDNA alleles (one allele: allele A). Furthermore, with the most common sequences (5′ COI mtDNA HI and ITS2 rDNA allele A, 54/121 and 13/47 individuals, respectively) being shared between molestus and pipiens individuals (between-ecotype p distance = 0%), these markers cannot be used to distinguish the two ecotypes due to the lack of a barcode gap. Furthermore, HIV has not been previously associated with either Cx. p. pipiens ecotype outside Serbia. This problem is particularly pronounced in the ITS2 rDNA data set since six out of seven Vojvodina alleles have not been described before. Therefore, the discovery of novel, regionally private sequences further limits the utility of both molecular assays. We considered the possibility that the assay's diagnostic power may have been obscured by the use of different ecotype-defining criteria. Namely, molestus and pipiens mosquitoes whose sequences were harvested from GenBank were defined using one or more of the following criteria: autogeny, habitat type, distributional data, siphonal index of larvae and 5′ COI mtDNA and CQ11FL molecular assays (table 1). This is particularly noteworthy in the case of the distributional data, where Cx. p. pipiens mosquitoes in South Korea, Taiwan, and Australia were ‘by default’ treated as molestus ecotype since it is assumed that only this ecotype was introduced to urban environments in these regions (Farajollahi et al., Reference Farajollahi, Fonseca, Kramer and Kilpatrick2011). The two criteria which have been considered reliable are autogeny and CQ11FL identification. Yet, when only a subset of mosquitoes defined by those two criteria is considered, 5′ COI mtDNA is still not sufficiently diagnostic to delimit the ecotypes. Namely, HI was shared by pipiens mosquitoes of the UK, Austria, Serbia, and Russia, and molestus mosquitoes sampled in the UK and Austria. Similarly, HII was represented by both molestus mosquitoes of Russia and pipiens mosquitoes sampled in Italy. Finally, the utility of the assay could be compromised by the high morphological similarity of Cx. p. pipiens and Cx. torrentium (Shaikevich & Zakharov, Reference Shaikevich and Zakharov2010). Therefore, we compared a 5′ COI mtDNA sequence from GenBank (AM403477) with the haplotypes retrieved from our sample and found that the Cx. torrentium sequence differed in 13 substitutions from HI (the most common Cx. p. pipiens haplotype). Thus, it is unlikely that sample misidentification affected the results.
In the transregional BAPS analyses of genetic structuring, both 5′ COI mtDNA- and ITS2 rDNA-based estimations retrieved five as the most probable number of clusters (fig. 2). In the case of 5′ COI mtDNA, the overall pipiens sample originating from different countries was predominantly represented by the members of cluster 1 (around 80% of individuals, Supplementary Table S1), similarly to Serbian SII and SIII pipiens samples (around 90% of individuals). When it came to the molestus sample, it was mostly represented by members of cluster 2 (around 60%), but also by the members of cluster 1 and cluster 3 (22 and 20%, respectively). Furthermore, molestus and pipiens sample groups were not always assigned to different clusters (fig. 2). Namely, the majority of samples (11/17), consisting of molestus, pipiens, and hybrid mosquitoes were placed in a single cluster, showing no 5′ COI mtDNA genetic structuring corresponding to molestus/pipiens division (fig. 2). Even in the case of locally separated groups corresponding to Cx. p. pipiens ecotypes, Italian pipiens ecotype was assigned to a cluster where the majority of samples belonged to the molestus ecotype (4/5 groups). In the cases of Italy, Turkey, Germany, and Russia, molestus and pipiens ecotypes were indeed placed in alternative 5′ COI mtDNA clusters. Conversely, both ecotypes sampled in the UK and Austria were placed in the same 5′ COI mtDNA cluster.
Membership of the Vojvodina Province Cx. p. pipiens urban (SII) and rural (SIII) individuals (on the left) and Cx. p. pipiens ecotypes of the total sample (on the right), in a number of presumed clusters, as determined in BAPS using COI mtDNA and ITS2 rDNA loci. An individual's probability of belonging to one of the clusters is represented by different colors.
In terms of ITS2 rDNA-based differentiation, no overall pattern of genetic structuring was apparent, with both molestus and pipiens ecotypes highly admixed by members of different clusters (fig. 2, Supplementary Table S1). The alleles retrieved in our study were distributed in four groups: (1) private for molestus, (2) private for pipiens, (3) private for Serbian pipiens urban and rural samples, and (4) shared among molestus, pipiens, and Serbian pipiens SII and SIII mosquitoes (the most common allele, A).
Culex p. pipiens aboveground urban and rural differentiation in Serbia
In the intraregional analyses of genetic differentiation, 26 5′ COI mtDNA (around 900 bp-long) and 21 ITS2 rDNA (around 325 bp-long) sequences were recovered from the Vojvodina Province since not all sequencing reactions were successful. The sequences yielded four 5′ COI mtDNA haplotypes (H1–H4) and seven ITS2 rDNA alleles (a–g, fig. 3), with three and ten variable positions in total, respectively. Regarding 5′ COI mtDNA, H2, H3, and H4 differed from H1 in one substitution each, while ITS2 rDNA alleles diverged up to four positions from allele a. No stop codons were registered in the 5′ COI mtDNA sequences.
Spanning network of mitochondrial cytochrome c oxidase subunit I DNA (COI mtDNA) sequence haplotypes and internal transcribed spacer 2 of ribosomal DNA (ITS2 rDNA) alleles derived from the Vojvodina Province Culex p. pipiens populations. Each circle represents one haplotype/allele, with black and white coloration indicating DNA sequences found in urban (SII) and rural (SIII) samples, respectively. Each dash indicates a single base difference between the sequences.
Genetic differentiation between aboveground urban and rural Cx. p. pipiens samples on the territory of northern Serbia was assessed using 5′ COI mtDNA and ITS2 rDNA by implementing individual- and population-based approaches. Genetic structure analyses performed in BAPS based on 5′ COI mtDNA and ITS2 rDNA sequences retrieved four and three genetic clusters in Vojvodina, respectively (fig. 2). 5′ COI mtDNA genetic clusters corresponded to each haplotype registered in Vojvodina (H1–H4), while seven ITS2 rDNA alleles were clustered in three distinct groups: alleles a and b were grouped together, allele e was singled out, and the remaining alleles c, d, f, and g were recognized as a separate group. In terms of the average membership coefficient of urban and rural Cx. p. pipiens samples, the 5′ COI mtDNA results showed no apparent distinction between mosquitoes of the SII and SIII samples as both groups were predominantly represented by the members of clusters 1 (around 70%) and 2 (around 20%) (Supplementary Table S1). However, in the case of ITS2 rDNA, there was a notable degree of genetic differentiation, with SII mosquitoes predominantly represented by members of cluster 1 (67%) and SIII mosquitoes mainly represented by cluster 2 (78%). Furthermore, no cluster 1 members were registered in SIII sample, while no cluster 2 members were registered in SII sample. However, a portion of both samples’ genetic structure was represented by members of cluster 3 (33% and 22% with respect to SII and SIII ecotypes).
Similarly, slightly discordant results of genetic differentiation were registered when population-based AMOVA analyses were implemented. The DNA barcode markers, 5′ COI mtDNA and ITS2 rDNA, showed no genetic differentiation at varying hierarchical levels, both when the samples were divided into urban and rural groups and when they were considered together (ungrouped) (Supplementary Table S2).
Therefore, neither individual- nor population-based approaches of genetic differentiation showed concordance with the correlation between the habitat type where the mosquitoes were sampled in Vojvodina, except in the case of the individual-based approach using ITS2 rDNA sequences.
Discussion
The diagnostic utility of 5′ COI mtDNA and ITS2 rDNA assays in molestus/pipiens delimitation
Sequence-based molecular approaches to Cx. p. pipiens ecotype designation predominantly used 5′ end of cytochrome c oxidase subunit I mtDNA (5′ COI mtDNA) as a DNA barcode marker (Shaikevich, Reference Shaikevich2007, Reference Shaikevich2017; Shaikevich & Zakharov, Reference Shaikevich and Zakharov2010; Danabalan et al., Reference Danabalan, Ponsonby and Linton2012; Gunay et al., Reference Gunay, Alten, Simsek, Aldemir and Linton2015). This locus was firstly proposed by Shaikevich (Reference Shaikevich2007) as diagnostic for molestus/pipiens delimitation when underground, autogenous populations in a European part of Russia were found to differ in one fixed substitution from aboveground, anautogenous populations (HII vs. HI labeled herein, respectively). A subsequent analysis of Cx. p. pipiens ecotypic differentiation on the same territory (Shaikevich & Zakharov, Reference Shaikevich and Zakharov2010) used the entire COI mtDNA sequence (around 1550 bp) to define three haplotypes (labeled A, B and C) associated with pipiens ecotype and a single, fixed haplotype (D) characteristic of molestus ecotype. When only 5′ COI mtDNA is considered, previously defined haplotypes A and C correspond to our HI, while B and D haplotypes correspond to our HV and HII. Since both HV and HII differ from HI in a single substitution, it is evident that there is no threshold value in 5′ COI mtDNA genetic distance which separates the two ecotypes.
DNA barcoding underwent conceptual developments since its inception as a distance-based molecular tool, assimilating the character-based framework of traditional taxonomy and recognizing taxon boundaries based on discrete nucleotide substitutions (character states) within a DNA sequence (DeSalle et al., Reference DeSalle, Egan and Siddall2005). In this regard, 5′ COI mtDNA would still stand as a valid DNA barcode marker since a diagnostic transition of guanine to adenosine at the position 270 of 5′ COI mtDNA sequence used by Shaikevich & Zakharov (Reference Shaikevich and Zakharov2010) distinguished molestus from pipiens mosquitoes. However, in our study, the transregional analysis of 5′ COI mtDNA barcode marker was not diagnostic for molestus and pipiens designation, with HI and HII being shared between them. Haplotype HI, which was proposed as specific for pipiens ecotype (Shaikevich, Reference Shaikevich2007), was registered in molestus mosquitoes in the UK (Danabalan et al., Reference Danabalan, Ponsonby and Linton2012), Austria (Zittra et al., Reference Zittra, Flechl, Kothmayer, Vitecek, Rossiter, Zechmeister and Fuehrer2016), and Australia (Batovska et al., Reference Batovska, Blacket, Brown and Lynch2016), while HII was, apart from molestus mosquitoes, registered in pipiens mosquitoes in Italy (Di Luca et al., Reference Di Luca, Toma, Boccolini, Severini, La Rosa, Minelli, Bongiorno, Montarsi, Arnoldi, Capelli, Rizzoli and Romi2016).
The practical outcome of this finding challenges studies which exclusively used the 5′ COI mtDNA assay to identify which Cx. p. pipiens ecotype was present in the studied region. Namely, Gunay et al. (Reference Gunay, Alten, Simsek, Aldemir and Linton2015) and Shaikevich (Reference Shaikevich2017) used only the 5′ COI mtDNA assay to identify Cx. p. pipiens ecotypes, recognizing representatives of HI as pipiens and those of HII as molestus. However, their findings could easily be erroneous in the light of the presented results, and even epidemiologically irrelevant since the assay does not actually imply any strict correlation to physiological and behavioral traits. Furthermore, 5′ COI mtDNA was used for molestus/pipiens delimitation along with the CQ11FL assay. Although Danabalan et al. (Reference Danabalan, Ponsonby and Linton2012) recognized 5′ COI mtDNA as a more robust molecular assay due to the unexpected presence of molestus mosquitoes in aboveground habitats in the UK, our results support Di Luca et al.’s (Reference Di Luca, Toma, Boccolini, Severini, La Rosa, Minelli, Bongiorno, Montarsi, Arnoldi, Capelli, Rizzoli and Romi2016) conclusion that the CQ11FL assay is more accurate for Cx. p. pipiens ecotype identification.
Although ITS2 rDNA proved successful for distinguishing the closely related species of the genus Anopheles (see Batovska et al., Reference Batovska, Cogan, Lynch and Blacket2017), we found it to be inefficient in Cx. p. pipiens ecotype designation. Similar to 5′ COI mtDNA, this marker cannot be used either as a genetic distance- or character state-based DNA barcode marker for identifying molestus and pipiens ecotypes. Whether ITS2 rDNA can reliably distinguish between Cx. p. pipiens and the closely related Cx. quinquefasciatus remains questionable. Wang et al. (Reference Wang, Tu, Huang, Chen, Chen and Yeh2017) successfully distinguished between the two taxa in Taiwan using this marker, but Batovska et al. (Reference Batovska, Cogan, Lynch and Blacket2017) failed to do so in Australia. Probably owing to its high mutational rate, these results also point out spatial ITS2 rDNA variability: its utility as a DNA barcode marker to separate Cx. p. pipiens from Cx. quinqefasciatus varied between the continents and 20/22 alleles (table 1) in our study were private for the region where the mosquitoes were sampled.
In summary, neither 5′ COI mtDNA nor ITS2 rDNA were useful in distinguishing Cx. p. pipiens ecotypes. Given that their success was not latitude-dependent, these assays should not be used for the identification of molestus and pipiens mosquitoes, especially since a better-suited marker exists (CQ11FL). A further advantage of this assay based on length polymorphism is its ability to distinguish molestus and pipiens hybrids from pure ecotype mosquitoes. However, a combination of several markers seems to give the best resolution for identifying Cx. pipiens complex members. Namely, Cx. torrentium, Cx. quinquefasciatus, and Cx. pipiens should be discriminated based on the polymorphism in intron 2 of the gene encoding acetylcholine esterase 2 (ace-2; Smith & Fonseca, Reference Smith and Fonseca2004), while CQ11FL assay (Bahnck & Fonseca, Reference Bahnck and Fonseca2006), despite the possibility of recombination obscuring the identification (Bahnck & Fonseca, Reference Bahnck and Fonseca2006; Shaikevich & Vinogradova, Reference Shaikevich and Vinogradova2014), should then be used for distinguishing between molestus and pipiens ecotypes.
The admixture of Cx. p. pipiens urban/rural samples in Serbia and its epidemiologic consequences
Concerning the differentiation between urban and rural samples of Cx. p. pipiens in Serbia, neither 5′ COI mtDNA nor ITS2 rDNA indicated genetic differentiation between urban and rural Cx. p. pipiens samples in the population-based approach. Furthermore, we analyzed 306 individual genotypes from a previous study (Krtinić et al., Reference Krtinić, Francuski, Ludoški and Milankov2016) using allozyme loci, performing both individual- (STRUCTURE) and population-based (AMOVA) analyses and obtaining the same urban/rural admixture pattern (data not shown). Such a concordance of markers which differ in evolutionary rate and inheritance mode suggests an intense gene flow between SII and SIII samples. Our results are compatible with the findings of Honnen & Monaghan (Reference Honnen and Monaghan2018) who tested isolation by distance in Cx. p. pipiens mosquitoes of urban, peri-urban, and rural habitats in Germany using microsatellite loci, concluding that all the sampled individuals were derived from a single population with continuous gene flow through overlapping ranges. Furthermore, using the differences in electrophoretic mobility of hexokinase allozymes, Krtinić et al. (Reference Krtinić, Francuski and Milankov2014) previously recognized gene flow as a highly relevant microevolutionary process shaping genetic connectivity of urban and rural samples in Vojvodina Province. Therefore, although urban area exposes mosquitoes to the effects of air, water, and noise pollution (Honnen & Monaghan, Reference Honnen and Monaghan2018), gene flow still appears to keep the urban and rural habitats of Cx. p. pipiens genetically connected, as evidenced by both sequence-based (5′ COI mtDNA and ITS2 rDNA) and multilocus markers (allozymes, microsatellite loci). The future studies should investigate whether intense gene flow characteristic of the supposedly neutral markers overpowers divergent selection acting upon the loci involved in specific adaptations to urban and rural habitats.
A well-recognized issue of Cx. p. pipiens mosquitoes is the hybridization of the two ecotypes since it can lead to an alteration of epidemiologically relevant traits (Fonseca et al., Reference Fonseca, Keyghobadi, Malcolm, Mehmet, Schaffner, Mogi, Fleischer and Wilkerson2004). Namely, the hybrid offspring of ornithophilic pipiens ecotype and mamophylic molestus ecotype could act as the bridge vector promoting an exchange of avian and mammalian pathogens. In Serbia, such a scenario is supported by the registered outbreaks of West Nile virus (Popović et al., Reference Popović, Milošević, Urošević, Poluga, Lavadinović, Nedeljković, Jevtović and Dulović2013; Kemenesi et al., Reference Kemenesi, Krtinić, Milankov, Kutas, Dallos, Oldal, Somogyi, Németh, Bányai and Jakab2014), clearly demonstrating the ability of Cx. p. pipiens mosquitoes to transfer an avian pathogen to human hosts. Furthermore, hybrids were detected in lower latitudes: Portugal (Gomes et al., Reference Gomes, Sousa, Novo, Freitas, Alves, Côrte-Real, Salgueiro, Donnelly, Almeida and Pinto2009, Reference Gomes, Sousa, Vicente, Pinho, Calderón, Arez, Almeida, Donnelly and Pinto2013), Greece (Papa et al., Reference Papa, Xanthopoulou, Tsioka, Kalaitzopoulou and Mourelatos2013), and Italy (Di Luca et al., Reference Di Luca, Toma, Boccolini, Severini, La Rosa, Minelli, Bongiorno, Montarsi, Arnoldi, Capelli, Rizzoli and Romi2016), but also in higher latitudes: Austria (Zittra et al., Reference Zittra, Flechl, Kothmayer, Vitecek, Rossiter, Zechmeister and Fuehrer2016), the Netherlands (Reusken et al., Reference Reusken, de Vries, Buijs, Braks, den Hartog and Scholte2010), and Germany (Rudolf et al., Reference Rudolf, Czajka, Börstler, Melaun, Jöst, von Thien, Badusche, Becker, Schmidt-Chanasit, Krüger, Tannich and Becker2013), therefore supporting that gene flow between molestus and pipiens mosquitoes is not as restricted by latitude-dependent effect as thought previously. Finally, with gene flow keeping urban and rural samples connected, these effects of hybridization are not necessarily restricted to urban areas. Taking into the account and the possible spread of resistance alleles (between different urban areas, across rural areas between them), it is evident that studying microevolutionary processes which shape Cx. p. pipiens genetic structure remains an issue of high relevance for public health.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0007485319000105.
Acknowledgements
The authors would like to thank the Associate Editor Thierry Backeljau and three anonymous reviewers for their constructive criticism and helpful comments on earlier drafts of the manuscript. Lj.F., N.G., and V.M. are supported by the Ministry of Education, Science and Technological Development of Serbia (Dynamics of gene pool, genetic, and phenotypic variability of populations, determined by the environmental changes, No. 173012). B.K. is supported by Ciklonizacija d.o.o. Novi Sad.
