Ecology of West Nile Virus in the Danube Delta, Romania: Phylogeography, Xenosurveillance and Mosquito Host-Feeding Patterns

The ecology of West Nile virus (WNV) in the Danube Delta Biosphere Reserve (Romania) was investigated by combining studies on the virus genetics, phylogeography, xenosurveillance and host-feeding patterns of mosquitoes. Between 2014 and 2016, 655,667 unfed and 3842 engorged mosquito females were collected from four sampling sites. Blood-fed mosquitoes were negative for WNV-RNA, but two pools of unfed Culex pipiens s.l./torrentium collected in 2014 were tested positive. Our results suggest that Romania experienced at least two separate WNV lineage 2 introductions: from Africa into Danube Delta and from Greece into south-eastern Romania in the 1990s and early 2000s, respectively. The genetic diversity of WNV in Romania is primarily shaped by in situ evolution. WNV-specific antibodies were detected for 19 blood-meals from dogs and horses, but not from birds or humans. The hosts of mosquitoes were dominated by non-human mammals (19 species), followed by human and birds (23 species). Thereby, the catholic host-feeding pattern of Culex pipiens s.l./torrentium with a relatively high proportion of birds indicates the species’ importance as a potential bridge vector. The low virus prevalence in combination with WNV-specific antibodies indicate continuous, but low activity of WNV in the Danube Delta during the study period.


Introduction
Emerging or re-emerging mosquito-borne viruses (moboviruses) are of growing concern in Europe [1]. Several moboviruses circulate on the European continent [2]. Thereby, West Nile virus Thus, in order to get comprehensive insight into the ecology of WNV in the DDBR, the mosquito fauna was studied in a longitudinal surveillance program over three years. Molecular assays were applied to (i) screen for WNV infections in mosquitoes, analyze the evolutionary mechanism of the virus and its dispersal patterns in Europe, in particular in Romania and the DDBR, (ii) detect WNV-specific antibodies in the blood meals from horses, dogs, humans and birds and (iii) identify potential vector species by analyzing the host-feeding patterns of the blood-fed mosquitoes.

Materials and Methods
Mosquitoes were collected within a longitudinal arbovirus surveillance program between 2014 and 2016 at two sampling sites in a rural/urban environment (Letea, Sulina) and two near-natural sampling sites (Dunărea Veche and Lake Ros , ulet , ) in the DDBR. Each year, on average, every tenth day between April and September, three to four (2014,2015) or one (2016) carbon dioxide-baited Heavy Duty Encephalitis Vector Survey trap(s) (Bioquip Products Inc., CA, USA) were installed at each site. A detailed description of the collection sites can be found in Török et al. [9]. The DDBR Authority issued research permits (9/25.04.2014, 10692/ARBDD/25.04.2014; 7717/ARBDD/28.04.2016, 11/28.04.2016). The collected specimens were transported on dry ice, stored in the freezer and identified by morphology on chill tables using a stereomicroscope (Olympus SZX12, Tokyo, Japan) [39]. Blood-fed mosquitoes were separated from unfed specimens. Furthermore, morphologically identified Culex pipiens specimens were typed to species level (Cx. pipiens pipiens f. pipiens, Cx. pipiens pipiens f. molestus or Cx. torrentium) using a molecular assay [40].
For the WNV screening, mosquito pools between 1 and 250 specimens were pooled per sampling site and sampling date. Mosquitoes were put in 2 mL safe-lock tubes (Eppendorf, Hamburg, Germany) or 50 mL centrifuge tubes (Sarstedt, Nümbrecht, Germany) with zirconia beads (2 mm, BioSpec Products, Bartlesville, OK, USA) and 0.5 or 3 mL chilled high-glucose (4.5g/L) Dulbecco's modified Eagle's medium (DMEM) (Sigma-Aldrich, St. Louis, MO, USA). Mosquitoes were homogenized in a TissueLyser or TissueLyser II (Qiagen, Hilden, Germany) for 2 min at 30-50 Hz. The suspension was clarified by centrifugation for 1 min at 8000 rpm and 4 • C. RNA was extracted with a KingFisher Flex 96 Deep-Well Magnetic Particle Processor using the MagMAX CORE Nucleic Acid Purification Kit (ThermoFisher Scientific, Waltham, MA, USA). Samples were tested with pan-flavivirus RT-PCR modified from Chao et al. [41] as described in detail by Becker et al. [42]. WNV-positive mosquito pools were subjected to Sanger sequencing (LGC Genomics, Berlin, Germany) for complete genome sequencing [43].
The blood-fed specimens were individually placed into 2 mL safe-lock tubes. Homogenization and extraction were conducted using the same protocol as described above. Thereby, 30 µL supernatant from each of ten specimens was pooled for WNV screening. Detection of WNV-RNA was conducted with the RealStar WNV RT-PCR Kit 1.0 (altona Diagnostics, Hamburg, Germany).
For the host identification, the supernatant of individual blood-fed specimens was heat-inactivated at 99 • C for 1 min in a Peqlab thermocycler (VWR International GmbH, Darmstadt, Germany) for the reduction of possible inhibitors. The PCR assay used the Phusion Blood Direct Master Mix (Thermo Fisher Scientific, MA, USA), 5 µL of the homogenate was used in a total of 30 µL reaction volume for PCR amplification of the cytochrome b gene [44,45]. Amplification was conducted by incubation for 5 min at 98 • C, followed by 40 cycles of 1 s at 98 • C, 5 s at 57 • C and 30 s at 72 • C, ending with incubation for 1 min at 72 • C. If the reaction with the first primer set yielded no result, the PCR reaction was repeated using another pair of vertebrate-specific primers targeting the 16S rRNA gene fragment [46]. The same applied to potential mixed blood meals, as indicated by double peaks in the sequence electropherograms at different positions, resulting in unreadable chromatograms. For this PCR, amplification was conducted by incubation for 5 min at 98 • C, followed by 40 cycles of 1 s at 98 • C, 5 s at 50 • C and 30 s at 72 • C, concluded by incubation for 1 min at 72 • C. The amplicons were sequenced (LGC Genomics, Berlin, Germany) and analyzed with Geneious v9.1.7 (Biomatters, Auckland, New Zealand). Sequences were compared to available sequences from GenBank database (https://blast.ncbi.nlm.nih.gov/). Host species were identified if the percentage identity was 95% or higher. The statistical computer program R [47] was used for all data analyses. Data manipulation and visualization was conducted with functions from the packages plyr [48], dplyr [49], magrrittr [50] and ggplot2 [51]. Spearman's rank correlation was used to analyze the statistical relationship between the number of analyzed specimens per mosquito species and the number of detected host species. For each mosquito species, higher order taxa (e.g., Anatidae, Bovidae, Chiroptera) were only considered for the calculations of host species, if no corresponding taxa of lower ranks were detected. The frequencies of detected birds, non-human mammals or humans between the six most abundant mosquito species and between the four sampling sites were compared with Chi-square tests with Bonferroni corrected p-values for multiple pairwise comparisons.
Horse-, human-, dog-and bird-derived blood meals were tested for WNV-specific IgG/IgY, using an indirect immunofluorescence (IIF) assay as described previously [52]. Host species were selected, which are important amplifying hosts (bird), known to become critically ill from WNV infections (human, horse) or were previously identified to be suitable sentinel species for WNV (dog, horse) [6,7,53,54]. In brief, Vero cells infected with WNV NY99 were seeded on microscope slides with 12 reaction wells (Marienfeld, Lauda-Königshofen, Germany). Slides were treated with Acetone (99%), 15 µL of each sample (single mosquito homogenized in 500 µL) was transferred into one reaction well. Cells were washed with PBS and stained with Alexa Fluor ® 488-conjugated Alpaca Secondary Antibodies (Jackson ImmunoResearch, West Grove PA, USA 1:200 in 1% Evans blue solution), namely goat anti-human IgG, goat anti-horse IgG, rabbit anti-chicken IgY and rabbit anti-dog IgG antibodies, depending on the identified blood-meal source. In order to test for cross-reactivity with heterologous flaviviruses potentially circulating in the sampling area, the WNV IgG positive samples were also tested for Usutu virus-(USUV) and tick-borne encephalitis virus-(TBEV) specific IgG using the same assay with the respective virus.
Genomes obtained for WNV strains from Danube Delta were compared with all complete and partial publicly available NS5 gene sequences from Europe and Africa. Phylogenetic trees were inferred using the Bayesian Markov chain Monte Carlo (MCMC) approach available in BEAST v1.10 [55]. Analyses were performed under the best fit nucleotide substitution model identified as the GTR +Γ for complete genome and TN93+Γ for partial NS5 datasets using jModelTest 2 [56] and a prior MCMC was chosen by testing all models and determining Bayes factors (log 10 BF). We employed TempEst for an interactive regression approach to explore the association between genetic divergence through time and sampling dates [57]. In order to assess the spatial temporal dynamics of WNV, the time to most recent common ancestor (tMRCA), and the effective population dynamics of WNV, we employed a relaxed uncorrelated log normal (UCLN) molecular clock, a flexible demographic model (coalescent Gaussian Markov Random field Bayesian Skyride model, GMRF) as the best demographic scenario detected. In all cases, each of the MCMC chain lengths was run for 10 8 generations (with 10% burn-in) and subsampled every 10 4 iterations to achieve convergence. The Bayesian maximum clade credibility (MCC) trees were visualized using FigTree v1.4.1 (http://tree.bio.ed.ac.uk/software/figtree/). To test the hypothesis that WNV is periodically imported from Africa into Europe, a phylogeographic analysis was conducted using a discrete model attributing state characters represented by the detection locality of each strain and the Bayesian stochastic search variable (BSSV) algorithm implemented in BEAST v1.10 [55].

Genome Characterization of WNV in the DDBR
Both WNV positive mosquito pools have been subjected to Sanger sequencing for complete genome sequencing as described elsewhere [43] and deposited in GenBank under the accession numbers MH939153 and MH939154. Sequence comparison between the two sequenced genomes revealed 51 nucleotide (identity rate 99.5%) and 8 amino acid (identity rate 99.8%) differences almost all of them distributed along the polyprotein ( Figure S1). Several structural and nonstructural genes of the WNV from Danube Delta exhibited unique or similar amino acid changes exclusively with African WNV strains ( Figure S2).

Phylogeography and Spatio-Temporal Dispersal Pattern of WNV
The Bayesian phylogenetic tree of the complete coding sequence of WNV showed that the strains from Danube Delta clustered in the Eastern European clade 1 together with the strains Hyalomma/Romania/2013, Volgograd/2007 and Italy 792/14 ( Figure 1). Given that the majority of available sequences from Romania are partial NS5 gene fragments, we have inferred Bayesian MCC phylogenies with similar topologies as for the complete genome-based tree. In addition, the phylogenies revealed that the Romanian WNV strains fell into two distinct monophyletic clades within WNV phylogeny, suggesting two distinct introductions into Romania (Figure 2b). One clade designated as Eastern European clade 1 (EEC1) included all WNV strains from Danube Delta and some from south-east Romania (Bucharest), while the second clade designated as Western European clade 1 (WEC1), which forms also a distinct monophyletic clade with WNV strains from south-east Romania, but not from DDBR ( Figure 1b, Figure 2b). To assess the viral migration and origin of the WNV in Romania, a discrete-trait phylogeography analysis [58] using the complete genome and NS5 datasets was used to reconstruct the WNV movements between continents/countries. Both datasets exhibited a strong temporal signal and the coefficient of rate variation supported the use of a relaxed clock model ( Figure 1a, Figure 2a). The phylogenetic analysis revealed that the long-distance movement pattern of WNV between Africa and Europe occurred. We estimated at least 6 intercontinental and 10 continental (Europe) viral migration events ( Figure 3). For Romania, we observed at least 2 distinct introduction events ( Figure 1b, Figure 2b). The limited number of available sequences from Romania and the lack of WNV data from several European and African countries make it difficult to infer with confidence the spatiotemporal pattern of WNV EEC1. The time to the most recent common ancestor (tMRCA) of the Romanian WNV strains from WEC1 and EEC1 clades indicates a very recent emergence. These strains were most likely introduced into Romania during the 1990s and 2000s as two distinct introduction events (Figures 1-3). The EEC1 clade seems to be a descendant of an ancestor that probably emerged in South Africa around 1910 (95% HPD for 1901-1920; posterior probability 0.96) (Figure 1b), while the WEC1 clade shares a common ancestor that probably emerged in Greece around 1999 (95% HPD for 1994-2003; posterior probability 0.99) (Figure 2b). The spatial origin and diffusion patterns of the WNV were reconstructed using a BSSV analysis. The earliest introduction and migration event of WNV lineage 2 in Romania (Danube Delta) was detected from South Africa (analysis based on complete genome) or Senegal (based on NS5) between 1992 and 2001, after which the virus dispersed to Russia, Italy and southeast Romania (Bucharest) (Figures 1-3). The second origin and introduction of WNV in Romania was detected to be from Greece (based on NS5) between 2001 and 2002 (Figures 1-3). Furthermore, the phylogeographic analysis also revealed the co-circulation of both EEC1 and WEC1 in south-east Romania, but not in Danube Delta (Figures 1-3).     The analysis of the complete genome sequences revealed nonsynonymous geographic and clade-specific mutations in all members of EEC1 including some African ancestral specific amino acid residues which further strengthen the African origin of the WNV circulating in Danube Delta ( Figure S2).

Screening for WNV-Specific IgG Antibodies
Nineteen blood meals (2.2%, n = 858 analyzed mosquito specimens) contained WNV-specific antibodies ( Table 1, Table S2). Seven of these samples originated from dogs (6.3%, n = 111) and 12 from horses (3.1%, n = 391). All blood meals from birds (n = 85) and humans (n = 271) were WNV IgY/IgG negative. Positive samples were detected for all four sampling sites. WNV IgG positive samples were also tested for USUV-and TBEV-specific IgG. Only one WNV IgG positive blood meal from a dog was also tested positive for USUV-specific IgG.
The success rate for the identification of the blood sources was 60.7% (2331 specimens), amounting to 2348 identified hosts from 43 species and three unspecified taxa (Anatidae, Bovidae, Chiroptera). The difference of 17 specimens between the number of detected hosts and analyzed mosquito specimens results from mixed blood-meals (Table 2). Hosts were detected for 12 (92.3%) out of the 13 analyzed blood-fed mosquito species, with no successful PCR amplification for a single specimen of Culex martinii. The largest number of host taxa was detected for Cq. richiardii (30 species, 827 successfully analysed specimens), followed by Cx. pipiens s.l./torrentium (17 species, 56 specimens), An. maculipennis s.l. (15 species, 280 specimens) and An. hyrcanus (13 species, 791 specimens). Both, Ae. vexans and Ae. caspius fed on a moderate number of host species (nine host species each, 230 and 106 specimens) (Table S3). The other mosquito species with two to 18 specimens fed on two to five host species. Not surprisingly, the number of collected specimens and detected host species were statistically significantly correlated (Spearman rho = 0.90, p < 0.001).
These three most common host groups were determined for all six most abundant mosquito species (Figure 4). The ratios of the three host groups were statistically significantly (χ2 = 252.72, df = 10, p < 0.001) different between the six most abundant mosquito species. All five taxa except Cx. pipiens s.l./torrentium showed similar host-feeding patterns with clear preference for non-human mammals (81.2-95.7% of detected blood meal sources), followed by the host groups human and bird with 3.5-16.7% and 0.9-3.9%, respectively. However, Ae. vexans had statistically significantly different host group proportions with a higher proportion of non-human mammals compared to An. hyrcanus, An. maculipennis s.l. and Ae. caspius (adjusted p values < 0.05, Table S4). Similarly, we observed a higher proportion of non-human mammals for An. maculipennis s.l. and An. hyrcanus compared to Cq. richiardii (adjusted p values < 0.01, Table S4). Furthermore, the host-feeding pattern of Cx. pipiens s.l./torrentium was found to be statistically different (adjusted p values < 0.001, Supplementary Table  S4) compared to all other five abundant species, with 53.2% non-human mammals, followed by 35.5% birds and 11.3% human.  Table S4). Similarly, we observed a higher proportion of non-human mammals for An. maculipennis s.l. and An. hyrcanus compared to Cq. richiardii (adjusted p values < 0.01, Table S4). Furthermore, the host-feeding pattern of Cx. pipiens s.l./torrentium was found to be statistically different (adjusted p values < 0.001, Supplementary Table S4) compared to all other five abundant species, with 53.2% non-human mammals, followed by 35.5% birds and 11.3% human. The sampling sites had statistically significant different compositions of detected host groups ( Figure 4, χ2 = 114.41, df = 6, p < 0.001). There were no statistical differences between the two nearnatural sites Dunărea Veche and Lake Roșuleț (adjusted p values > 0.05, Table S5), but all other combinations were statically significantly different (adjusted p values < 0.05, Table S5, Figure S3). We observed higher percentages of non-human mammals and lower percentages of humans for the two more anthropogenic influenced sites, Letea and Sulina. Furthermore, we detected higher relative proportions of birds and humans for Sulina compared to Letea.

Discussion
In this study, we elucidated the possible origin, pattern of spatial-temporal dynamics, and ecoepidemiological factors of WNV in the ecosystem of DDBR. Our phylogeographic analysis identified at least two distinct introduction events of WNV lineage 2 to Romania. It circulates under a number of different virus variants (EEC1 and WEC1) with South Africa/Senegal and Greece as a possible hub for the progenitor of WNV strains involved in the outbreaks in Romania. The presence of a The sampling sites had statistically significant different compositions of detected host groups (Figure 4, χ2 = 114.41, df = 6, p < 0.001). There were no statistical differences between the two near-natural sites Dunărea Veche and Lake Ros , ulet , (adjusted p values > 0.05, Table S5), but all other combinations were statically significantly different (adjusted p values < 0.05, Table S5, Figure S3). We observed higher percentages of non-human mammals and lower percentages of humans for the two more anthropogenic influenced sites, Letea and Sulina. Furthermore, we detected higher relative proportions of birds and humans for Sulina compared to Letea.

Discussion
In this study, we elucidated the possible origin, pattern of spatial-temporal dynamics, and eco-epidemiological factors of WNV in the ecosystem of DDBR. Our phylogeographic analysis identified at least two distinct introduction events of WNV lineage 2 to Romania. It circulates under a number of different virus variants (EEC1 and WEC1) with South Africa/Senegal and Greece as a possible hub for the progenitor of WNV strains involved in the outbreaks in Romania. The presence of a geographically distinct WNV clade (WEC1) is likely due to very recent introduction, adaptation to the local ecological conditions and some geographic barriers such as climate, vegetation, and vector species. Furthermore, the long-term circulation (EEC1) and adaptation of the virus to the host populations and its enzootic maintenance lead to spread into new geographic regions and local virus variants (in situ evolution).
Although the overlap between the phylogenetic and geographical clustering of the Romanian and Russian members of the Eastern European clade of WNV lineage 2 was expected, it is interesting to note that the clade also contains an Italian strain. This suggests a new, independent introduction of the EEC1 in the south-central part of the continent. Similarly to Eastern Europe, the Italian Peninsula is crossed by major Afro-European bird migration routes. To date, the dispersion pattern of WNV into temperate Eurasia can be best explained by bird migration [59][60][61][62], with short-distance migratory species as potential mode of WNV spread within Europe [62]. Interestingly, we found evidence of adaptive evolution in the WNV from Danube Delta also in non-structural genes, which likely indicates that the host immune selection pressure does not cause increases in viral fitness [63]. Mutations observed at amino acid positions T108I (C), S196P and R361K (E), I1192V (NS2A) and G2932R (NS5) have been found to be involved in the formation of EEC1. Although the impact of these mutations mostly from the nonstructural genes is unclear (likely occurred due to introduction of WNV in this country), similar changes modulated the host antiviral response by inhibition of interferon signaling [64]. The residue alternations R851K (NS1), I1462M (NS2B), R1516K (NS3), T2296A (NS4), N2305S (NS4) and R2719K (NS5) are specific for African variants. Similar patterns of convergent evolution have been described for WNV and suggest that a limited number of residue changes are permitted due to functional constraints [65].
This study successfully used a xenosurveillance approach to monitor the presence of WNV-specific antibodies in different host species. As demonstrated previously [66], mosquito-based surveillance allows non-invasive blood-sampling from free-roaming vertebrate hosts (e.g., feral horse) and from species which are rare or have a cryptic behavior (e.g., raccoon dog (Nyctereutes procyonoides), European mink (Mustela lutreola), Eurasian otter (Lutra lutra) or Golden jackal (Canis aureus)). This study supports previous studies, which identified horses and dogs as suitable sentinel species for WNV [53,54]. WNV seroprevalence in dogs (6.3%) and horses (3.1%) was similar to conventional sampling of the species in different areas with WNV activity [67][68][69]. Nevertheless, the seroprevalence in horses was markedly lower than in southeastern Romania (15.1%) [70]. WNV-specific antibodies were found in blood meals from horses and dogs in all four sampling sites, but not in mosquito blood-meals from human or bird. Thus, this indicate widespread, continuous WNV circulation, but probably only on a low level. In addition, due to potential cross-reactivity of the applied serological assay, we cannot exclude the possibility that one of the samples was also positive for USUV, a virus with a similar transmission cycle to WNV.
Feral horses and free-ranging cattle were the two most commonly detected host species, accounting for more than 50% of the detected hosts. These animals have their origins in pre-1990 state-owned collective farms and private homesteads from where they were released in more recent years. There is no official census published, but it is estimated that 4000 horses and a few other thousand cattle roam and reproduce freely in the DDBR. The high abundance in combination with the relatively huge body size [71] might explain why both host species are facilitated so often. However, there were differences between the sampling sites. Non-human mammals dominated the detected hosts for the two sampling sites located in the interface between anthropic and natural landscapes (Sulina and Letea), i.e., homesteads in direct proximity of livestock. In contrast, the other two near-natural sampling sites (Lake Ros , ulet , and Dunărea Veche) were both located deep inside the Danube Delta and only insignificantly anthropogenically influenced. However, humans are commonly present in a fishing cabin and an agricultural holding. In the absence of high abundances of cattle and horse, mosquitoes might rather select other available hosts, e.g., birds and humans. Thus, host availability probably is a decisive factor for the host-selection of mosquitoes, influencing the risk of local pathogen transmission.
At the same time, this study highlights the importance of Cx. pipiens s.l./torrentium as a WNV vector in Europe [72]. The only two WNV positive pools belonged to this taxon. Previous studies described the species complex as predominantly ornithophilic [30][31][32]60,73]. However, this study again demonstrates its catholic host-feeding pattern. As the other five most abundance taxa, Cx. pipiens s.l./torrentium predominantly fed on non-human mammals and humans but had the highest proportion of birds, i.e., a more than nine times higher proportion, making the species a potential bridge vector. This is in line with studies from Africa [74], Middle East [75], Europe [37], and North America [76]. Although several other collected mosquito species (e.g., An. hyrcanus) have been found positive for WNV-RNA in Romania [11,13,24], members of the Culex pipiens complex are considered the main vectors for WNV in both urban and rural/natural transmission cycles [13].
Active WNV circulation in the DDBR is strongly implicated by WNV infections of unfed mosquito specimens and serological evidence of WNV-specific antibodies in the hosts. However, this study collected only few specimens of the most competent vector Cx. pipiens s.l./torrentium. In previous studies in the DDBR [24,77], >95% of mosquitoes captured with bird-baited traps belonged to this species complex, while the most abundant species in our study (Cq. richiardii) was absent. The usage of a single type of trap in only a few sampling sites across the delta's heterogenous landscape is likely to have contributed to a biased sampling outcome [78,79]. In addition, mobovirus transmission is generally restricted to small foci in non-epidemic years [80]. Therefore, the presence of WNV can be underestimated depending on the number, type and location of traps. Further studies are needed to identify and further understand the driving factors of landscape and time.

Conclusions
The detection of WNV-RNA and WNV-specific antibodies confirms the circulation of this important mobovirus in the DDBR. Serological evidence for WNV circulation confirms the applicability of mosquito-based surveillance in sero-epidemiological studies. Host identification for blood-fed mosquitoes allows the usage of host-specific conjugates. In addition, host-feeding patterns of Cx. pipiens s.l./torrentium underly the relevance of the taxon as an enzootic and bridge vector for WNV in Europe, which was further confirmed by the detection of WNV lineage 2 RNA in two pools of unfed specimens from the same taxon. Local overwintering or reintroduction of the virus could be considered decisive factors for the evolution, dispersal and endemisation of WNV in temperate Europe. Thus, to better understand the impact of ecological/immunological factors on WNV evolution, studies based on more comprehensive genetic data, including those from previously unsampled geographic regions, are required.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4915/11/12/1159/s1, Figure S1: Schematic representation of the nucleotide (a) and amino acid (b) differences between WNV strains from Danube Delta generated during this study, Figure S2: Schematic representation of the specific amino acid replacements in the Danube Delta WNV genomes, Figure S3: Percentage of host-feeding groups (birds, human, non-human mammals) for the four sampling sites in the DDBR, Romania (2014-2016), Table S1: Frequency and percentage of mosquitoes collected at four sampling sites in the DDBR, Romania (2014-2016), Table S2: Samples of blood-fed mosquito species positive for West Nile virus-specific IgG antibodies with information on the host species, mosquito species, sampling year collected at four sampling sites in the DDBR, Romania (2014-2016), Table  S3: Number and percentage (in brackets) of detected host taxa for the six less abundant species, Table S4: Results of Chi-square tests for the comparison of host-feeding groups between the six most abundant mosquito species with adjusted p-values collected at four sampling sites in the DDBR, Romania (2014-2016), Table S5: Results of Chi-square tests for the comparison of host-feeding groups of mosquitoes between the four sampling sites in the DDBR, Romania