Avian Feeding Preferences of Culex pipiens and Culiseta spp. Along an Urban-to-Wild Gradient in Northern Spain

Mosquitoes (Diptera: Culicidae) are regarded as annoying biting pests and vectors of disease-causing agents to humans and other vertebrates worldwide. Factors that affect their distribution and host choice are not well understood. Here, we assessed the species abundance, community composition, and feeding patterns of mosquitoes in an urban-to-wild habitat gradient in northern Spain. Adult mosquitoes from four habitats (urban, periurban, rural, and wild) were collected by aspiration from mid-July to mid-September, 2019. Thirteen species were represented among the 268 specimens (132 females and 136 males) trapped, including six new records reported for the first time in the region. Culex pipiens was the most abundant species in all habitats except in the wild, where Culex territans was dominant. The highest mosquito diversity was recorded in the wild habitat [species richness (S) = 10 and Shannon/Margalef-Diversity Indices (H’/MI = 1.51/1.36)] and the lowest in the urban habitat (S = 3; H’/MI = 0.24/0.41). Blood-engorged specimens (n = 65) represented 49.2% of the total female collections. Eighty percent of the blood-meals (n = 52) were successfully identified based on cytochrome c oxidase I subunit (COI) DNA barcoding. Nine species of birds were identified in blood meals from the three ecological forms of Cx. pipiens (n = 48), Culiseta fumipennis (n = 3), and Culiseta morsitans (n = 1) collected along the four sampling habitats. Four dominant bird species were recorded in Cx. pipiens, i.e., Parus major (35.4%), Turdus merula (18.7%), Pica pica (18.7%), and Passer domesticus (10.4%). Despite the availability of dog and human hosts in the sampling sites located in the urban habitat, Cx. pipiens seemed to have a preference to feed on birds. Culiseta fumipennis blood-meal host records are reported for first time in Europe. These findings on mosquito blood-feeding preferences and habitat community changes will help to better understand vector-host associations and pathogen transmission paths.


INTRODUCTION
Mosquitoes (Diptera: Culicidae) are recognized as the most important arthropod group negatively impacting human and animal health worldwide. Numerous mosquito species serve as enzootic, bridge, and/or epidemic vectors of human pathogens. Although tropical countries are more exposed to mosquito-borne diseases, Europe is experiencing an increasing number of human cases. In southern Europe, local autochthonous transmission of important arbovirus diseases has been recorded [dengue virus (DENV), Chikungunya virus (CHIKV), and Zika virus (ZIKV)]. Other arboviruses [West Nile virus (WNV); Usutu virus (USUV); Sindbis virus (SINV); Tahyna virus (TAHV); Batai virus (BATV); Inkoo virus (INKV); and Snowshoe Hare virus (SSHV)] but also protozoans (Plasmodium spp.) and nematodes (Dirofilaria spp.) have been notified in the continent with variable distribution and effect on humans (Calzolari, 2016). Specifically in Spain, WNV, USUV, and DENV viruses have been detected in native and non-native mosquito species (Vázquez et al., 2011;Aranda et al., 2018). Zoonotic parasitic diseases, such as the canine heartworm and other filarioid nematodes, have also been found in native mosquito species (Bravo-Barriga et al., 2016). Protozoan parasites of avian species are also prominent in Spain, i.e., Haemoproteus and Plasmodium are commonly found in mosquitoes fed on birds (Gutiérrez-López et al., 2016;Martínez-de la Puente et al., 2016).
Among native mosquito species, the genus Culex has deserved great attention, due to their biting nuisance and potential as bridge vector for different pathogens (Vázquez et al., 2011;Martínez-de la Puente et al., 2016;Brugman et al., 2018). Mosquitoes of the Culex pipiens complex have a cosmopolitan distribution, and include several species, subspecies, forms, races, physiological variants, or biotypes (Becker et al., 2012). No consensus exists on the taxonomic status of the members of this complex. In Europe, Cx. pipiens is presented into three intraspecific forms: Cx. pipiens molestus, Cx. pipiens pipiens, and Cx. pipiens hybrids (Brugman et al., 2018). Hybrids between Cx. pipiens and Cx. quinquefasciatus have also been reported from the Mediterranean Basin (Shaikevich et al., 2016). Both Cx. pipiens forms (pipiens and molestus) and their hybrids are present in the Iberian Peninsula (Martínez-de la . Cx. pipiens has a good tolerance to human-altered environments and is an opportunistic feeder on birds and mammals, including humans (Brugman et al., 2018). Another native genus, namely Culiseta, includes some species that have been moderately implicated in the transmission of important arboviruses, such as USUV and WNV, in Europe (Martinet et al., 2019). Nevertheless, up to date no report of virus transmission has been noticed in Spanish mainland.
Mosquitoes inhabit diverse types of habitats in urban, periurban, and wild environments, where they undertake hostseeking activity and search for suitable substrates for the development of their progeny. Several factors are involved in the mechanisms that affect spatial-temporal distribution and host selection. The urban-to-wild gradient is composed of distinct environmental habitats that influence distribution, diversity, and abundance of mosquito species . In Europe, urbanization processes created suitable habitats for a small number of anthropophilic species, mostly Cx. pipiens (Becker et al., 2012), and Aedes albopictus (Li et al., 2014). Conversely, permanent water bodies in natural environments serve as favorable breeding sites for a wide range of species, such as Anopheles maculipennis s.l., Aedes vexans, Aedes sticticus, Aedes caspius, and Aedes detritus (Medlock and Vaux, 2015). The elucidation of how mosquito communities change along urban-to-wild habitats provides valuable information to characterize the environmental and ecological processes that significantly define the mosquito populations and hostrelated interactions.
Feeding patterns of hematophagous arthropods are also a critical component of the ecological cycles of pathogen transmission. Identifying the host affinity of mosquito vector species under field conditions is particularly important given the role of wild and domestic animals as reservoir hosts of several important arboviruses (Failloux et al., 2017). In addition, since not all avian species are competent hosts of arboviruses (Kilpatrick et al., 2007), the identification of blood meals is a first step for the detection of key reservoir species. Consequently, an enhanced understanding of the factors involved in the distribution and hostchoice of mosquitoes, particularly disease vectors, would improve the efficiency of surveillance and control programs .
Mosquito surveillance and screening efforts have been focused on the southern and eastern regions of Spain (Muñoz et al., 2011;Ferraguti et al., 2016;Martínez-de la Puente et al., 2016) while the northern regions have been comparatively overlooked due to their less favorable climatic conditions (harsh winters and short summers) for the proliferation of mosquitoes, bites and their diseases (Eritja et al., 2005;Ruiz-Arrondo et al., 2019). However, climate change and anthropogenic factors are beginning to affect the distribution and incidence of mosquito vectors and their pathogens, an example being the recent arrival and establishment of the exotic species Ae. albopictus  and Ae. japonicus in northern Spain (Eritja et al., 2019).
The main goal of the present study was to ascertain the species diversity including the screening of potential exotic species, abundance, and feeding patterns of mosquito species along an urban-to-wild gradient in northern Spain. In addition, we determined the success of identification of female mosquitoes based on digestion status, and the gonotrophic stages distribution. Finally, the efficacy of active suction trapping as monitoring tool is discussed.

Study Area and Sampling Sites
This study was conducted in the municipality of Vitoria-Gasteiz (ca. 260,000 inhabitants), capital of the Autonomous Community of the Basque Country (northern Spain). Four types of habitats, located 1.0, 2.5, 4.5, and 12.0 km from the city center, were selected following a gradient from a highanthropogenic site (urban) to a highly conserved natural site (wild). The human population density for each habitat was an estimation calculated from demographic data from 2018 ( Gobierno Vasco-Eusko Jaurlaritza, 2019) and the final selection of the study site was based on environmental characteristics ( Figure 1A). The urban habitat (42.853397, -2.682756, and ca. 860 habitants/km 2 ) was a densely populated area with buildings, schools, and multiple gardens containing hedges of Prunus cerasifera and Chamaecyparis lawsoniana, and several tree species, including Aesculus hippocastanum, Magnolia grandiflora, Liriodendron tulipifera, and Paulownia tomentosa. The periurban habitat (42.855776, -2.710002, and ca. 150 habitants/km 2 ) consisted of an industrial environment mixed with buildings, wastelands and rows of Quercus ilex and Fraxinus excelsior. Both urban and periurban settings had sewer systems in the vicinity. The rural habitat (42.841991, -2.717480, and ca. 40 habitants/km 2 ) comprised of a green patch (with paths, small lagoons, and diverse vegetation) surrounded by country homes and agricultural fields. The wild habitat (42.890751, -2.532650) where human activities were very limited or nonexistent, was located in an ornithological park at the edge of a large water reservoir (swamp) which contained aquatic plants, thick vegetation, grasslands, and small woodlands.

Mosquito Collection and Host Counting
We employed a battery-powered Prokopack aspirator (John W. Hock, Florida, United States) to collect outdoor-resting mosquitoes across the four habitats. This quiet, low-cost, and easy handling aspirator allowed working in both open and densely vegetated areas without disturbing or attracting attention of humans and animals. Collections were conducted between 9:00 a.m. and 12.00 p.m. along a 100 m line transect, and lasted approximately 30 min per habitat. Each habitat was sampled fortnightly from mid-July to mid-September 2019, giving a total of five sampling periods. Transects were selected according to the presence of natural resting sites, such as humid and shadow shelters (tree-lines, hedgerows) and corners in concrete walls. Insect collection was carried out by passing the aspirator up and down and from side to side, so that the entire surface was sampled. Insects escaping from their resting place were also collected by aspirating them out of the air. Collection cups that became full or clogged were replaced as necessary during each collection event. Immediately after collection, the cups were stored at −30 • C until further taxonomical and molecular analyses.
The number of hosts encountered in each sampling transect was scored concurrently to the mosquito aspiration along each sampling transect. In the case of mammals (i.e., humans and dogs), count was done at the beginning, middle and end of each transect, during 5-min point counts with a limited observation radius (30 m), unless the features of the landscape minimized the visual capacity of the operator. For avian species, counting was performed in a similar way by determining the number of individuals that were seen either walking, perched on vegetation/trees or took off after the walking activity of the operator within the visual radius (30 m). In this regard, no determination was done at species level. Those birds flying up in the sky were excluded from the counting ( Figure 1B).

Morphological Identification
Prior to taxonomical identification, specimens were separated by sex, and females sorted by their physiological status based on abdomen condition (blood-fed, gravid, and unfed). Species identification was done using appropriate taxonomic keys (Schaffner et al., 2001;Becker et al., 2010) based on morphological features of females and male genitalia. In order to identify blood-meal host species, abdomens of blood-engorged females were removed and individually stored in 2 ml screw-top vials after being classified according to the scale of Sella (II-VII) by observation of the degree of digestion of the blood and developing eggs (gravid stage = Sella stage VII) (Martínez-de la Puente et al., 2013).

Molecular Identification of Blood-Fed Mosquitoes and Their Hosts
Cytochrome c oxidase I subunit (COI) DNA barcoding was used to identify damaged or morphologically indistinguishable mosquito specimens (n = 2), as well as vertebrate host species of gravid (n = 11), and blood fed females (n = 65) at the Centre for Biodiversity Genomics, University of Guelph (Guelph, ON, Canada). Standard barcoding methods were used for mosquito identification (Hebert et al., 2003) while a modified version that accounts for DNA degradation, as well as the possibility of mixed feeding, was used for blood-meal vertebrate host identification. In all cases, abdomens were removed using sterile forceps and placed individually in 2-ml screw-tubes. DNA extraction was performed by immersing the abdomen in a lysis buffer containing proteinase-K followed by incubation at 56 • C for 18 h. Bloodengorged abdomens were purposefully disrupted during transfer into screw-tubes to facilitate DNA extraction from the blood meal. DNA purification employed a glass fiber-based bind-washelute method (Ivanova et al., 2006). Briefly, the lysate was mixed with two volumes (120 µL) and a binding mix, and then applied to a well of a 96-well silica membrane plate. The plate was centrifuged at 5000 g for 5 min, and washed two times. Traces of wash buffer were removed by incubating the plate at 56 • C for 30 min, after which 40 µL of elution buffer was applied directly to the membrane and allowed to incubate at room temperature for 1 min. DNA was eluted from the membrane into a clean plate via centrifugation at 5000 g for 5 min.
For mosquito identification, the entire barcode region of COI was amplified using the primers C_LepFolF + C_LepFolR (Hernández-Triana et al., 2014) and the products were used directly for Sanger sequencing (Hajibabaei et al., 2006). For and VR1i_t1 (CAGGAAACAGCTATGACTAGACTTCTGGGTGICCIAAIAA ICA). The forward primers were designed to bind to vertebratebut not mosquito -DNA, so that only vertebrate blood meal host DNA is amplified. These primers targeted a 185 bp fragment of the COI barcode region, allowing them to amplify highly degraded blood meal DNA. The thermocycling profile consisted of initial denaturation at 95 • C for 2 min, followed by 60 cycles of 95 • C for 40 s, 56 • C for 40 s, and 72 • C for 30 s, and a final extension at 72 • C for 5 min. After confirming amplification on a 2% E-gel (Invitrogen), the PCR products were diluted two-fold with sterile water and used as a template for a second round of PCR. This second PCR (PCR2) served to index the PCR1 amplicons with IonXpress UMI tags and Ion Torrent sequencing adapters. The M13F and M13R tails of the PCR1 primers served as universal binding sites for PCR2. PCR2 reactions consisted of the same components as PCR1. The PCR2 thermocycling consisted of an initial denaturation at 95 • C for 2 min, followed by 5 cycles of 95 • C for 40 s, 45 • C for 40 s, and 72 • C for 30 s, and then 35 cycles of 95 • C for 40 s, 51 • C for 40 s, and 72 • C for 30 s, and a final extension at 72 • C for 5 min. The PCR2 products were pooled in equal volumes and purified using SpeedBeads (Sigma Aldrich) by incubating 400 µL of pooled PCR2 product with 400 µL of beads. The DNA-bead mixture was incubated at room temperature for 8 min, after which the beads were collected on a magnetic rack. The bead pellet was washed three times with 1 mL of freshly prepared 80% ethanol and air dried for 10 min. DNA was released from the beads by mixing the pellet with 200 µL of sterile water for 1 min at room temperature. The beads were collected on a magnetic rack and 180 µL of the supernatant was transferred to a clean 1.5 mL tube. The purified library was quantified using a Qubit 2.0 fluorometer, adjusted to 26 pM with sterile water, and loaded onto an Ion Chef automated platform (Thermo Scientific) following manufacturer's instructions. Sequencing was performed using 400 bp chemistry on an Ion Torrent S5 (Thermo Scientific) using a 530 chip. Raw sequence reads were automatically demultiplexed following sequencing using the Torrent Browser sequencing software. The demultiplexed reads were further processed using standalone bioinformatics tools (see section "Bioinformatic, sequence information and statistical analysis" below for details). The efficacy of the technique to amplify both mammalian and avian hosts was validated prior to the analysis of the samples by testing the specificity of the primers by the inclusion of positive (DNA extracted from blood-fed insects known to have fed on certain vertebrate hosts), and negative controls. The amplified products were sequenced to confirm that they originated from the target vertebrate species and not the insect.

Molecular Identification of Sibling Species and Forms
The identification of the members of the Anopheles maculipennis complex (n = 5) and blood fed/gravid Cx. pipiens (a total of 64 out of 76 gravid/blood-fed female mosquitoes submitted to the molecular methods described in the previous section) ecological forms were carried out from DNA extracted from the thorax of individual specimens in the Animal Health Department of NEIKER (Spain). In both cases, genomic DNA extraction was carried out by NZY Tissue gDNA isolation kit (NZYTech, Lisboa, Portugal). Identification of the specimens of the An. maculipennis species complex was carried out by a PCR-RFLP assay targeting polymorphisms in the Internal Transcribed Spacer 2 (ITS-2) of the ribosomal DNA (Vicente et al., 2011). Amplicons were digested first with HhaI and secondly with HpaII restriction enzymes (Invitrogen/Thermo Scientific, Vilnius, Lithuania). Identification of the different Cx. pipiens female forms (Cx. pipiens pipiens, Cx. pipiens molestus, and its hybrids) was performed by PCR amplification of the flanking region of the CQ11 microsatellite, following the protocol already described (Bahnck and Fonseca, 2006). Amplicons from both PCR methods were separated by 2% agarose gel electrophoresis using GelRed Nucleid Acid Stain (Biotium, Hayward, CA, United States) as staining solution, and with a 100 bp DNA ladder (Thermo Scientific, Vilnius, Lithuania) as a molecular weight marker.

Bioinformatic, Sequence Information and Statistical Analysis
Raw sequence reads were filtered (based on a minimum length of 100 bp and minimum quality of QV20) and cleaned (removal of adapter and primer tails) to ensure that only high-quality reads were included in the final dataset. Reads passing these filters were clustered into operational taxonomic units (OTUs) with 97% identity and a minimum of 10 reads per OTU. Each OTU sequence was compared using BLAST tool 1 and BOLD Systems 2 , which are composed of global vertebrate COI sequences. Identifications were considered valid only if the query sequence matched a reference sequence with at least 95% identity with 100 bp of coverage between the queried sequence and the reference sequence. Furthermore, mosquito identifications were manually vetted in order to ensure that they conformed with local fauna sighting records.
Detailed specimen records and sequence information (including trace files) were uploaded to the Barcode of Life Database (BOLD-http://www.boldsystems.org) and can be found within the Working Group 1.4 Initiative "Human Pathogens and Zoonoses" container "MCBCS-Surveillance of mosquitoes and Culicoides in the Basque Country, Spain." The Digital Object Identifier (DOI) for the publicly available projects in BOLD is doi:dx.doi.org/10.5883/DS-MQBMBC. All generated sequences have been submitted to GenBank (accession numbers: MT519609-MT519682).
The diversity indices of Shannon and Margalef (H'/MI), and species richness (S) were used as a measure of mosquito community heterogeneity. Host counting was obtained from the sum of the five trapping periods. Both mean abundances of mosquitoes within habitats and gonotrophic stages were compared by non-parametric Kruskal-Wallis test followed by Mann-Whitney U test pairwise comparisons (adjusting significance levels). All tests were conducted with IBM SPSS statistics v 23.0 software package.

RESULTS
A total of 268 mosquitoes (132 females and 136 males; mean ± SEM = 13.4 ± 2.9 specimens/habitat) belonging to four genera and 13 species were identified during the five sampling periods (

Mosquito Community Composition and Abundance
Cx. pipiens was the most predominant species (200 specimens, 74.6% from the total collections), followed by Cx. territans (32, 11.9%), Cs. longiareolata (11, 4.1%), and other less frequent species (25, 9.3%; Table 1 and Figure 2A). Cx. pipiens was the only species collected along the four habitats and the most predominant in all sampled settings except in the wild habitat, where Cx. territans was the dominant species. Despite mosquito abundance being greater in the urban habitat, the difference was not significant (χ 2 = 6.4, df = 3, and p = 0.094). Higher mosquito S and biodiversity was recorded in the wild habitat (S = 10, H'/MI = 1.51/1.36) compared to the other habitats (S ≤ 4; H'/MI ≤ 0.40/0.75; Table 1).
Molecular identification of blood-fed/gravid females by means of PCR (COI gene) and sequencing, confirmed the identity of Cx. pipiens specimens, thus this species is referred in the text as Cx. pipiens (100% match identity to bank identification number in BOLD: AAA4751, accession number GU908075). Note that molecular identification allowed the separation of Cx. torrentium females from the sibling Cx. pipiens females.
The pattern of Cx. pipiens ecological forms was successfully determined in 89.1% (57/64) of the blood-fed/gravid Cx. pipiens analyzed, showing the presence of the three ecological forms in different proportions [Cx. pipiens pipiens (n = 50; 87.7%; in urban, periurban and rural habitats); Cx. pipiens molestus (n = 3; 5.3%; in urban and rural habitats); and Cx. pipiens hybrids (n = 4; 7.0%; solely in urban habitat)]. Cx. pipiens pipiens was predominant in all the habitats but the wild, in which no Cx. pipiens blood-fed specimens were collected.
Morphologically similar members of the subgenus Culiseta were separated into three species (Cs. fumipennis voucher specimen = 100% identity, accession number KM258139; Cs.

Analysis of Gonotrophic Stages
The distribution of the gonotrophic stage categories for each habitat and species is represented in Figures 1C, 2B, respectively. Overall, blood-engorged females represented the highest percentage of the catches, accounting for 49.2% (65/132) of the total (Figure 3A). However, no statistically significant differences were detected in the abundance of the three gonotrophic stages (χ 2 = 0.002, df = 2, and p = 1.2) ( Figure 3A). In contrast, significant differences were found in the number of blood-fed (χ 2 = 10.64, df = 3, and p = 0.014) and unfed females (χ 2 = 11.84, df = 3, and p = 0.008) within habitats (Supplementary Table 2 and Figure 3B). The statistical analysis revealed that significantly more bloodengorged females were aspirated in the urban habitat, while unfed individuals constituted the dominant gonotrophic stage in the wild habitat. Regarding gravid females, no differences were detected among sampling sites.

Host-Derived Blood-Meals Identification
Blood-meal host identification was successful for 80.0% (n = 52) of the 65 engorged females analyzed in this study. All the stages within the scale of Sella (II-VI) yielded high number of successful identifications regardless their blood-digestion degree [II = 9/13 FIGURE 3 | Distribution of gonotrophic stages in female mosquitoes in the four sampling habitats (A) Total frequency of gonotrophic stages (blood-engorged, gravid, and unfed). (B) Abundance (Mean + SEM) of gonotrophic stages per habitat. Columns with different letters (a,b) within each gonotrophic stage are statistically different, while "ns" denotes non-significant differences in the total number of females between gonotrophic stages (Kruskal-Wallis test followed by Mann-Whitney U test).
Culiseta spp. (n = 4) fed upon the European stonechat Saxicola rubicola, the song thrush Turdus philomelos and P. major for Cs. fumipennis, and T. merula for Cs. morsitans ( Figure 1D). Blood-engorged specimens from three other species (Cx. hortensis, Cs. litorea, and An. claviger s.l.) tested negative to amplification. Mixed host feeding was not detected in any of the analyzed specimens.

DISCUSSION
The results of this study revealed that Cx. pipiens mosquitoes were the dominant species in urban, periurban, and rural habitats, where they showed a preferred tendency to feed on avian hosts. This feeding pattern was clearly observed in the urban setting, where a higher number of blood-fed Cx. pipiens females were trapped. Taking into account the predominance of nonavian hosts (humans and dogs) in the urban environment, it is highly remarkable that the blood-feeding pattern recorded for Cx. pipiens was strongly biased to bird hosts. In contrast, the wild habitat did not yield blood-fed specimens of the two predominant Culex species (i.e., Cx. pipiens and Cx. territans), though blood-fed Culiseta spp. were demonstrated to have fed also on bird hosts.
In this study, neither human nor dog sequences were obtained in blood from Cx. pipiens mosquitoes from the urban habitats despite the relative high abundance of these hosts (64% and 26%, respectively). Similarly, humans were highly present in periurban and rural habitats (70% and 40%, respectively). On the basis of relative host abundance measures and mosquito diet, our results suggest a prevailing affinity of Cx. pipiens to feed on avian hosts in these habitats, however, it is not possible to definitively prove that the mosquitoes were not feeding on certain mammalian hosts (i.e., a failure to detect a host could be due to DNA degradation, primer binding issues for certain taxa, etc.). In Europe, the ornithophilic preference of Cx. pipiens has been observed in some studies Roiz et al., 2012;Gomes et al., 2013;Radrova et al., 2013;Brugman et al., 2017), however, other authors have also highlighted the importance of mammals as a relevant feeding source (Alcaide et al., 2009;Muñoz et al., 2011;Rizzoli et al., 2015;Börstler et al., 2016;Martínez-de la Puente et al., 2016). Both human and non-human vertebrates (i.e., Capreolus capreolus, Sus scrofa, Lepus europaeus, Canis lupus familiaris, Felis silvestris, Equus ferus caballus, Bos taurus, Oryctolagus cuniculus, Rattus rattus, and some reptiles) from Cx. pipiens have been reported in variable proportions depending on the study and regardless of the pipiens forms (Alcaide et al., 2009;Muñoz et al., 2011;Rizzoli et al., 2015;Börstler et al., 2016;Martínez-de la Puente et al., 2016). Five out of the eight species of birds identified in the blood-meals from the urban setting were included in the bird census undertaken by ornithologists in the same park (Ayuntamiento de Vitoria-Gasteiz, 2019). They noted Columbia livia and P. domesticus as the most frequent birds, followed by Apus apus, P. pica, Cyanistes cauruleus, P. major, T. merula and, to a lesser extent, three other species. This marked feeding tendency of Cx. pipiens to feed on T. merula and P. domesticus was also observed in other European studies Roiz et al., 2012;Gomes et al., 2013;Rizzoli et al., 2015). It is noteworthy that synanthropic and abundant bird species such as the rock pigeon C. livia could have been underestimated, since their abundance in the urban setting did not reflect the proportion of blood meals taken, which is consistent with other studies carried out in Europe and United States (Rizzoli et al., 2015;Kothera et al., 2020). Analyses of gravid/blood-fed Cx. pipiens showed that they fed mostly on birds regardless of the ecological form, which is in agreement with similar studies carried out in Spain, United Kingdom, and Portugal (Gomes et al., 2013;Martínezde la Puente et al., 2016;Brugman et al., 2017). Consequently, mosquito populations that mostly feed on birds would have a larger capacity to amplify WNV due to the fact that humans and other mammals do not support viremias high enough to infect mosquitoes (Platt et al., 2007).
In the wild habitat, the few blood-meals recorded corresponded to bird hosts from Cs. fumipennis and Cs. morsitans, which is in agreement to the relative bird density observed (89% of visible host). According to literature, Culiseta species feed on birds and mammals, including humans (Brugman, 2016), but in United States they have also been found feeding on reptiles (Blosser et al., 2017). To best of our knowledge, the Cs. fumipennis blood-meal identification represents the first report of the feeding pattern of this species in Europe. In Germany and United Kingdom, Cs. morsitans have been recorded feeding on birds but also on mammals, including humans (Börstler et al., 2016;Brugman et al., 2017). Further studies are necessary to establish solid conclusions about the feeding preferences of both species.
According to our findings, the highest S was found when the distance was increased from the urban sampling site, with the wild site harboring the highest number of species. This relatively low S is tightly linked to the negative effect of the urbanization processes on the diversity of mosquito communities . Non-disturbed wild habitats, in contrast, are constituted by a plethora of suitable heterogeneous microhabitats, leading to a favorable scenario for enhancing mosquito diversity. Environmental factors including land use, vegetation, and hydrological characteristics are known to affect mosquito abundance and community composition . Among them, the spatial heterogeneity gradient is reported as one of the most important key indicators with relevant influence in mosquito species diversity and distribution (Overgaard et al., 2003;Leisnham et al., 2014). For instance, urbanized environments are organized in small patches (buildings, green areas, and streets) which are highly structured with many interfaces between them. Conversely, natural environments habitats are structured in larger patches and less interfaces, resulting in high species diversity. As stated by McDonnell and Pickett (1990), suburban areas are usually prone to congregate higher S and diversity. Indeed, we observed that fewer species were found in urbanized environments compared to rural/wild environments. However, the number of catches did not follow a well-defined pattern, as the highest abundance was observed in the anthropogenic habitat. This contrasts with previous reports (Ibañez-Justicia et al., 2015;Ferraguti et al., 2016;de Valdez, 2017;Honnen and Monaghan, 2017) but the collection method employed was different among surveys (i.e., CDC-mini light traps, BG-Sentinel traps baited with CO 2 , and aspiration). One possible explanation to the observed discrepancies is that mosquitoes could be more concentrated on the sampling site of the urban setting where fewer resting sites were available to use compared to wild areas, which in turn offer a wide variety of features for mosquito activity (Burkett-Cadena et al., 2013). These discrepancies are the result of the presence of substantial numbers of Cx. pipiens (mostly the pipiens form) and Cs. longiareolata, highly anthropophilic and container-breeding species. Both species find in urban areas their preferred oviposition sites, i.e., artificial habitats such as urban sewage and ground water systems (Bueno Marí, 2010), and other temporary water bodies (gardens, terraces, etc.). Cx. pipiens abundance decreased as the distance from the urban setting increased, being replaced by other species, such as Cx. territans, Cs. fumipennis, and An. maculipennis s.s. The lower abundance of mosquitoes in the natural habitat might have been caused by the presence of a dense vegetated barrier between the sampling transect and swamp land, which could have prevented the dispersion of mosquitoes from breeding sites. Indeed, mosquito community composition and habitat selection are highly variable among countries (Möhlmann et al., 2017). The combination of various collection methods would allow a better representation of the mosquito and host community fauna, however, the use of other routine traps, i.e., suction traps are inappropriate as they are subject to vandalism and logistic constraints in public urban/periurban areas.
Identification rates of host-derived blood-meals were not consistent with Sella's scale, since we did not observe a marked reduction in the amplification success as the degree of digestion in the blood increased (Martínez-de la Puente et al., 2013;Brugman et al., 2017;Santos et al., 2019). These discrepancies may be attributable to the technique implemented which has been improved for amplification of degraded DNA. In addition, the technique used did not allow the identification of host-derived blood meals from gravid females, thus limiting the determination of host-feeding patterns exclusively to specimens containing signs of blood on their abdomens.
Our study has demonstrated the effectiveness of the sampling performance of the Prokopack aspirator for collection of outdoor-mosquitoes in comparison to other superior models (i.e., CDC-BP aspirator) (Vazquez-Prokopec et al., 2009;Maia et al., 2011). The Prokopack aspirator is simple to use, easier to maneuver, cheap, and lightweight. In addition, aspiration yielded larger numbers of blood-engorged specimens when compared to other standard methods (Maia et al., 2011;Brugman et al., 2017;Santos et al., 2019), which makes this technique suitable for studies of host-feeding preferences and/or pathogen detection. Thus, captures of blood-fed females are usually scarce (0.5-3.0%) using CDC traps Radrova et al., 2013). Moreover, direct aspiration from mosquito resting sites provides a more representative information and estimation of physiological condition, richness, abundance, sex ratio, and age structure (Silver, 2008) in comparison to other trapping methods. Conversely, the main constraint is the damage inflicted on the specimens after being aspired and retained in the trap mesh for a long time, which frequently hampers proper identification of females. The differences in the gonotrophic stage frequencies (i.e., proportion of blood meals taken by mosquitoes) of female mosquitoes among habitats may be influenced by multiple factors difficult to be analyzed, such as the composition, abundance, and availability of vertebrate hosts, which in turn is linked to the landscape anthropization level and ecological conditions. Other variables such as collection methods (e.g., CDC traps baited with CO 2 capture preferentially host-seeking females) and availability of breeding and resting places with different accessibility (e.g., in the urban setting lower resting places were available compared to the wild) (Maia et al., 2011;Barrera et al., 2012;Muñoz et al., 2012;Ferraguti et al., 2016), could explain the high proportion of blood-fed and unfed mosquitoes in urban and wild habitats, respectively. For example, in our study collections might be biased toward mosquitoes fed on avian hosts as most sampling sites were bird resting sites (bushes, tree hedges, etc.). A more controlled design of the sampling sites selection and a more accurate host counting method might help to solve these issues.
Faunistic studies focused on mosquito communities in northern Spain are scarce and mostly devoted to the detection of those species with relevant medical interest Goiri et al., 2020). Here, we report for the first time the presence of six new species in the region, increasing from 14 to 20 the total number of species present in the Basque Country (Bueno-Marí et al., 2012;González et al., 2015González et al., , 2020Goiri et al., 2020), including new records of Anopheles maculipennis s.s., a major vector of several pathogens in Europe (Kampen et al., 2016). Interestingly, in nearby territories, only An. atroparvus had been identified within this species complex . This study also contributes to the differentiation of Culiseta sibling species by barcoding. The separation of some female species of Culiseta remains challenging by both DNA barcoding and morphological features (Ruiz-Arrondo et al., 2020). The female sequences obtained had enough genetic divergence to be separated into individual OTUs (data not shown).
In conclusion, Cx. pipiens was the most widespread species along the studied habitats according to the aspiration-based sampling methodology. Cx. pipiens, Cs. fumipennis, and Cs. morsitans exhibited an overwhelming affinity for avian hosts. These results on mosquito blood meals and community composition contribute to a better understanding on the transmission risk of pathogenic agents of medical and veterinary importance. The study also highlights the importance of implementing monitoring programs that include mosquito species and host-feeding surveillance but also vector ecology.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
MG and AG-P conceived and planned the experiments. MG carried out the field experiments and morphological identification of mosquitoes. SL performed statistical analyses. SP, PH, and LH-T carried out the molecular analysis. FG made molecular identification of sibling species and forms. MG, PA-E, SL, and IR-A wrote the manuscript. AG-P and FG reviewed and edited the final version of the manuscript. All authors provided critical feedback and helped shape the research, analysis and manuscript.

ACKNOWLEDGMENTS
The authors thank "Unidad de Anillo Verde y Biodiversidad del Ayuntamiento de Vitoria-Gasteiz" and "Centro de Estudios Ambientales del Ayuntamiento de Vitoria-Gasteiz" for their approval and consent to perform this study at bird-relevant observation spots. Bird pictures were taken with a Canon (EOS-1Ds Mark II + 600 mm f4) by the ornithologist Mr. Brain Webster (Instituto Alavés de la Naturaleza) and used here under his consent and approval.