First Detection and Molecular Characterization of Usutu Virus in Culex pipiens Mosquitoes Collected in Romania

Usutu virus (USUV) is an emergent arbovirus in Europe causing mortality in bird populations. Similar to West Nile virus (WNV), USUV is maintained in sylvatic cycles between mosquito vectors and bird reservoirs. Spillover events may result in human neurological infection cases. Apart from indirect evidence provided by a recent serological study in wild birds, the circulation of USUV in Romania was not assessed. We aimed to detect and molecular characterize USUV circulating in mosquito vectors collected in South-Eastern Romania—a well-known WNV endemic region—during four transmission seasons. Mosquitoes were collected from Bucharest metropolitan area and Danube Delta, pooled, and screened by real-time RT-PCR for USUV. Partial genomic sequences were obtained and used for phylogeny. USUV was detected in Culex pipiens s.l. female mosquitoes collected in Bucharest, in 2019. The virus belonged to Europe 2 lineage, sub-lineage EU2-A. Phylogenetic analysis revealed high similarity with isolates infecting mosquito vectors, birds, and humans in Europe starting with 2009, all sharing common origin in Northern Italy. To our knowledge, this is the first study characterizing a strain of USUV circulating in Romania.


Introduction
Usutu virus (USUV) is a single-stranded, positive-sense RNA virus, belonging to Flavivirus genus, Flaviviridae family. Similar to other flaviviruses, this emerging pathogen is maintained in sylvatic cycles between ornithophilic mosquitoes as vectors and birds as amplifying hosts. Mammals, including humans, are considered "dead-end" hosts, as they can be incidentally infected through mosquito biting, but do not develop a sufficient viremia to sustain transmission [1,2].
USUV was first isolated in 1959 from Culex neavei mosquitoes collected in South Africa [3]. In Europe, USUV was first detected in Austria, in 2001, when it produced increased mortality in wild birds in Vienna [4]. However, retrospective study based on archival bird tissues presented evidence of a much earlier circulation of the same strain in 1996, in Tuscany, Italy [5]. Since the first report, the virus was confirmed as cause of mortality in wild and captive birds in Europe, especially in Passeriformes (e.g., common blackbirds-Turdus merula) and Strigiformes populations [6][7][8][9]. As reviewed [2], virus isolation, nucleic acid detection, and seroconversion studies detected USUV in a wide range of other vertebrates in Europe, such as horses, dogs, bats, squirrels, wild boars, roe deer, bats, and lizards. Similar to the West Nile virus (WNV) transmission cycle, Culex pipiens mosquitoes act as vectors of USUV in Europe [10], although invasive species such Aedes albopictus [11] and Aedes japonicus [12] were also found infected.
Worldwide diversity of USUV is reflected in the existence of eight genetic lineages (Africa 1-3 and Europe 1-5) [13,14]. Despite being recently detected in Europe, phylogeny study concluded that USUV was introduced regularly from Africa into Europe in the last 50 years [13]. Co-circulation of both African and European lineages was documented in Europe or sometimes even in the same country. However, Europe 2 lineage is considered to be the most prevalent there [2].
To date, only approximately one hundred symptomatic and asymptomatic USUV human infection cases were recorded globally, mostly in Europe. Thirty of those cases were neuroinvasive infections [15].
Seroprevalence study on wild birds conducted in South-Eastern Romania between 2018 and 2019 presented indirect evidence of USUV circulation in this country [16]. However, serological survey conducted between 2019 and 2020 on healthy blood donors showed no evidence of USUV circulation in human population from the North-Western region of Romania [17].
Knowledge on USUV circulation in Romania is extremely scarce and there is no sequence data regarding the circulating strains. To overcome these gaps, we set up an entomological study aiming to characterize the USUV detected in pools of mosquitoes collected in Bucharest and Danube Delta, during four transmission seasons.

Sampling Sites and Mosquito Collection
Adult mosquitoes were collected each year between July and September. Nine sites located in Bucharest metropolitan area and four sites located in Danube Delta (the left bank of the Sulina branch, opposite to Vulturu village, Tulcea county) were investigated in 2012 and 2013. During 2019 and 2020, the study was carried out in thirteen sites located in Bucharest metropolitan area ( Figure 1). Collections were performed using Centers for Disease Control and Prevention (CDC) light and gravid traps (John W. Hock Company, Gainesville, FL, USA), BG-Sentinel traps (Biogents, Regensburg, Germany), bird-baited traps, and mammal-baited traps. Chickens and guinea pigs were used as baits after obtaining the approval of Cantacuzino Institute Ethics Committee. Animal-baited traps were used only in 2012-2013 seasons in the Danube Delta. CDC light traps and BG-Sentinel traps were used in 2012-2013 seasons both in Danube Delta and Bucharest metropolitan area. In 2019-2020 seasons we exclusively used CDC gravid traps. All these traps proved good collection rates and can operate continuously for several days. CDC gravid traps attract mainly pregnant mosquito females, thus increasing the chances of finding arbovirus positive samples. Mosquitoes were transported to the laboratory in liquid nitrogen containers and there were kept at −70 • C freezer until further processing.

Mosquito Processing and RNA Extraction
Mosquitoes were identified to species level on chill table under stereomicroscope, using dichotomous keys described by Becker et al. [18] and grouped according to species, sex, physiological age, site and date of collection, into pools never exceeding 50 individuals. Viral RNA was extracted from pool homogenates using QIAamp Viral RNA Mini Kit (Qiagen, Hilden, Germany), as previously described [19].

Mosquito Processing and RNA Extraction
Mosquitoes were identified to species level on chill table under stereomicroscope, using dichotomous keys described by Becker et al. [18] and grouped according to species, sex, physiological age, site and date of collection, into pools never exceeding 50 individuals. Viral RNA was extracted from pool homogenates using QIAamp Viral RNA Mini Kit (Qiagen, Hilden, Germany), as previously described [19].

Detection of USUV by Real-Time RT-PCR
Mosquito samples collected between 2012 and 2013 were tested using primers and probe described by Del Amo et al. [20] and GoTaq ® Probe 1-Step RT-qPCR System (Promega, Madison, WI, USA). Briefly, reaction mix contained five microliters of RNA, 0.25 µM of each primer, 0.2 µM probe, 0.4 µL GoScript™ RT Mix for 1-Step RT-qPCR, ten µL of GoTaq ® Probe qPCR Master Mix with dUTP, and RNAse-free water up to 20 µL. Thermal amplification was carried out following the profile recommended by the real-time RT-PCR kit manufacturer on an Mx3005P Real-Time PCR System (La Jolla, CA, USA). For testing the samples collected between 2019 and 2020 a real-time RT-PCR assay developed by Nikolay et al. [21] and SuperScript™ III Platinum™ One-Step qRT-PCR Kit (Invitrogen, Waltham, MA, USA) were used. Primers and probe concentrations were 0.5 µM and 0.2 µM, respectively. A total of 10 microliters of RNA extract were tested in a final reaction volume of 25 µL. Thermal amplification was carried out following the profile recommended by the real-time RT-PCR kit manufacturer on an Mx3005P Real-Time PCR System (La Jolla, CA, USA). Each PCR run was validated by including a positive control (USUV strain SAAR-1776, 1959, GenBank acc. no. AY453412) and a negative control (PCR grade water).

Sequencing and Phylogenetic Analysis
A primer walking strategy [22] was employed for generating amplicons spanning genomic regions preM-M-E and NS4B-NS5-3′UTR, respectively. PCRs were carried out using five microliters of RNA extract, 0.8 µM of each primer, and OneStep RT-PCR Kit

Detection of USUV by Real-Time RT-PCR
Mosquito samples collected between 2012 and 2013 were tested using primers and probe described by Del Amo et al. [20] and GoTaq ® Probe 1-Step RT-qPCR System (Promega, Madison, WI, USA). Briefly, reaction mix contained five microliters of RNA, 0.25 µM of each primer, 0.2 µM probe, 0.4 µL GoScript™ RT Mix for 1-Step RT-qPCR, ten µL of GoTaq ® Probe qPCR Master Mix with dUTP, and RNAse-free water up to 20 µL. Thermal amplification was carried out following the profile recommended by the real-time RT-PCR kit manufacturer on an Mx3005P Real-Time PCR System (La Jolla, CA, USA). For testing the samples collected between 2019 and 2020 a real-time RT-PCR assay developed by Nikolay et al. [21] and SuperScript™ III Platinum™ One-Step qRT-PCR Kit (Invitrogen, Waltham, MA, USA) were used. Primers and probe concentrations were 0.5 µM and 0.2 µM, respectively. A total of 10 microliters of RNA extract were tested in a final reaction volume of 25 µL. Thermal amplification was carried out following the profile recommended by the real-time RT-PCR kit manufacturer on an Mx3005P Real-Time PCR System (La Jolla, CA, USA). Each PCR run was validated by including a positive control (USUV strain SAAR-1776, 1959, GenBank acc. no. AY453412) and a negative control (PCR grade water).

Sequencing and Phylogenetic Analysis
A primer walking strategy [22] was employed for generating amplicons spanning genomic regions preM-M-E and NS4B-NS5-3 UTR, respectively. PCRs were carried out using five microliters of RNA extract, 0.8 µM of each primer, and OneStep RT-PCR Kit (Qiagen, Germany) according to the kit insert. Amplicons were sequenced with BigDye™ Terminator v3.1 Cycle Sequencing Kit on a SeqStudio™ Genetic Analyzer (Applied Biosystems, Waltham, MA, USA). Consensus sequences were obtained with BioEdit version 7.2.5 [23]. Phylogenetic trees were generated using MEGA11 version 11.0.13 [24] after choosing the fittest nucleotide substitution model with the same software. The tree annotations were added using EvolView v2 [25].

Results
A total number of 13,823 mosquitoes were collected during the study. Urban collections comprised specimens belonging to C. pipiens s.l. and invasive A. albopictus species, whereas Danube Delta collections contained specimens belonging to four species (Table 1).  Figure 1).
First partial genomic sequence (1857 nt) contained almost the entire coding sequence of preM, whole coding sequence of M protein, and partial E protein (nucleotide positions 591-2447, reference sequence GenBank acc. no. NC_006551). The second genomic fragment (3154 nt) spanned positions 7356-10,509 and contained partial coding sequence of NS4B protein, entire coding sequence of NS5 viral polymerase, and partial 3 UTR. The two sequences were deposited in GenBank as a gapped format submission under the acc. no. OQ414983. The two sequences were used for assigning genetic lineage of the detected USUV [13,14,22,26]. Phylogenetic analyses placed the USUV detected in this study into Europe 2 lineage, group A/sub-lineage EU2-A. Romanian sequences clustered with USUV sequences obtained from Culex sp. mosquito vectors, birds (common blackbird and great tit), and human samples collected starting with 2009 in Italy, Central Europe (Austria, Czech Republic, and Hungary), and Germany. The topology of the two phylogenetic trees constructed was consistent (Figures 2 and 3).
We investigated our sequences for the presence of phylogenetic informative substitutions (i.e., cluster or country specific substitutions) [13,26], host specific mutations [13], and markers linked to neuroinvasiveness [27]. Substitutions G595S and E3425D specific for EU2-A sub-lineage and associated with neuroinvasiveness were identified in Romanian sequences. M2645I substitution characteristic to Italian isolates was also detected.  (EU1-EU4) showing the Romanian sequence clustering into sub-lineage EU2-A. White dot: sequence obtained in this study. Lineages were assigned according to [13,14,22,26] and similar sequences retrieved by NCBI BLAST were also included in the analysis. The evolutionary history was inferred from a preM-M-E fragment (nucleotides 591-2447, reference sequence GenBank acc. no. NC_006551) by using the Maximum Likelihood method and Kimura 2-parameter model, 1000 bootstrap replicates. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 0.1000)). Numbers at nodes represent the bootstrap percentages (values <70% are not shown).  (EU1-EU4) showing the Romanian sequence clustering into sub-lineage EU2-A. White dot: sequence obtained in this study. Lineages were assigned according to [13,14,22,26] and similar sequences retrieved by NCBI BLAST were also included in the analysis. The evolutionary history was inferred from a preM-M-E fragment (nucleotides 591-2447, reference sequence GenBank acc. no. NC_006551) by using the Maximum Likelihood method and Kimura 2-parameter model, 1000 bootstrap replicates. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 0.1000)). Numbers at nodes represent the bootstrap percentages (values <70% are not shown).  (EU1-EU4) showing the Romanian sequence clustering into sub-lineage EU2-A. White dot: sequence obtained in this study. Lineages were assigned according to [13,14,22,26] and similar sequences retrieved by NCBI BLAST were also included in the analysis. The evolutionary history was inferred from a NS4B-NS5 fragment (nucleotides 7372-10,398, reference sequence GenBank acc. no. NC_006551) by using the Maximum Likelihood method and Tamura-Nei model, 1000 bootstrap replicates. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 0.1000)). Numbers at nodes represent the bootstrap percentages (values <70% are not shown).
We investigated our sequences for the presence of phylogenetic informative substitutions (i.e., cluster or country specific substitutions) [13,26], host specific mutations [13], and markers linked to neuroinvasiveness [27]. Substitutions G595S and E3425D specific  (EU1-EU4) showing the Romanian sequence clustering into sub-lineage EU2-A. White dot: sequence obtained in this study. Lineages were assigned according to [13,14,22,26] and similar sequences retrieved by NCBI BLAST were also included in the analysis. The evolutionary history was inferred from a NS4B-NS5 fragment (nucleotides 7372-10,398, reference sequence GenBank acc. no. NC_006551) by using the Maximum Likelihood method and Tamura-Nei model, 1000 bootstrap replicates. A discrete Gamma distribution was used to model evolutionary rate differences among sites (5 categories (+G, parameter = 0.1000)). Numbers at nodes represent the bootstrap percentages (values <70% are not shown).

Discussion
In this study we present the molecular characterization of USUV detected in mosquito vectors collected in South-Eastern Romania. To our knowledge, this is the first sequencebased evidence of USUV circulation in our country. However, seroprevalence study showed the circulation of USUV in wild birds sampled in localities in Tulcea county, including Danube Delta Biosphere Reserve, and in Constanta county [16] (Figure 1), but not in human population in other part of the country [17]. The sampling sites chosen for our four-year study were Danube Delta and Bucharest metropolitan area. Danube Delta is the largest wetland in Europe and a major stopover for wild birds migrating between Europe and Africa. It was suggested that the migratory birds might act as potential long-distance dispersal vehicles of USUV and that the bird flyways could predict the continental and intercontinental dispersal patterns of this virus [13]. Moreover, Danube Delta was shown to be an area of circulation and introduction for WNV which displays common ecology and epidemiology traits with USUV [28][29][30][31]. Bucharest, the second sampling region in our study, is also known for the intense circulation and human transmission of WNV, significant outbreaks occurring there in the last three decades [32][33][34], with C. pipiens mosquitoes acting as vectors [29,35,36].
We detected USUV by real-time RT-PCR in C. pipiens s.l. female pool collected in 2019, in Bucharest. Our phylogenetic analyses showed that the virus belong to Europe 2 lineage (EU2). As reviewed elsewhere [2,37], this lineage is the most prevalent on our continent and circulates or co-circulates with other European and African lineages in countries from Central Europe (Austria, Hungary, Czech Republic, and Slovakia), Western Europe (France and Germany) and Southern Europe (Serbia and Croatia). Large scale studies aiming to reconstruct the evolutionary history of USUV have shown that European lineages share a common ancestor and are the result of local evolution [13] and placed the origin of EU2 between 1993 [13] and 2003 [26]. Italy is considered a hot spot for USUV evolution and a dispersal source of European lineages throughout the continent. Recent study on the Italian USUV isolates demonstrated that the EU2 can be split in two groups/sub-lineages, EU2-A and EU2-B, respectively [26]. Our sequence clusters within sub-lineage EU2-A which is considered to have emerged in February 2007 and contains samples collected between 2009 and 2018 from Northern and Central Italian regions, Austria, and Hungary [26]. Specific mutations have been described for the European lineages and sub-lineages. G595S (E protein) and E3425D (NS5 protein) substitutions are considered EU2-A amino acid signatures [26] and are present in the Romanian sequence, too. Substitution G595S was also considered a hallmark of human isolates [13] and along with E3425D was first detected in the viral isolate infecting the first human patient presenting neurological symptoms. G595S substitution is located in DIII domain of E protein that is recognized by neutralizing antibodies. It was hypothesized that the two substitutions play a role in the neuroinvasive capacity of USUV [27]. The relatedness of the Romanian sequence with the Italian cluster is indicated also by the amino acid substitution M2645I (NS5 protein) considered characteristic for the Italian isolates [26].
It was demonstrated that USUV spread in Europe was from Northern Italy to Austria and from there to Hungary and subsequently to Serbia. The study proposes bird populations migration-namely of Turdus spp.-as mechanism for virus dispersal across the continent [26]. Due to the high similarity between the Romanian and the cluster of sequences originating in Italy, we can assume that a potential spread from neighboring Hungary or Serbia could have been the source of USUV in Romania.
Although we could not find evidence of USUV circulation in mosquito vectors collected between 2012 and 2013 in Danube Delta, a latter seroprevalence study (2018)(2019) documented the presence of USUV in wild birds sampled in the area, mainly belonging to Passeriformes and Bucerotiformes orders. Of note, the species presenting the highest seroprevalence in this study was the common blackbird [16]. However, the same study failed to detect USUV specific antibodies in birds sampled in Bucharest in the same period when the C. pipiens s.l. mosquitoes were found positive in our study, most probably because of the limited number of sera tested.
Only one out of the 419 mosquito pool tested in our study was found positive for USUV in four transmission seasons. USUV might be a poor competitor in this WNV endemic region. The two viruses are antigenically related and the antibodies against WNV may cross-react and prevent USUV infection in reservoir birds previously exposed to WNV, as it was shown in vaccinated mice [38]. As well, it was shown in a laboratory study that WNV is selectively transmitted by C. pipiens mosquitoes after being co-exposed to both viruses via an infectious blood meal [39]. The circulation of WNV in 2019, the year when the USUV positive pool was detected, was documented in the Bucharest metropolitan area by 15 confirmed WNV human infection cases [40].
The current available data, although limited, suggested that the level of USUV circulation in South-Eastern Romania was low and might have been overwhelmed by the very active transmission of WNV, which has been causing significant human outbreaks in this region. However, as shown for WNV, in the eco-epidemiological conditions in Romania the emergence of a mosquito-borne virus may occur within a transmission season [36]. Therefore, mosquito and bird surveillance, blood donor screening, and increased awareness regarding the possible involvement of USUV in neuroinvasive infection cases are permanently needed.

Conclusions
In this study we present the first sequence-based evidence for USUV circulation in Romania. The detected strain belonged to Europe 2 lineage, EU2-A sub-lineage, clustering with viruses originating from Northern Italy, Central Europe, and Germany.