Genomic monitoring to understand the emergence and spread of Usutu virus in the Netherlands, 2016–2018

Usutu virus (USUV) is a mosquito-borne flavivirus circulating in Western Europe that causes die-offs of mainly common blackbirds (Turdus merula). In the Netherlands, USUV was first detected in 2016, when it was identified as the likely cause of an outbreak in birds. In this study, dead blackbirds were collected, screened for the presence of USUV and submitted to Nanopore-based sequencing. Genomic sequences of 112 USUV were obtained and phylogenetic analysis showed that most viruses identified belonged to the USUV Africa 3 lineage, and molecular clock analysis evaluated their most recent common ancestor to 10 to 4 years before first detection of USUV in the Netherlands. USUV Europe 3 lineage, commonly found in Germany, was less frequently detected. This analyses further suggest some extent of circulation of USUV between the Netherlands, Germany and Belgium, as well as likely overwintering of USUV in the Netherlands.

health impact of USUV emergence in Europe at this stage. The circulation of USUV in birds can be linked to an increased human exposure to zoonotic risk. Indeed, the detection of viremic blood donations in the Netherlands co-incided with months of active USUV circulation in birds 22 , human cases described in Italy were from an area with demonstrated concomitant circulation of USUV in mosquitoes and birds [16][17][18][19] , and in Austria an increase in human cases was observed alongside the increased USUV activity in birds 23 .
Various USUV lineages are co-circulating in Europe 4,24 . In Germany, circulation of five different USUV lineages (Europe 2, Europe 3, Europe 5, Africa 2, and Africa 3) has been described 25,26 and in France two lineages (Africa 2 and Africa 3) were detected in mosquitoes from the same region 27,28 . The co-circulation of different USUV lineages is thought to be due to independent introduction events to Europe by long-distance migratory birds from Africa, where different USUV lineages are presumed to be circulating, followed by local amplification and evolution leading to a geographical signal in phylogenetic analyses 29 . The assignment and nomenclature of USUV lineages are not standardized, and it is unclear if the lineages differ in their potential to be transmitted by mosquitoes and to infect or cause disease in different wild bird species and/or humans 29 .
In the Netherlands, USUV was first detected in 2016 when it was identified as the cause of an outbreak among blackbirds and captive great grey owls (Strix nebulosa) 28,30 . The virus circulated also in 2017 and 2018, and each late summer to autumn was associated with an increased die-off in blackbirds. Through a national wildlife disease surveillance programme, dead birds reported by citizens were collected, dissected and tissues were submitted for USUV diagnostics. To date, only a limited number of USUV genomes are available in the public domain, and it is unknown if USUV was already present in the Netherlands before 2016. In addition, it is unknown if there was a single introduction event or several independent introduction events, and -if so -what the geographic origin of the viruses is. Recent advances in third generation sequencing technologies have opened up new opportunities for infectious disease research. Pathogen genomics can be used to resolve crucial questions regarding origin, modes of transmission and ecology of emerging viral disease [31][32][33][34][35] . Therefore, we have sequenced 112 USUV genomes from blackbirds brain tissues on the Nanopore sequencing platform to analyse the emergence and spread of USUV throughout the Netherlands in 2016, 2017 and 2018.

Results
USUV genomic sequencing. Between September 2016 and September 2018, 165 dead blackbirds were screened for the presence of USUV by RT-PCR. This screening resulted in 118 USUV positive birds and genomic sequences were generated with an USUV specific multiplex PCR using the Nanopore sequencing technology. Successful sequencing was defined as obtaining at least 80% of the USUV genome with a read coverage threshold of 100x per amplicon 36 . Near complete or complete viral genomes were recovered from 112 of 118 samples ( Table 1). The median Ct value of the USUV positive brain tissues was 20, and the threshold for successful sequencing directly from brain tissue samples was shown to be around Ct value 32. Above this Ct value, only some amplicons were sequenced with a coverage >100x. An overview of the number of reads generated, the percentage of USUV reads and the proportion of successfully sequenced amplicons is displayed in Supplementary Table 1. phylogenetic analysis of USUV in the netherlands. Phylogenetic analysis revealed the co-circulation of USUV lineages Europe 3 and Africa 3 in the Netherlands (Fig. 1), with Africa 3 lineage viruses most frequently detected. The USUV lineage Africa 3 has been previously detected in Germany but not as major variant (9 of the 108 published whole genomes). USUV from blackbirds in the Netherlands are at this stage considered the most numerous within this lineage. In contrast, while the Europe 3 lineage was most commonly found to be associated with blackbird deaths in Germany between 2011 and 2016, this lineage appears to be less frequently detected in the Netherlands. Viruses belonging to the Europe 3 lineage formed several groups which were interspersed with viruses found in neighbouring countries. The two lineages were shown to co-circulate in the Netherlands during 3 consecutive years. In all three years, Africa 3 lineage viruses were more frequently detected (Table 2).
These differences in phylogenetic signal strongly suggests differences in the emergence of these two lineages in the Netherlands, with continuous circulation of the Africa 3 lineage and repeated introductions of viruses from the Europe 3 lineage. However, systematic sequencing of a representative set of samples from birds from other parts of Europe in the same time period is needed to resolve this question. PhyCLIP 37 , a statistically-principled approach to delineate phylogenetic trees into clusters, was used to resolve the phylogenetic clustering of USUV (Fig. 1). PhyCLIP showed that the Africa 3 lineage can be divided into 3 different clusters: (1) a cluster comprising 3 sequences from Germany from 2016 and 1 sequence from the Netherlands from 2018 (2) a cluster comprising 32 USUV sequences detected in the Netherlands between 2016-2018, as well as 3 USUV sequences detected in Germany in 2014 and 2016, and (3) a cluster comprising 55 sequences detected in the Netherlands in 2016-2018 and 4 sequences detected in Germany and Belgium in 2015-2016, in this cluster 10 sequences could further be delineated into a nested sub-cluster. Delineation of these 3 clusters and the nested sub-cluster of the Africa 3 lineage does not appear to be structured over time (Fig. 2). Furthermore, sequences from the Netherlands from clusters Africa 3.2 and Africa 3.3 have been detected from across the country each year, while cluster USUV Africa 3.1 was only detected once in 2018. PhyCLIP delineated the Europe 3 USUV lineage as one clade with two different sub-clusters: (I) a sub-cluster with sequences detected in Belgium, Germany and the Netherlands in 2011-2017 and (II) a sub-cluster detected in Belgium, Germany and the Netherlands in 2011-2018. Also for the Europe 3 lineage, delineation of the sub-clusters does not appear to be structured over time. Four sequences from blackbirds and mosquitos detected in Italy that were previously classified as a distinct clade, Europe 4 24  www.nature.com/scientificreports www.nature.com/scientificreports/ resolution clusters from the Europe 5 lineage, the first consisting of sequences from mosquitos in Israel from 2004, the second consisting of sequences from mosquitos in Israel from 2015. Several USUV sequences could not be classified into clusters. This suggests that more lineages might be circulating, but that the number of sequences available at the moment is too limited to enable delineation into clusters by PhyCLIP.
An overview of the geographical distribution of the USUV clusters detected in dead blackbirds in the Netherlands is displayed in Fig. 1 and in the Supplementory Movie. The first set of dead blackbirds positive for USUV collected in September 2016, were found in 10 of 12 different provinces, indicating already widespread circulation at the time of detection. The USUV identified from this first set were genetically diverse and each of the four USUV clusters identified in the Nertherlands were already present at that time. The different USUV clusters all have a widespread geographic distribution throughout the Netherlands.

Molecular clock phylogeny of USUV lineages detected in the netherlands. BEAST analysis was
performed to determine the time to most recent common ancestor of the USUV Africa 3 and Europe 3 lineages, both detected in the Netherlands. Separate BEAST analyses were performed for each of the two lineages. The

Discussion
The emergence of USUV in Europe is causing massive die-offs of mainly common blackbirds. Besides being of concern for wildlife conservation, the occurrence of numerous and sustained outbreaks in wildlife should be considered as a serious warning signal from the perspective of surveillance for zoonotic pathogens. The increasing numbers of reports of asymptomatic human USUV infections as well as cases of mild to severe neuroinvasive USUV infections in humans may be due to changes in awareness and surveillance but may also be an effect of increased human exposure to this zoonotic risk. Dense genomic monitoring of viral pathogens supports outbreak investigations by providing insights to paterns of transmission. Since USUV was only recently recognized in the Netherlands, we used genomic sequencing to gain insight in how this emerging arbovirus spread and evolve in a previously naïve population. We describe the genetic characterization of USUV genomes from tissue samples from dead blackbirds in the Netherlands throughout 3 consecutive years. An amplicon-based sequencing approach using Nanopore sequencing was used to monitor the genetic diversity of USUV 36 . This protocol proved to be sensitive enough to sequence the vast majority of tissue samples from dead blackbird surveillance. This study shows that genomic sequencing on the Nanopore platform is a powerfull approach for monitoring and tracing ongoing arbovirus infections in dead wildlife.
The USUV Africa 3 lineage was found to be predominantly circulating in the Netherlands. Furthermore, divergence inside this lineage was shown to justify its division into 3 higher resolution clusters. The absence of clear temporal delineation in the phylogeny of the USUV Africa 3 lineage could indicate that subsequent to the first identification of USUV in the Netherlands in 2016, the outbreaks in 2017 and 2018 were more likely caused by USUV lineages that persisted during winter in the Netherlands, rather than by repeated introduction of USUV strains from areas outside the Netherlands. It remains unclear if and how USUV can persist in the Netherlands during winter time, whether through vertical tranmission in mosquitoes, maintainance in birds or overwintering mosquitoes, or if it remains present in other animal species. Therefore, extending the diversity of species included in USUV surveillance efforts would be a useful consideration. A recent study in Germany also showed an increase in the detection of the USUV Africa 3 lineage in 2017 and 2018 26 , however only partial sequence data is available www.nature.com/scientificreports www.nature.com/scientificreports/ and more genetic information is needed before it can accurately be used in phylogenetic analysis. The results differed for the USUV Europe 3 lineage: this lineage was less frequently detected in the Netherlands, and sequences from the Netherlands were interspersed with viruses from Germany and Belgium. The reason for this is unknown but might suggest that the USUV Europe 3 lineage -unlike the USUV Africa 3 lineage -is not enzootic in the Netherlands but periodically introduced from neighbouring countries.
USUV was detected in the Netherlands through bird mortality in 2016, but the most recent common ancestors of both the Africa 3 lineage and Europe 3 lineage are estimated to be well before the initial detection. The two time windows identified are broad, and more data from other regions and previous years has to be produced to generate more precise estimates. However, this information, taken together with the described diversity of USUV detected in the Netherlands since the beginning of the outbreak in September 2016, suggests limited circulation of different USUV lineages in the Netherlands before its initial detection. USUV circulation may have been boosted in 2016 by favourable environmental conditions for mosquitoes 28,30 . Alternatively, several introductions from neighbouring countries could have occurred close to the date of first detection and spread efficiently with these favourable conditions. Phylogenetic and molecular clock analyses conducted for the reconstrution of the Zika virus epidemic in the Americas have estimated the time of introduction of Zika virus more than a year before it was detected through human-disease surveillance 35 . It is unclear whether blackbird mortality associated with USUV infection did not occur prior to 2016 or whether it occurred at levels that did not prompt detection through citizen science-based alerting system. In October 2012, 66 songbirds, including 34 blackbirds, were found dead throughout the Netherlands but all tested negative for USUV 38 . Serological testing of a small number of serum samples from birds in the Netherlands in 2015 did show a seroprevalence of 2.8% for USUV; however, the positive birds consisted of species that are partially migratory and therefore this result does not necessarly indicate local circulation of the virus at that time 39 . Retrospective testing of other samples for the presence of USUV may help to elucidate this. In addition, care has to be taken with the interpretation of the BEAST analysis since recently it was shown that by sequencing an ancient hepatitis B virus the evolutionary rate of the virus has shown to be different than previously thought 40 . To date, only very recent data on USUV is available and we do not know how the virus behaves over a prolonged period of time.
Several USUV sequences remained unclassified by the statistically-principled approach for phylogenetic clustering (PhyCLIP 37 ). This indicates that more lineages might be present, but that the amount of genetic data available is at the moment too limited to enable accurate phylogenetic clustering. In addition, by adding a large amount of sequences to the public database, the genomic information currently available is very skewed towards the USUV strains circulating in the Netherlands, Germany and Italy. Recent identification and generation of USUV genomic data outside Europe (Senegal 5 and Israel 41 ) have had an important impact on the understanding of the global USUV diversity. To get a better overview of USUV diversity in Europe and a better understanding of USUV evolution, it is important that more sequence data is generated from other countries. Given the detection of USUV in the Netherlands, it is important to identify the vector competence in local mosquitos, although a previous studie has already shown the vector competence in the Culex pipiens mosquito which is present in the Netherlands 3 .
This study describes a proactive investigation of USUV emergence and spread in the Netherlands through a national wildlife disease scanning surveillance programme associated with the application of a protocol enabling fast and accurate complete genome sequencing of USUV on the Nanopore platform. We add knowledge relating to the USUV epidemiology and describe the genomic diversity profile of USUV circulating in the country. Our analysis suggests that USUV is likely circulating between neighboring countries in Western Europe, where it has been established and is overwintering. We report an important genomic diversity of the viruses circulating in the Netherlands, observable already at the time of the first identifications of USUV in the country. We highlight issues of bias in surveillance, as well as the possibility of eventual silent circulation of the virus in European countries preceding epizootics detection. Another virus circulating in Europe, West-Nile virus, has antigenic cross-reactivity, similar transmission cycle and may interact at population level with USUV 42 . With the first report of West Nile virus in Germany in 2018 43 , the detection of authochthonous West Nile virus in the Netherlands is likely a matter of time. The readily established dead wildlife disease surveillance programme presented here and the expansion of the described protocol for sequencing on the Nanopore platform to West Nile virus should in this scenario allow for the real time monitoring of eventual reciprocal interaction in the dynamics of West Nile virus.

Methods
Dead blackbirds surveillance. Mortality in blackbirds was reported through a citizen science-based alerting system. For a selection of the cases reported, dead blackbirds were collected, autopsy was performed, and brain tissues were sampled for USUV diagnostics if autopsy indicated USUV as the possible cause of death until the first detection of USUV in 2016 30 , and systematically afterwards. The selection was based on freshness of the carcasses and on geographic location, with oversampling at locations where USUV activity had not been identified by then. Therefore, the sampling reflects the edges of observed bird mortality rather than the local evolution of the virus over time. Locations where dead blackbirds were found were registered by municipality, and date of death was based on the date of the sample collection.
Usutu virus diagnostics. Tissue from dead blackbirds was homogenized using the Fastprep bead beater (4.0 m/s for 20 seconds). Samples were spun down for 10 minutes at 10.000 xg after which Phocine distemper virus (PDV) was added as internal NA extraction control to the supernatant and total NA was extracted from the supernatant using the Roche MagNA Pure. The NA was screened for the presence of USUV using real-time PCRs described by Nicolay et al. 44 and Jost et al. 45  www.nature.com/scientificreports www.nature.com/scientificreports/ Multiplex pcR for nanopore sequencing. The multiplex PCR for MinION sequencing was performed as previously described 36 . In short, random primers (Invitrogen) were used to perform reverse transcription using ProtoScript II (NEB, cat. no. E6569) after which USUV specific multiplex PCR was performed in 2 reactions using Q5 Hot Start High-Fidelity DNA Polymerase (NEB, cat no. M0493). Nanopore sequencing was performed according to manufacturer's instructions using the 1D Native barcoding genomic DNA Kit (Nanopore, EXP-NBD103 and SQK-LSK108) on a FLO-MIN106 flowcell. A total of 12 or 24 samples were multiplexed per sequence run.
Data analysis Minion sequencing. Raw sequence data was demultiplexed using Porechop (https:// github.com/rrwick/porechop). Primers were trimmed and reads were quality controlled to a minimal length of 150 and a median PHRED score of 10 using QUASR 46 . A reference based alignment against an arbitrary chosen Usutu virus genome was performed in Geneious 47 or in CLC Genomic Workbench 11.0 (https://www. qiagenbioinformatics.com/). The consensus genome was extracted and compared to the non-redundant database using Blastn 48 . The most closely related sequence was selected and used for a second reference based alignment using the quality controlled reads. The consensus genome was extracted and positions with a coverage <100 were replaced with an "N". This threshold was set based on a recent publication demonstrating that by using this threshold the quality of full genomes is very high (one position per 106 USUV whole genomes sequences sequenced is called erroneously) 36 . Homopolymeric regions were manually checked and resolved consulting the closest reference genome. phylogenetic analysis. All available full length USUV genomes were retrieved from GenBank 49 on 12 February 2019 and aligned with the newly obtained USUV sequences using MUSCLE 50 . Sequences with >20% "Ns" were not included in the phylogenetic analysis. The alignment was manually checked for discrepancies after which IQ-TREE 51 was used to perform maximum likelihood phylogenetic analysis under the GTR + I + G4 model as best predicted model using the ultrafast bootstrap option with 1,000 replicates.
phycLip. PhyCLIP 37 was used to delineate the phylogenetic tree in clusters. The default settings for intermediate-resolution clustering were used with the recommended range of input parameters to determine the optimal parameters: S = 3-10 (increasing by 1), FDR = 0.05-0.20 (increasing by 0.05) and gamma = 1-3 (increasing by 0.5).
BeASt analysis. Bayesian phylogenetic trees were inferred using BEAST v1.10.3 52 . The HKY site model was used with 4 gamma categories with 3 partition into codon positions to generate an uncorrelated relaxed molecular clock. The tree prior was set to exponential growth and random sampling for the Africa 3 lineage and to constant size for the Europe 3 lineage. MCMC was set to 80,000,000 generations for both lineages. Log files were analyzed in Tracer v1.7.1 to check if ESS values were beyond threshold (>200). Tree annotator v1.8.4 was used with 10% burnin and median node heights. The tree was annotated and visualized using FigTree 53 .
Geographical distribution of USUV strains detected in dead blackbirds in the netherlands. The map was created using ArcGIS 10.6 software (ESRI Inc., Redlands, CA, USA) (http://www.esri.com/) and uses the datasets Gemeentegrenzen 2019 and Provinciegrenzen 2019 by ESRI Nederland (services.arcgis.com) for administrative boundaries.

Data availability
The genomic sequences of the Usutu viruses sequenced in this study have been deposited in the GenBank database under the accession numbers MN122145 -M122256.