# Relevance of the synergy of surveillance and populational networks in understanding the Usutu virus outbreak within common blackbirds (Turdus merula) in Metropolitan France, 2018

Malika Bouchez-Zacria, Clément Calenge, Alexandre Villers, Sylvie Lecollinet, Gaëlle Gonzalez, Benoît Quintard, Antoine Leclerc, Florence Baurier, Marie-Claire Paty, Éva Faure, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Malika Bouchez-Zacria, Clément Calenge, Alexandre Villers, Sylvie Lecollinet, Gaëlle Gonzalez, et al. Relevance of the synergy of surveillance and populational networks in understanding the Usutu virus outbreak within common blackbirds (Turdus merula) in Metropolitan France, 2018. Peer Community Journal, 2025, 5, pp.e9. 10.24072/pcjournal.513. hal-04903038

HAL Id: hal-04903038

https://ofb.hal.science/hal-04903038v1

Submitted on 21 Jan 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research article

Published
2025-01-21

**Cite as**
Malika Bouchez-Zacria, Clément Calenge, Alexandre Villers, Sylvie Lecollinet, Gaelle Gonzalez, Benoit Quintard, Antoine Leclerc, Florence Baurier, Marie-Claire Paty, Éva Faure, Cyril Eraud and Anouk Decors (2025)
*Relevance of the synergy of surveillance and populational networks in understanding the Usutu virus outbreak within common blackbirds (Turdus merula) in Metropolitan France, 2018*,
Peer Community Journal, 5: e9.

**Correspondence**
clement.calenge@ofb.gouv.fr

**Peer-review**
Peer reviewed and
recommended by
PCI Infections,
https://doi.org/10.24072/pci.infections.100236

(CC)

BY

This article is licensed
under the Creative Commons
Attribution 4.0 License.

# Relevance of the synergy of surveillance and populational networks in understanding the Usutu virus outbreak within common blackbirds (*Turdus merula*) in Metropolitan France, 2018

Malika Bouchez-Zacria<fnref n="1" />, Clément Calenge<fnref n="1" />,², Alexandre Villers³, Sylvie Lecollinet<fnref n="1" />,⁴,⁵, Gaelle Gonzalez⁵, Benoit Quintard⁶, Antoine Leclerc⁷, Florence Baurier⁸, Marie-Claire Paty⁹, Éva Faure¹⁰, Cyril Eraud<fnref n="1" />,³, and Anouk Decors¹¹

Volume 5 (2025), article e9

https://doi.org/10.24072/pcjournal.513

## Abstract

Usutu virus (USUV) was first isolated in Africa in 1959 and has since spread to and through Europe with a typical enzootic mosquito-bird cycle. In France, it was first detected in birds in 2015, but during summer 2018 the spread of USUV was particularly significant throughout the country, killing mainly common blackbirds (*Turdus merula*) and to a lesser extent great grey owls (*Strix nebulosa*), among other captive and non-captive wild bird species. Previous studies of USUV in France have focused on reconstructing pathways of introduction, but not on structural aspects of virus spread within the country. Data (RT-PCR of geolocated dead birds) on this 2018 outbreak were collected through both an event-based wildlife network named SAGIR and the health surveillance of the French-speaking Association of Zoo Veterinarians (AFVPZ). In addition, common blackbird populations could be monitored through another network (REZOP). Statistical analysis (spatial, temporal, spatiotemporal and environmental determinants) of the SAGIR and AFVPZ network data helped to highlight the early appearance of separate large clusters of USUV cases in mid-July 2018, the subsequent diffusion into smaller and secondary clusters at the end of August 2018, and a meanwhile enlargement of the first clusters with an increase in the number of cases. High human density (top 10.5% densest areas in France) and wetland concentration (top 19.3% most likely wetland areas) were significant factors in USUV case locations. Using generalised additive mixed models on REZOP data, we also highlighted the decline in common blackbird population trends in areas with medium and even more with high USUV pressure (areas defined based on SAGIR-AFVPZ data) following the 2018 outbreak (respectively -7.4% [-11.4; -3.9]⁹⁵% and -15.7% [-16.2; -9.1]⁹⁵%). A large area (radius ~150 km) in the centre and centre-west of France, and smaller areas in the south-east, north and north-east of France (each with a radius ~ 50 km) were particularly affected. We conclude on the importance to work with synergistic networks to assess infection spread in wild bird species, as well as the negative impact of an emerging arbovirus. The responsiveness of such a network system could be improved by automating alerts.

<footnote n="1">¹Independent researcher - Audincthun, France, ²Office Français de la Biodiversité - Direction Surveillance, Évaluation, Données - Unité Données et Appui Méthodologique - Saint Benoist, Le Perray en Yvelines, France, ³Office Français de la Biodiversité Direction de la Recherche et de l'Appui Scientifique (DRAS) Service Conservation et Gestion des Espèces à enjeux Villiers-en-Bois France - Villiers-en-Bois, France, ⁴CIRAD-INRAE, UMR ASTRE - Petit-Bourg, Guadeloupe, France, ⁵ANSES, INRAE, Ecole Nationale Vétérinaire d'Alfort, UMR VIROLOGIE, Laboratoire de Santé Animale - Maisons-Alfort, France, ⁶Parc Zoologique et Botanique de Mulhouse - Mulhouse, France, ⁷ZooParc de Beauval - St Aignan, France, ⁸Laboratoire Départemental du Cher - Bourges, France, ⁹Santé Publique France - Saint-Maurice, France, ¹⁰Fédération Nationale des Chasseurs - Issy les Moulineaux, France, ¹¹Office Français de la Biodiversité, Direction de la Recherche et de l'Appui Scientifique - Orléans, France</footnote>
Introduction

The Usutu virus (USUV), so called because of its isolation in *Culex neavei* mosquitoes in 1959 (McIntosh, 1985) near the Usutu river in Swaziland (Eswatini since 2018) is an arbovirus of the *Flaviviridae* family, *Orthoflavivirus* genus (Clé et al., 2019), belonging to the Japanese encephalitis virus serocomplex, and phylogenetically close to Japanese Encephalitis Virus (JEV) and West Nile Virus (WNV) (Calisher & Gould, 2003). After South Africa, the virus was detected in other African countries: Central African Republic, Senegal, Ivory Coast, Nigeria (Nikolay et al., 2011), Uganda (Nikolay et al., 2011; Mossel et al., 2017), Burkina Faso (Nikolay et al., 2011), Mali, Madagascar (Chevalier et al., 2020), Kenya (Ochieng et al., 2013), Tunisia (Ben Hassine et al., 2014; M'ghirbi et al., 2023), Morocco (Durand et al., 2016), and Israel (Mannasse et al., 2017). The first known occurrence in Europe dates back to 1996 in Italy (a retrospective finding in dead birds, mainly blackbirds) (Weissenböck et al., 2013). Since then, first detections of USUV followed in other European countries: Austria (Weissenböck et al., 2002) in 2001, Czech Republic (Hubálek et al., 2008) in 2004, Hungary (Bakonyi et al., 2007) in 2005, Poland (Hubálek et al., 2008), Spain (Busquets et al., 2008) and Switzerland (Steinmetz et al., 2011) in 2006, Serbia in 2009 (Lupulovic et al., 2011), Germany (Jöst et al., 2011), Greece (Chaintoutis et al., 2014) and Slovakia (Csank et al., 2018) in 2010, Croatia (Barbic et al., 2013) in 2011, Belgium (Garigliany et al., 2014) in 2012, France (Lecollinet et al., 2016) in 2015, Netherlands (Rijks et al., 2016) in 2016 and United Kingdom (Folly et al., 2020) in 2020 (Vilibic-Cavlek et al., 2020).

The analysis of 92 complete USUV genomes, including 77 genomes from mosquito, bird and bat species of Germany (2010-2014) highlighted that i) USUV can be classified in different lineages (Engel et al., 2016), (eight in number – three African and five European (Cadar et al., 2017; Clé et al., 2019)), ii) the most common ancestor emerged in Africa at the beginning of the 16th century, iii) USUV was regularly introduced from Africa during the last 75 years (a first introduction from 1950 through 1960's in Western Europe (Spain), a second one between 1970 and the 1980's in Central Europe (Austria) and a third one around 1996 in Western Europe (Spain)) and iv) in situ evolution could explain the genetic diversity of European lineages, while extensive gene flow drove African ones. In addition, flyway networks of migratory birds were strongly consistent with spatial movements observed through genetic data and suggested the possibility of USUV exportation from Africa i) to Spain via the east Atlantic and/or Black Sea/Mediterranean flyways and ii) to Central Europe through the Black Sea/Mediterranean flyway (Engel et al., 2016).

The natural life cycle of USUV involves passeriform and strigiform birds as amplifying hosts and ornithophilic species of mosquitoes as vectors (Clé et al., 2020). The virus has been detected in several Culicidae (*Culex*, *Aedes*, *Culiseta*, *Mansonia*) (Clé et al., 2019) but *Cx. pipiens* (experimentally shown to be USUV competent (Martinet et al., 2023)) is considered as the main vector in Europe (Becker et al., 2012; Martinet et al., 2023). To date, several incident or dead-end hosts (which certainly develop a low viremia, insufficient to re-infect mosquitoes (Martinet et al., 2023)) have been identified (Vilibic-Cavlek et al., 2020) : horses (Durand et al., 2016; Bażanów et al., 2018; Csank et al., 2018), dogs (Durand et al., 2016), rodents (Diagne et al., 2019), squirrels (Romeo et al., 2018), wild boar (Escribano-Romero et al., 2015; Bournez et al., 2019) and roe deer (Bournez et al., 2019). The virus has also been detected in bats, whose epidemiological role is as yet undetermined (Cadar et al., 2014). The virus has a zoonotic potential as humans (incident hosts) have been found infected with different clinical pictures: fever and rash in Africa (Ashraf et al., 2015), neuroinvasive infections (Gaibani & Rossini, 2017) in Europe (meningoencephalitis (Cavrini et al., 2009; Pecorari et al., 2009); idiopathic facial paralysis (Simonin et al., 2018)). However, these manifestations remain anecdotal and the asymptomatic form appears to predominate, as demonstrated by incidental findings during screening of asymptomatic blood donors and seroprevalence studies (Angeloni et al., 2023).

In Africa, USUV seems not to be pathogenic for local bird populations, either because it is not a naturally virulent virus, or because of immunity or genetic resistance of birds, which have been in contact with related flaviviruses for a long time (Bakonyi et al., 2004). In Europe, common
blackbirds are particularly affected, and to a lesser extent some birds of prey (Giglia et al., 2021), such as the great grey owl (Clé et al., 2019). The lesion patterns of 160 common blackbirds in the Netherlands included hepatosplenomegaly (major symptom), coagulative necrosis, lymphoplasmacytic inflammation and vasculitis. Reported symptoms included non-specific ones (immobility, apathy and ruffled plumage) and neurological ones (depression, stiff neck, inability to fly and epileptic seizures) (Giglia et al., 2021).

In several European regions, epidemiological patterns increasingly suggest endemic profiles rather than repeated introductions from other endemic countries (Constant et al., 2022), without clear explanation of efficient overwintering and subsequent amplification (Clé et al., 2019). Moreover, co-circulation with WNV is frequent (Beck et al., 2013; Fros et al., 2015; Constant et al., 2020). These situations have led European countries such as Italy to institute annual surveillance (entomological, human and veterinary), rather than focusing on a specific season (Constant et al., 2022).

In France, the first virus detection (in August and September 2015) took place in the North-east (Haut-Rhin department) and Central-east of France (Rhône department), respectively in two and three dead common blackbirds collected through SAGIR, the French event-based surveillance network collecting and analysing dead free-ranging wildlife animals (see details in next paragraph). Phylogenetic analysis helped demonstrating proximity between isolates i) from Haut-Rhin and Germany and ii) from Rhône and Spain. In addition, isolates identified in the two French departments belonged to different USUV lineages (Europe 3 in Haut-Rhin, Africa 2 in Rhône). Those results suggested a different source of introduction for the two French departments (Lecollinet et al., 2016). Since 2015, outbreaks in French avifauna have been recorded yearly in 2016, 2017 and 2018, with new departments involved (Johnson et al., 2018). In 2018, the Usutu virus circulated earlier and more extensively than it had in preceding years (Épidémiosurveillance Santé animale (ESA), 2018). In Camargue (South of France), which is located at the crossroad of several bird migration routes (Vittecoq et al., 2013)) two lineages (Africa 2 and Africa 3) were found in 2016 among female *Cx. pipiens*, suggesting two independent introductions of the virus: i) one from West to East and probably explained by migratory birds following the East Atlantic flyways (Africa 2 lineage) and ii) one explained by migratory birds following a central Mediterranean flyway and spreading after to Northern Europe, or originating from Germany (Africa 3) (Eiden et al., 2018). The USUV detection, in the southern region of Camargue, of same lineage, in same region and for consecutive years (2015, 2018 and 2020) in mosquitoes highly suggests endemicity (Constant et al., 2022).

Until now, studies on USUV in France focused on emergences and genetic characterisation to reconstruct introduction paths (Vittecoq et al., 2013; Lecollinet et al., 2016; Eiden et al., 2018; Constant et al., 2022), but not on spatial, temporal, spatio-temporal and populational aspects of the virus diffusion within the country. The aim of our study was twofold: i) from an epidemiological perspective, we first analysed the USUV diffusion among related bird populations to understand the way the disease had spread throughout mainland France and thus to define an epidemiologic pattern of this 2018 bird outbreak and its relationships with environmental variables (in particular wetlands, supposed to be a proxy for the density of mosquitoes); ii) from an ecological perspective, we tried to understand how this diffusion could have impacted common blackbird populations.

For the epidemiological perspective (both infectious and spatial approach), we used data from the SAGIR network which is dedicated to epidemiological surveillance, through the collection of dead or moribund birds and mammals (opportunistic sampling method). This participatory surveillance program aims at detecting as early as possible abnormal mortality or morbidity signals (Millot et al., 2017; Decors et al., 2022).

Moreover, in 2018, abnormal mortalities due to USUV have also been observed on captive birds in zoos and were reported in the database of the French Association of Zoo Veterinarians (AFVPZ). The French Agency for Food, Environmental and Occupational Health & Safety (Anses) combined SAGIR and AFVPZ data. This dataset offered an opportunity to study the epidemiology of USUV (date, location, and size of clusters), by identifying the timing and geography of reporting of dead birds positive to USUV in France.
We expected that in clusters with high density of dead birds identified from this dataset, the high mortality would have an impact on the population dynamics of wild birds. Common blackbirds, in particular, are of interest as this species is the most affected species by USUV in France. The impact of USUV on common blackbird populations can be assessed thanks to the data collected by another French observation network named "REZOP" (Réseau Oiseaux de Passage, passing birds network), which is dedicated to the evaluation of bird population dynamics and focuses on non-captive wild Alaudidae, Columbidae, Turdidae, Phasianidae, Corvidae and Sturnidae. In this program, birds belonging to the focus species are counted every year, through a systematic sampling method and collected data are statistically analysed to define population trends (Villers et al., 2021).

We hypothesised that USUV circulation in 2018 resulted in decreased common blackbird populations in the years following the 2018 epizootics, based on a previous study on blackbirds dead following USUV infection in Germany between 2011 and 2015 (Lühken et al., 2017). Although the approach of this German study was slightly different, as we will discuss in the dedicated part of this paper, we also expected that USUV circulation in birds would be associated with some environmental variables and would have a negative impact on common blackbird population trends.

# Material and methods

## Study design

This study employed an observational ecological design, using data from existing surveillance networks to assess the spatial and temporal distribution of USUV infections among bird populations in France. Data were collected opportunistically from mortality events reported through the SAGIR and AFVPZ surveillance networks, encompassing wild and captive birds. The study design allowed for the examination of environmental and ecological variables influencing USUV distribution, such as proximity to wetlands and human population density, which serve as proxies for mosquito habitats and bird sampling pressures respectively. We also identified spatial patterns (i.e., clusters) and temporal trends of USUV outbreaks at the population level. We used REZOP data to assess the correlation between estimated USUV infection levels and observed trends in blackbird populations across various geographic areas.

## Data collection and surveillance networks

The birds in our study were collected through the SAGIR network and zoological gardens across mainland France (Figure 1A) from July 15th to August 31st, 2018 (after this latter date, a message was sent to the local coordinators of the SAGIR network, asking them to limit data collection to only one sample per department, thus avoiding the additional cost of analyses linked to the increased epizootic outbreak; data collected by the network were thus defined at the department level from September onwards). Before August 31st, the most precise geographical information for each reported case was the municipality (a municipality being the smallest French administrative subdivision, corresponding to town and village areas, with a median area of 9 km²).

Data were provided i) mainly by the event-based surveillance network SAGIR and ii) by the AFVPZ surveillance of dead animals in zoos (with systematic necropsy and analysis of additional samples depending on observed lesions and ante-mortem clinical picture).

SAGIR is based on a collaboration between the French Office of Biodiversity (OFB), who administrates the network, hunters, local and national federations of hunters and the French Association of directors and executives of public veterinary analysis laboratories (Decors et al., 2022). In each *department* (a department being a French administrative district covering in average 5,700 km² [SD = 2,500 km²]), two departmental technical contacts (one belonging to the departmental service of OFB, the other belonging to the departmental federation of hunters) (Decors et al., 2022) coordinate a network of volunteers (professionals, as well as hunters, naturalists and farmers) who can report abnormal events (mortality or morbidity) that they
accidentally find in the field. The reported events can also include dead animals of local rescue centres on rare occasions. The sampling of these events is therefore opportunistic. For traceability purposes, a standardised individual digital form is filled for each mortality event collected, reporting epidemiological, agricultural and ecological circumstances surrounding the discovery of carcasses. When available, clinical signs are also provided (Millot et al., 2017), deduced from observation of the carcass and signs left in the environment. Sometimes animals are seen alive before they die (e.g. 12% of common blackbirds in the database from 2014 to 2023, regardless of cause of death) (A. Decors, personal communication).

<figure><img src="image_3.png" /><figcaption>**Figure 1** - Study area (grey) (A) with regions and towns mentioned in the study, wetlands (B) and human population densities (C) (B: from black to pink: from the least likely to the most likely wetlands; C: from black to pink: from the least to the most densely populated areas; B and C: yellow surrounded by a yellow line: main towns and regions mentioned in the study)</figcaption></figure>

## Laboratory procedures

In SAGIR, carcasses (fresh, chilled or frozen) accompanied by their form are addressed to the local administrative laboratory of veterinarian analyses where necropsies are realized. Following a gross pathologic examination, relevant aetiology tests are conducted: parasitology, bacteriology, virology, mycology, toxicology and/or histology (Millot et al., 2017).

In AFVPZ network, dead zoo animals are systematically necropsied within a maximum of 20 hours (usually less than six hours after death). The zoo veterinarian performs the necropsy and takes biological samples when lesions are observed, for analysis (bacteriological, mycological, parasitological, or virological if relevant). The samples are sent to the nearest departmental veterinary laboratory or, in some cases, to a laboratory specialised in the species concerned. Histological analysis is also carried out almost systematically.

Birds collected from both networks were also tested by the National Reference Laboratory (NRL, Anses), if i) their ante-mortem clinical picture was suggestive of an USUV infection (e.g. neurological symptoms) and/or ii) the epidemiologic context was consistent with a USUV suspicion (close mortality among sensitive species, and especially great grey owls) and/or iii) lesions at post-mortem examination were suggestive of an USUV infection (liver, spleen, cerebral necrosis). Spleen, liver, brain and lung (sometimes kidney and heart) samples were then sent to the NRL and tested for USUV and WNV through Real-Time Polymerase Chain Reaction (hereafter RT-PCR). Tissues samples were stored at -80°C before RT-PCR analyses and were grinded in Dulbecco modified Eagle's minimal essential medium (DMEM) (Lecollinet et al., 2016) with ceramic beads (MP Biomedicals, Illkirch, France) and FastPrep ribolyzer in BSL3 facilities. Automated total RNA extraction, amplification and detection of USUV genome were carried out as described in (Moutailler et al., 2019).
## Epidemiological data analysis

In our study, we focused only on birds reported USUV infected, and we discarded other reported birds. We found it easier to analyse the space and time structure with point patterns (availability of a large variety of methods to deal with such data) so that we transformed our dataset into a point pattern. Thus, for each reported bird, we randomly attributed a precise point location by randomly selecting a point within the municipality where it was collected. Hereafter, the terms 'points' or 'cases' will refer to these birds. All USUV-infected birds were considered, whatever their species or origin (captive versus non-captive).

We analysed the temporal, spatial and spatio-temporal distribution of reported cases, as well as the relationship with spatial clusters of cases with wetlands and with human density (Figure 1B and C). Wetlands and human density were supposed to be a proxy for the density of *Culex pipiens* mosquitoes (Becker et al., 2010; Haba & McBride, 2022). Human density was also supposed to be a proxy for common blackbirds sampling pressure, as those birds are highly commensal with humans and easy to detect in gardens. A map of potential wetlands was derived from a raster map of continental France built by INRAE (Orléans) (http://geowww.agrocampus-ouest.fr/web/?p=1538): the INRAE map estimates the probability that each 1km x 1km pixel of the map harbours a wetland and then discretize this probability in three classes (0 = no wetlands, 1 = rather strong probability, 2 = strong probability, or 3 = very strong probability). This map also identifies lakes and foreshores in France. We downloaded this map, and redefined pixels with lakes and foreshores as "confirmed wetlands" (new class 4). Then, we smoothed this map using a sliding window, by calculating for each pixel the sum of this ordered factor (0 to 4) for the focus pixel as well as the pixels immediately on the left, the right, the top and the bottom of the focus pixel (leading to a score comprised between 0 and 20 giving a smoothed probability of wetland in a given place; see Figure 1B).

A map of human density was calculated by finding the municipality in which the centre of each 1km x 1km pixel was located and attributing the 2005 human density to it (number of inhabitants divided by municipality area; see Figure 1C). We then log-transformed this map (transformation: $x \to \log(x + 1)$ where x was the number of inhabitants per km²). The 2005 human density was provided by Institut Geographique National (France) (IGN) (https://geoservices.ign.fr).

We then carried out a statistical analysis of the SAGIR and AFVPZ data aiming at identifying spatial and spatio-temporal patterns. All our analyses were carried out with the R software (R Core Team 2022). We have programmed an R package named usutuFrance, available at https://github.com/ClementCalenge/usutuFrance, containing all the code and data used to fit the model. It can be installed in R with the package devtools, using the function `devtools::install_github("ClementCalenge/usutuFrance", ref="main")`. This package includes a vignette describing how the user can reproduce the calculations carried out in this paper (vignette available with the command `vignette("usutuFrance")` once the package has been installed and contains supplementary analyses that helped to understand the structure of our data. This vignette serves as the supplementary material of our paper.

## Spatial analysis

We first used a Ripley's K function to identify the scales at which clusters of points could be identified (Diggle, 2013). This function allows the characterisation of second order point pattern properties: it is proportional to the mean number of other cases expected in a x km radius around a typical case of this point pattern (so that the larger it is, and the more there are cases in the neighbourhood of a case in average – in other words, clusters can be observed at this scale). The K function was computed for all distances *t* comprised between 0 and 250 km on the observed point pattern ($\hat{K}(t)$); the theoretical value of this function was also calculated for the same range of distances under the hypothesis of complete spatial randomness (CSR) ($K_{CSR}(t)$). Thus, under the hypothesis that the cases are randomly distributed in space, the difference $\hat{L}(t)$ between $\hat{K}(t)$ and $K_{CSR}(t)$ should be equal to 0 (Baddeley et al., 2014), so that positive deviations of $\hat{L}(t)$ from zero reflect clustering in the data (and conversely, negative values reflect repulsion mechanisms).
We also built confidence envelopes around simulated $L_r(t)$ functions expected under CSR: we repeated $n=100$ simulations of the CSR (by randomly distributing $N$ cases in France, where $N$ is the observed number of cases in our study), and for each simulation $r$, we estimated the function $L_r(t)$ for all values of $t$ from 1 km to 250 km. We plotted the observed value $\hat{L}(t)$ and the distribution of values expected under the CSR hypothesis on the same graph to compare them ($\hat{L}(t) = L_r(t)$: complete random reporting of cases; $\hat{L}(t) > L_r(t)$: more points within a $t$ radius around a typical case than under CSR hypothesis; $\hat{L}(t) < L_r(t)$: less points within a $t$ radius around a typical case than under CSR hypothesis i.e. avoidance process).

To identify areas of high density of cases in the point pattern, we carried out a kernel smoothing of the point pattern, using a smoothing parameter of 125 km (Wand & Jones, 1995), to visualize them on a map. The results (see Results section) suggested that the point pattern could be adequately described by a Thomas process (Diggle, 2013), i.e. a process in which cases are clustered in the following way (Baddeley et al., 2014): i) a random number of “parents” is generated from a Poisson distribution parameterized by a parameter $\kappa$ controlling the density of ‘parent’ cases per km². These parents are randomly placed in the study area. ii) For each “parent”, a random number of “children” (corresponding to the cases in our study) is generated from a Poisson distribution parameterized by a parameter $\lambda$ controlling the mean number of “children” per “parent”, and iii) for each “parent”, the children locations are randomly distributed around the parent’s location. More precisely, the locations of the children are supposed to be located at a distance from the parent randomly drawn from a semi-normal distribution with a standard deviation $\sigma$, with an angle from west direction drawn from a uniform distribution bounded by $\pi$ and $-\pi$. We fitted this process to our dataset using the method of minimum contrasts based on the function $K$. Estimated parameters of the fitted model helped defining the point pattern typical structure and the product of $\lambda$ and $\kappa$ provided the cases density per km². Note that our spatial analysis revealed a clustering of the dataset at very small scale (see results), possibly due to the tendency of observers to look for other cases in a municipality where a first case was identified. Consequently, and to improve model fit, we decided to thin the point pattern prior to the fit to remove this artefactual structure, keeping only the first bird found when two birds were collected at a distance lower than 5 km (see Supplementary material for more details).

## Environmental determinants

Then, we used this fitted model to test, using a randomization approach, whether the mean value of the two environmental variables (wetlands and population density) were greater in the places where the dead birds infected by USUV were found than expected by chance (Manly, 1991): i) we first calculated the means of the density of wetlands and the mean of the log of human population density within the observed points pattern, ii) we simulated the fitted Thomas process 999 times over France, iii) we calculated for each simulation the mean of the density of wetlands and the mean of the log of human population density (simulated distribution), and iv) we compared the observed values (means) with the simulated distribution (means and standard errors) and determined p-values with the help of a bivariate randomization test.

## Temporal analysis

To define the temporal evolution of the number of cases throughout the study period, we computed the mean number of cases per day during this period and smoothed the curve using moving average approach (Diggle, 1990): for each $d$ day of the study period, we calculated the mean number of cases per day, over a week centred on $d$. We plotted the obtained values on a graph to visualize temporal fluctuations.

## Spatio-temporal analysis

We computed a space-time $K$ function ($\hat{K}(u, v)$) (Gabriel et al., 2013) giving the expected number of cases within a distance of $u$ km around a randomly sampled case of the observed point pattern (for all values of $u$ comprised between 0 and 250 km), and which occurred at most $v$ days after this case (for all values of $v$ comprised between 0 and 10 days). To test whether the spatio-
temporal patterns identified by this function could have been obtained by chance, we also used a randomization approach. More precisely, for each one of the 1,000 simulations of this approach, we randomly permuted the dates associated to the reported birds (i.e. the cases order was randomised), thus conserving the spatial and temporal point pattern structures but randomising the association between space and time. For each simulation, we calculated the related $K(u, v)$ on the randomized dataset. We then compared, for each distance $u$ and each duration $v$, the observed value of $\hat{K}(u, v)$ to the distribution of simulated values $K(u, v)$. Probabilities that a simulation could lead to $K(u, v) > \hat{K}(u, v)$ were plotted, with duration on the x-axis, distances on the y-axis (probabilities near zero indicating clustering) and a contour limit identifying the set of distance/duration pairs for which the proportion of simulation with $K(u, v) > \hat{K}(u, v)$ was lower than 5%. This allowed to identify the space-time clusters in our dataset.

## Impact on common blackbird population densities

### Bird survey and data collection

The abundance of common blackbirds was monitored in France over the 1996-2022 using a large-scale survey designed to assess population trends of common game species (*ACT Survey* (Boutin et al., 2001, 2003)). The spatial coverage of this survey is based on 1,070 grid cells measuring 28 x 20 km (1/50 000e IGN map; see Supplementary material). Each cell includes a sampling route (~4 km in length) randomly selected outside highly urbanized areas, with five points counts spaced by ~1km. Each year, these points are surveyed twice during the breeding season by hundreds of observers (affiliated to OFB or departmental federation of hunters): from April 1st to 30th and from May 15th to June 15th. The census protocol adopts the point count methodology: each point is surveyed for ten minutes in early morning. During this period, observers record the number of different singing males within a 500 m radius circle around the point (see Villers et al. 2021 for more details). Over the 1996-2022 period, each route was surveyed on average for 24.6 years (SD +/- 3.8 years; range: 1-27) and each year, the number of surveyed route was on average 977.1 routes (SD +/- 69.5; range 750-1038).

### Population trends

For a given route and survey date, the number of singing birds was summed across the five points. We ensured that for each observation, all the following conditions were met: i) time of the count was available, ii) the start time of the survey was between -60 mn and +180 mn relative to sunrise time and iii) the date of the observation was posterior to the 31st of March and anterior to the 30th of June. We used these date and time intervals as presumed windows that minimize biases related to the detection rate of singing males. Our final dataset included 1,070 different routes (see Supplementary material), totalizing 50,454 counts over the study period.

We used the framework of Generalized Additive Mixed Models (GAMM) implemented in the R package mgcv (Wood, 2006) to calculate annual population indices and underlying population trends. The number of singing males per route and per date (*Counts*) was modelled as a function of a smoothed function (using thin plate regression splines) of the counting *Date* (expressed as a Julian date, i.e. the number of days elapsed since the beginning of the year), a smoothed function of *Time-of-day* (expressed in minutes relative to local sunrise thanks to the R package suncalc (Thieurmel & Elmarhraoui, 2022), and a random effect associated to *Route ID* (a 1,070-levels factor). We used *Date* and *Time-of-day* as smooth terms to capture some components of the detection process (Knape, 2016), including plausible change in singing rates in relation to breeding phenology and/or singing behaviour during the morning. In order to straightforwardly compare trends in the abundance index according to different levels of USUV incidence, we also created a factor *UsutuArea*, indicative of the severity of the 2018 USUV episode in space. For each point of the *ACT Survey* we extracted the value of the density of USUV cases estimated by the kernel smoothing carried out previously (cf. supra, section USUV diffusion process, Spatial analysis). These values were averaged for each *route* (mean value of the five sampling points), and converted into a 3-levels factor, according to the 50th and 75th percentiles of the distribution of USUV cases density: *routes* displaying kernel values below the median were coded as *low*, over
the 75th percentile as high, and medium otherwise. In order to model the temporal variation of the abundance index in each area, we included the UsutuAre a as a fixed effect, a smooth function of years (therefore considered as continuous covariate, Year, with a thin plate regression spline), and an interaction term allowing to fit a temporal trend to each level of the UsutuAre a factor (with independent smoothing parameters (Pedersen et al., 2019)). This interaction allowed to identify different population trends corresponding to different levels of USUV cases density. Finally, we also included in the model a random effect associated to the year considered as a 27-levels factor (Year.as.Factor). This statistical approach allowed separating short-term fluctuations (from one year to the next) from longer-term (over a decade) population changes, through the inclusion of temporal random effects in simulation of smooth trends (Knape, 2016). The number of degrees of freedom of the basis function of the smoother of the year effect was set equal to k=10.

Models were fitted with a Tweedie distribution. Tweedie distributions, a special cases of exponential dispersion models, are probability distributions encompassing the continuous (normal, gamma and inverse Gaussian), discrete (Poisson) and compound (Poisson-gamma) distributions. In the R package mgcv, the value of the power parameter p of the Tweedie distributions (with 1 < p < 2) can be estimated simultaneously with the estimation of other parameters, offering flexibility in the type of response that can be modelled.

To compute the trends and their associated confidence intervals from the fitted model, we adopted a posterior simulation sensu (Wood, 2017) (chapter 6.10, p. 293), in order to obtain the distribution of predicted values of abundance for each level of the UsutuAre a factor. A similar approach is implemented in the R package poptrend (Knape, 2016), but it does not offer the possibility to compute trends for specific models, e.g. aiming at comparing trends for different levels of a grouping variable. The general idea is to combine prediction matrices and simulation (replicate) from the statistical distribution of the parameters β of the model, in order to obtain distribution estimates (i.e. including uncertainty of parameter estimate) of any quantity of interest, in our case the relative abundance index of common blackbirds for an average route and year, hereafter referred to as the variable Index. Details of the computations can be found in the supplementary material.

Trends were expressed as:

$$ (1) \qquad \%change = 100 \times \frac{(Index_{endsurvey} - Index_{startsurvey})}{Index_{startsurvey}} $$

To estimate population trends accounting for the variability of the abundance index, we computed, for each simulated replicate, the mean abundance index for the four years before (2015-2018) and after (2019-2022) the USUV episode of the 2018 summer. We then plotted for each level of the UsutuAre a factor the value of the abundance index for the two periods (2015-2018 and 2019-2022) and compared the population trends for each area thanks to Equation (1).

# Results

## USUV diffusion process

### Epidemiological surveillance network and provided data

The dataset (60 % SAGIR, 34% AFVPZ and 6% rescue centres) included 60 birds tested for USUV of which 50 were detected infected in 29 French departments between July 15th and August 31st, 2018. Wild common blackbirds and captive great grey owls were particularly affected even if several other species were concerned (Table 1).
<table><caption>**Table 1.** Distribution of birds collected and detected infected by USUV within the two epidemiological surveillance networks of the study between July 15th and August 31st, 2018 (data from rescue centres were stored in the SAGIR database)</caption><thead><tr><th rowspan="2">Species</th><th colspan="3">Origin</th></tr><tr><th>SAGIR</th><th>AFVPZ</th><th>Rescue centres</th></tr></thead><tbody><tr><td colspan="4"><strong>Wild birds</strong></td></tr><tr><td>Common blackbird (Turdus merula)</td><td>26</td><td>1</td><td>3</td></tr><tr><td>Blue tit (Cyanistes caeruleus)</td><td>2</td><td>0</td><td>0</td></tr><tr><td>Robin (Erithacus rubecula)</td><td>1</td><td>0</td><td>0</td></tr><tr><td>Song trush (Turdus philomelos)</td><td>1</td><td>0</td><td>0</td></tr><tr><td colspan="4"><strong>Captive birds</strong></td></tr><tr><td>Great grey owl (Strix nebulosa)</td><td>0</td><td>12</td><td>0</td></tr><tr><td>Snowy owl (Bubo scandiacus)</td><td>0</td><td>2</td><td>0</td></tr><tr><td>Strigidae (unspecified species)</td><td>0</td><td>1</td><td>0</td></tr><tr><td>Western capercaillie (Tetrao urogallus)</td><td>0</td><td>1</td><td>0</td></tr></tbody></table>

## Spatial analysis

The comparison of observed $\hat{L}(t)$ and the distribution of simulated $L(t)$ (n=100) showed that infected birds were not randomly distributed and were concentrated in specific geographic areas: it allowed to identify a clustering of the point pattern at 5 km, 125 km and 200 km (Figure 2A).

The kernel smoothing carried out with a smoothing parameter equal to 125 km showed a high density of reported cases in the Brenne area (largest cluster located at the southwest of Paris on Figure 2B – Paris being the northernmost cluster delimited by contour lines on the figure (see Figure 1A for the names of towns and regions). Other high cases density areas could be identified, particularly in Paris, towards Orléans (south of Paris) and Nantes (west of France; Figure 1A, Figure 2B). The distance between Paris and Orléans high USUV cases density areas is around 200-250 km, as well as between those of Brenne and Nantes. This last result was consistent with the third peak observed in Figure 2A (clustering at 200 km).

The fitting of the Thomas process model helped determining the density of clusters (“parents”) per 10,000 km² (κ parameter: 0.17) and the mean cluster size (λ parameter: 4.36, mean number of “children” per “parent”). Thus, the mean density of detected cases was estimated to be 0.74 cases per 10,000 km². The standard deviation of the distribution of distances between the cases and the centroid of the cluster (the parent case) was of 77 km. This last result was consistent with the K function’s second peak, which indicated the presence of clusters at a scale of about 130 km. Indeed, this standard deviation measures the typical distance between cases and the centroid of the cluster, while the K function is based on the distance between two cases. By simulation, we checked that the mean distance between points sampled in a bivariate Gaussian distribution with a standard deviation of 77 km was close to 130 km (see Supplementary material).

## Environmental determinants

The observed mean of wetland probability score in USUV places was equal to 4.81, while the one expected under the absence of effect of the environmental variables (for n=999 simulations of the fitted Thomas process) was of 2.69 (SD=0.65). The difference was significant (p-value=0.007), indicating that USUV cases were more often reported in or close to wetland than expected by chance. This mean wetland probability score in USUV places corresponded to the top 19.3% most likely wetland areas. Similarly, the observed mean of human log-density population in USUV places was of 5.02, while the one expected under the hypothesis of absence of effect of the environmental variables was of 3.50 (SD=0.25). The difference was also significant (p-value=0.001), indicating that places where dead birds positive to USUV were reported were more
often located in places with high human density. This mean log-density population in USUV places corresponded to the top 10.5% densest areas in France.

## Temporal analysis

There was an initial increase in the number of cases between July 15th and 20th 2018 (Figure 2C). Then from July 20th to August 15th the number of cases reported per day stagnated at about one per day. From August 15th onwards, there was a sharp increase, with an average of more than two cases reported per day towards the end of the month.

<figure><img src="image_4.png" /></figure>

**Figure 2** - Spatial (**A**, **B**), temporal (**C**) and spatio-temporal (**D**) characteristics of the point pattern of birds detected infected to USUV during the study period (**A**: Comparison between $\hat{l}(t)$ (observed, grey line) and $L_r(t)$ (envelope of simulated point pattern under CSR hypothesis, grey area); **B**: Kernel smoothing of the point pattern (grey lines) and birds locations (black crosses); **C**: Mean number of cases (mobile average approach) reported per day and depending on date; **D**: Probabilities that spatio-temporal K function $K(u, v)$ simulated under the hypothesis of independence between spatial and temporal patterns, was higher than observed spatio-temporal K function $\hat{K}(u, v)$, depending on duration and distance (white line: contour limit of the set of distance/duration pairs for which the proportion of simulation with $K(u, v) > \hat{K}(u, v)$ was lower than 5%)
## Spatio-temporal analysis

The plot of probabilities that a simulation could lead to $K(u, v) > \hat{K}(u, v)$ showed a strong clustering at small distances of about 10 to 20 km as soon as two or three days after a case, expanding to 50-60 km five days after a case. But we also observed a clustering on larger distances (150 km) starting as soon as two or three days after a case, which increased to reach a large range of distances after a week (90 to 150 km) (Figure 2D). The large distances clustering was consistent with the first phase of the diffusion (between July 15th and 20th) (Figure 2C), with the emergence of the main clusters of Brenne, Paris, Orléans and Nantes (Figure 2B). The small distance of clustering was consistent with the emergence of secondary clusters and the expansion of preexisting main clusters (Figure 2B) during the third phase of the diffusion (from August 15th) (Figure 2C).

## Impact on common blackbird population densities

The means and confidence intervals depicting cells of respectively low, medium and high USUV cases density have been reported in Table 2, the distribution of those cases in Figure 3A and their spatial distribution in Figure 3B.

The estimation of common blackbird population trends estimated through variations of the mean population indexes for 2015-2018 and 2019-2022 periods was not statistically different from 0 for the low USUV area but presented a significant decline for both medium and high USUV areas, being more pronounced for the latter than for the former (Table 2). Two main results were highlighted: i) a decreasing trend for common blackbird populations that have been a priori exposed to USUV since the 2018 USUV outbreak and ii) the greater the pressure of infection was, the greater the negative effect on common blackbird populations was (Figure 3C).

<table><caption>**Table 2.** Trends computed between 2015-2018 and 2019-2022 for the three different levels of density of USUV common blackbird cases (CI: Confidence Interval).</caption><thead><tr><th>USUV area</th><th>Mean trend</th><th>2.5%CI</th><th>97.5%CI</th></tr></thead><tbody><tr><td>low</td><td>+0.8%</td><td>-0.4%</td><td>+2.1%</td></tr><tr><td>medium</td><td>-7.4%</td><td>-11.0%</td><td>-3.9%</td></tr><tr><td>high</td><td>-12.7%</td><td>-16.2%</td><td>-9.1%</td></tr></tbody></table>

## Discussion

Birds affected by USUV between July 15th and August 31st, 2018, in France were mainly common blackbirds and to a lesser extent captive great grey owls. Cases were clustered at i) 5 km, ii) 125 km (which was the average distance between two typical cases which we highlighted with a Thomas process model) and iii) 200 km (which was the distance between main big clusters which appeared simultaneously at the beginning of the outbreak and expanded thereafter). Three periods of USUV diffusion have been featured: i) a first increase of USUV cases between July 15th and 20th (with clustering at large distances of about 150 km between two and three days after a case and between 90-150 km a week after a case, corresponding to the emergence of main clusters – Brenne, Paris, Orléans and Nantes), ii) a stagnation between July 20th and August 15th, and iii) an important increase from August 15th (with strong clustering at small distances of about 10-20 km, between two and three days after a case and expanding to 50-60 km five days after a case, corresponding to the emergence of secondary clusters and the expansion of preexisting main clusters). Detected case locations were statistically associated with wetlands and high human log-density population areas. Common blackbird population trends estimated suggested an impact of the USUV infection on them, with a decrease of the mean population index over 2015-2018 compared to 2019-2022. Areas of a priori greater USUV pressure presented more pronounced negative trends than the two other areas.
<figure><img src="image_5.png" /><figcaption>**Figure 3** - Common blackbird population trends in response to the USUV episode (A: Histogram of USUV cases density, from a kernel smoothing, for the routes where common blackbird were monitored (solid line: median of the distribution, dashed line: 75% percentile of the distribution); B: Maps of the kernel smoothing of USUV cases and locations of the centroids of REZOP road (pale green: low, orange: medium and dark purple: high densities of USUV cases as estimated by the kernel method); C: Predicted abundance index (mean values and their 95% confidence intervals) for the three classes of USUV cases density (pale green: low, orange: medium and dark purple: high densities) before (2015-2018) and after (2019-2022) the 2018 USUV outbreak. Each solid line represents changes in abundance value for a given simulation (1,000 simulated values were generated).</figcaption></figure>

Species involved in the 2018 outbreak in the present study were those also found in other European countries such as Germany (Becker et al., 2012; Lühken et al., 2017). Several hypotheses might be put forward regarding the overrepresentation of common blackbirds: i) they are not trans-Saharan migratory species and as a consequence probably deprived of an immunity they could have acquired in Africa, ii) they might be more genetically susceptible than other birds species, iii) African strains introduced in Europe might have adapted to common blackbirds and become more virulent for this host (Bakonyi et al., 2004), iv) co-infections with *Plasmodium* might
exist with a possible interplay of the two agents, as highlighted in the Netherlands (Agliani et al., 2023), v) a feeding preference has been identified for them (when they were sufficiently abundant) by *Cx. pipiens* (Rizzoli et al., 2015) and vi) this species is abundant, ubiquitous, human commensal (Lühken et al., 2017) and thus easy to detect by the SAGIR network.

Any pattern identified by the analysis of data collected by a participatory program relying on non-controlled sampling pressure can be caused by the opportunistic observation process (characterized by a very heterogenous sampling pressure) as well as by the ecological or epidemiological process. In particular, the data collection process differs across departments (differences in detectability, in probability to report a carcass to the network, in probability to collect this carcass by the local contact, mode of transport, mode of storage, etc.). Thus, when a pattern was identified in the data (e.g. a cluster of cases), we could not decide between (i) the fact that this pattern was caused by the observation process (e.g. higher local sampling pressure due to higher human density and thereby detectability) and (ii) a biological cause (e.g. local factors favouring the diffusion of USUV in this place). It would have been interesting to model the data collection process to draw firmer conclusions on the biological process. However, the lack of information about the data collection process within the SAGIR network has been a long-standing issue, hindering such modelling (as is often the case with participatory networks). Our approach was therefore exploratory, and our results could only help formulating hypotheses, and could not allow confirming them. Keeping in mind this limit, we must stress that our results pointed out that such a network and statistical analysis of admittedly imperfect data could however allow the understanding of diffusion pattern and epidemiological process. Indeed, the network was alerted to the possible presence of USUV cases in France and the level of attention given to surveillance was greater for this disease than for other diseases. We considered the possibility that any pattern identified in our data could be caused either by the observation process (e.g. the heterogenous sampling pressure) or by ecological/epidemiological processes.

The short distance clustering (at 5 km) identified by our analysis could be explained by local sampling effort. We could not exclude that observers would have been more proactive in detecting dead birds around first detections. However, this short distance of 5 km is also close to flight distances of *Cx. pipiens* estimated in a study carried out in Chicago, Illinois, on *Culex* female mosquitoes. Indeed the mean dispersal distance was 1.15 km, and 90% of individuals stayed within 3 km from their larval habitat (Hamer et al., 2014). Another explanation for the 5 km clusters we observed would be a USUV diffusion mainly due to mosquitoes. Sufficient amplification of the virus in its various hosts would also have allowed a large number of mosquitoes to be infected.

The occurrence of different clusters simultaneously at the beginning of this 2018 USUV outbreak, far apart from each other, is in favour of USUV endemicity throughout the country. The endemicity has already been proven in Southern (Constant et al., 2022) and North-Eastern France (Johnson et al., 2018), but the 2018 outbreak we studied and the diffusion pattern we showed raised the hypothesis of a national endemicity. Previous outbreaks, even of lesser magnitude, might have enabled certain lineages to settle. Indeed, in France, USUV Europe 3 strain has been identified in common blackbirds in two north-eastern departments between 2015 and 2017. USUV Africa 2 strain has been isolated in common blackbirds, mosquitoes and humans in three south-eastern departments and in Rhone between 2015 and 2016, while USUV Africa 3 strains have been identified in mosquitoes and birds in two central and four south-eastern French departments between 2015 and 2017 (Johnson et al. 2018, NRL data). Favourable conditions (e.g. sizes, naive immune status and/or turnover of bird populations, temperature, rainfall, human population movements) could then lead to new outbreaks thereafter. The early detection of this 2018 outbreak was made possible by the well-informed SAGIR network, which vigilance has been reinforced since the 2015 outbreak for USUV relevant species. Throughout Europe, endemic circulation and spread of USUV infections have been reported, e.g. Northern Italy, Hungary or Austria (Constant et al., 2022). The study of lineages involved within this 2018 outbreak could help validate endemicity.

The different phases of the 2018 USUV diffusion we observed remained difficult to explain: a first increase of cases, with clusters far from each other, a stagnation and then an important
increase of cases, with secondary clusters and the enlargement of the first ones. If endemicity
might explain the cases pattern of the first phase, the origins of the third phase remain unclear.
Meteorological variables, as mean temperature (Hamer et al., 2014), could explain this result.
Indeed, the duration of Cx. pipiens eggs development depends on temperature (hatching after only
one day at 30°C, three days at 20°C, ten days at 10°C and incomplete embryonic development
below 7°C (Becker et al., 2010)). Temperatures were particularly high in July and August 2018 in
France (respectively +2.18°C and +1.1°C higher than the 1991-2020 baseline average
(https://www.meteocontact.fr/climatologie/france/bilans-climatiques)). Cx. pipiens density might
have increased in August due to these favourable conditions and promoted an increased possibility
of contact between common blackbirds and the arthropod. Moreover, laboratory experiments have
shown that higher temperatures resulted in higher infection rate of Cx. pipiens (Fros et al., 2015)).
These conditions could have allowed an important development of Cx. pipiens and the third phase
of USUV diffusion we could observe. Further studies collecting local relevant meteorological data
would be useful to assess this hypothesis. Nevertheless, we could not exclude that the diffusion
of information and knowledge of the presence of USUV in France would have also led to better
sampling effort from August 15th.

In addition to the above-mentioned hypothesis related to the population dynamics of the main
vector of the disease, other non-mutually exclusive hypotheses could explain the appearance of
USUV clusters in late summer (after August 15th). This period coincides with the initiation of
movements in common blackbirds, particularly juveniles that disperse outside their parents'
territories (Snow, 1958), and to a lesser extent, with the first mentions of post-breeding migrants
(Toulotte et al., 2022). Similar movements have been shown as the main drivers of seasonal
epizootics of the avian influenza in mallards (van Dijk et al., 2014), and the same processes might
be at play in the USUV episode reported here. USUV naïve birds such as juveniles, or individuals
whose original populations have not yet been exposed to USUV, might be more prone to develop
and spread the virus, easing the apparition of clusters.

In resident blackbirds, this period also overlaps with resource-demanding activities including
breeding for some late breeders (Sauvage, 2016), post-breeding (complete) moult in adults (Snow,
1969; Morrison et al., 2015) and partial moult in juveniles (Cramp, 1988). Previous works have
shown that such activities may induce trade-offs with immunity in birds (e.g. Moreno et al. 1999;
Sanz et al. 2004; Martin II 2005; Moreno-rueda 2010). Accordingly, late summer may coincide with
a period of greater susceptibility to infections in blackbirds. Unfortunately, the lack of available data
on age, sex, or moulting status, both for the recorded USUV cases, but also for individuals in the
exposed populations, did not allow us to confirm these assumptions. This constitutes a relevant
area of investigation for future studies and highlights the need for a deeper description of
specimens found dead.

The fitted Thomas process was a simple model, with the assumption that all clusters were of
the same size, that the mean number of cases per clusters did not vary from one cluster to another
and that the density of "parents" generating secondary cases was uniform across the study area.
However, we could not highlight any lack of fit of this model to our dataset, and by simulating this
process we were able to simulate a realistic point process and thereby test the effect of
environmental variables on the occurrence of USUV. This test showed that the density of USUV
cases varied as a function of wetlands and human density, two variables that we supposed to
reflect mainly mosquito abundance.

USUV is assumed to be mainly transmitted among birds by mosquitoes during their blood meal (Lühken et al., 2017), but a transmission from bird to bird could not be excluded in the diffusion process and the endemicity. Oral route transmission was suggested by damage of cells of the gastrointestinal tract observed as part of necropsies of infected common blackbirds and possible vertical transmission could not be excluded with antigen detection in the genital system (Giglia et al., 2021). An oral transmission could explain distances between cases higher than maximum flight distances of mosquitoes, while vertical transmission could be involved in endemicity mechanisms and perhaps even in the resurgence of cases after the reproduction period (from August 15th).
The association between wetlands and USUV was consistent with the life cycle of USUV, with mosquitoes as vectors and wetlands known as suitable habitat for *Cx. pipiens* (Becker et al., 2010; Haba & McBride, 2022). The association with high human population densities could be more surprising at first sight. In an event-based network, the number of observers could partly explain such a result. Indeed, the greater the number of people able to collect birds affected by the disease would be, the greater the potential for the disease to be detected in densely populated areas would be. In addition, common blackbirds are familiar to urbanised areas (Partecke & Gwinner, 2007), and their size and colour (especially in males) are striking (Lühken et al., 2017): collecting USUV dead common blackbirds would be easier in such areas. We therefore could not exclude that the association between human high-density areas and USUV places could be explained by a sampling bias, although this hypothesis is difficult to confirm with our opportunistic data. Actually, as for wetlands, a selection for high human density areas by mosquitoes might be a part of the explanation for this association with the USUV cases. Indeed, *Cx. pipiens* can also find favourable breeding grounds in urban (Haba & McBride, 2022) or peri-urban (Vogels et al., 2016) environments (e.g. rainwater collection containers, ponds in gardens, flooded cellars or construction sites (Becker et al., 2010)).

In our study, data from the SAGIR surveillance network were used to determine USUV pressure areas on common blackbird population sizes, which have been estimated using GAMM in each class of USUV pressure. The combination of data from these two networks helped assessing the impact of USUV on common blackbird populations. The negative impact of USUV on common blackbird populations that we highlighted is consistent with what has been showed in Germany, where a difference of about 15.7% has been highlighted between the means of populations indices in 2016 in common blackbird populations (lower population index compared to the baseline year) and statistically significant lower population index in the USUV-suitable areas compared with the USUV-unsuitable areas (Lühken et al., 2017). Even if in our study the cut in three classes of USUV pressure helped us to highlight gradual effects on common blackbird populations tendencies (the greater was the infection pressure, the greater was the negative effect on common blackbird populations), our results and the German one converge in the conclusion that common blackbird populations are affected by the USUV circulation. This result points out the impact of USUV circulation on the balance of ecosystems including common blackbirds and the services they provide (e.g. a role in seed dispersal) (Whelan et al., 2008).

The presence of USUV infection in zoological gardens also raised the question of conservation of foreign species. As previously said, SAGIR network is a precious tool for early detection. But in case of endemicity, control measures must also be considered (e.g. elimination of artificial container habitats, reduction of shade sources over aquatic habitats (Tuten, 2011), use of biological larvicides such as *Bacillus thuringiensis* subsp. *israelensis* (Virgillito et al., 2022), indoor caging during periods of high density of mosquitoes, mosquito nets). For the first time in 2018, the AFVPZ and SAGIR networks were really coordinated to face the USUV outbreak. Data collected and related to both captive and non-captive wild birds allowed the description and the characterization of the outbreak within mainland France (time of onset and duration). Non-captive wild common blackbirds could be seen in this context as USUV sentinel for captive wild birds and would be of great importance for species with conservation issues. This role could also be extended to non-captive other sensitive species (e.g. capercaillie or eagle-owl), and more specifically Passeriformes, as the birds collected in rescue centres have highlighted.

Another issue of the diffusion of USUV in France in wild birds as common blackbirds, living close to human populations, and in zoological species is contact with humans, who are known as USUV incident hosts (the asymptomatic form being predominant and fever illness, or neurological symptoms as encephalitis and meningitis (Vilibic-Cavlek et al., 2014) being rare (Angeloni et al., 2023)). Zoological parks are particularly favourable to the circulation of mosquitoes which can feed on both humans and animals (Martínez-de la Puente et al., 2020). In addition, screening of asymptomatic blood donors and seroprevalence studies have also revealed incidental findings of the virus in asymptomatic human populations (Cadar et al., 2017; Bakonyi et al., 2017; Zaaijer et al., 2019; Angeloni et al., 2023). Their link to epizootic outbreaks has raised the question of the
need to strengthen USUV detection in blood donors, particularly during such periods (Cadar et al., 2017). Long-term epidemiological studies would be necessary to provide a better understanding of viral dynamics in potential wild reservoirs (Vittecoq et al., 2013). The USUV antigen detection in feather (follicle shafts and bulbs) of infected birds raised the possibility for live bird testing (Giglia et al., 2021). Risk mapping of competent vectors for USUV (Cx. pipiens and Ae. albopictus) could also help developing appropriate prevention and control measures (Martinet et al., 2023).

In a context of climate change (Walther et al., 2002; Parmesan & Yohe, 2003), the impact on arboviruses is complex to define (Franklinos et al., 2019), but the perspective of a shorter wintering period and a longer period of activity cannot be ruled out. Such changes would not be without effect on the distribution patterns of arboviruses such as USUV and recent research has shown the role of climate change in the increase of the risk associated with West Nile Virus circulation through Europe (Erazo et al., 2024). Moreover, the analysis of the origins and the within-France spatial distribution of common blackbirds have highlighted that i) they mainly came from Western and Central Europe and ii) were strongly segregated in autumn and winter in France (Lahournat et al., 2021). The coming seasonal disturbances in a context of climate change therefore raise the question of the impact of USUV in France on its epidemiology on a European scale. The contributions of networks such as SAGIR and REZOP are therefore more necessary than ever.

Following this analysis of surveillance data, one could question the potential impact of the reactivity of a network such as SAGIR, coordinated with the other networks mentioned above, and what strategy should be adopted in the event of a future outbreak. Indeed, we highlighted here the relevance of combining data collected by different networks primarily dedicated to different purposes (e.g., SAGIR and REZOP). For example, an automatization of alerts could enhance the reactivity of the entire system of networks. The results of spatial analysis could be quickly accessible and help the definition of areas of risk, in near real time, and where surveillance (and if possible, control) efforts could be concentrated. Such an automatization could also be interesting in a financial perspective: as reported above, the number of cases had to be limited to one per department after August 31st, because of too high analysis costs associated with the growing outbreak. This drawback could be circumvented with better reactivity and if necessary, control of the virus diffusion. However, additional data (i.e. additional funding to analyse more samples even during the peak of the outbreak) would make it possible to reconstruct the spread of the virus over time by following different lineages more closely.

## Acknowledgements

We would like to thank all the observers and collectors of the SAGIR and REZOP networks, the local services of the French Agency for Biodiversity (OFB), the hunters and their departmental and national federations, the departmental laboratories for analysis and the national reference laboratory. We also thank the rescue centres and the zoological parks of the AFVPZ. We finally warmly thank Lorette Hivert and Nicolas Toulet for the Epifaune database managing and Jean-Philippe Martinet (Université de Reims Champagne-Ardenne, Faculté de Pharmacie, UR7510 ESCAPE-USC ANSES PETARD, Reims, France) for his valuable help in discussing entomological aspects and his review of these points.

Preprint version 3 of this article has been peer-reviewed and recommended by Peer Community In Infections (https://doi.org/10.24072/pci.infections.100236; Liégeois, 2024).

## Funding

The Ministry of the Environment financed operators of the SAGIR (OFB operators) and REZOP networks. Collection, transport of carcasses, necropsies, transfer of samples and additional tests other than USUV were paid for by the network - either OFB or departmental federations of hunters with a participation of the Ministry of Agriculture. Zoological parks paid for the transport of their samples for analysis. The Anses financed USUV analyses. The National Reference Laboratory was financed by the DGAL and the Ministry of Agriculture.
**Conflict of interest disclosure**

The authors declare that they comply with the PCI rule of having no financial conflict of interest in relation to the content of the article.

**Data, scripts, code, and supplementary information availability**

All the data, script, code and supplementary information have been packaged in an R package named usutuFrance, available on Github at: https://github.com/ClementCalenge/usutuFrance. We have also stored this package on Zenodo (https://doi.org/10.5281/zenodo.10992191; Calenge et al., 2024). The raw dataset used in this paper has also been stored as a text file on Zenodo (https://doi.org/10.5281/zenodo.10992555; Bouchez-Zacria et al., 2024).

**References**

Agliani G, Giglia G, de Bruin E, van Mastrigt T, Blom R, Sikkema RS, Kik M, Koopmans MPG, Gröne A, Van den Brand JMA (2023) The pathology of co-infection with Usutu virus and *Plasmodium* spp. in naturally infected Eurasian blackbirds (*Turdus merula*). *One Health*, **16**, 100534. https://doi.org/10.1016/j.onehlt.2023.100534

Angeloni G, Bertola M, Lazzaro E, Morini M, Masi G, Sinigaglia A, Trevisan M, Gossner CM, Haussig JM, Bakonyi T, Capelli G, Barzon L (2023) Epidemiology, surveillance and diagnosis of Usutu virus infection in the EU/EEA, 2012 to 2021. *Eurosurveillance*, **28**, 2200929. https://doi.org/10.2807/1560-7917.ES.2023.28.33.2200929

Ashraf U, Ye J, Ruan X, Wan S, Zhu B, Cao S (2015) Usutu virus: an emerging flavivirus in Europe. *Viruses*, **7**, 219–238. https://doi.org/10.3390/v7010219

Baddeley A, Diggle PJ, Hardegen A, Lawrence T, Milne RK, Nair G (2014) On tests of spatial pattern based on simulation envelopes. *Ecological Monographs*, **84**, 477–489. https://doi.org/10.1890/13-2042.1

Bakonyi T, Erdélyi K, Ursu K, Ferenczi E, Csörgo T, Lussy H, Chvala S, Bukovsky C, Meister T, Weissenböck H, Nowotny N (2007) Emergence of Usutu virus in Hungary. *Journal of Clinical Microbiology*, **45**, 3870–3874. https://doi.org/10.1128/JCM.01390-07

Bakonyi T, Gould EA, Kolodziejek J, Weissenböck H, Nowotny N (2004) Complete genome analysis and molecular characterization of Usutu virus that emerged in Austria in 2001: Comparison with the South African Strain SAAR-1776 and other flaviviruses. *Virology*, **328**, 301–310. https://doi.org/10.1016/j.virol.2004.08.005

Bakonyi T, Jungbauer C, Aberle SW, Kolodziejek J, Dimmel K, Stiasny K, Allerberger F, Nowotny N (2017) Usutu virus infections among blood donors, Austria, July and August 2017 - Raising awareness for diagnostic challenges. *Euro Surveillance: Bulletin European Sur Les Maladies Transmissibles* = *European Communicable Disease Bulletin*, **22**, 17–00644. https://doi.org/10.2807/1560-7917.ES.2017.22.41.17-00644

Barbic L, Vilibic-Cavlek T, Listes E, Stevanovic V, Gjenero-Margan I, Ljubin-Sternak S, Pem-Novosel I, Listes I, Mlinaric-Galinovic G, Di Gennaro A, Savini G (2013) Demonstration of Usutu virus antibodies in horses, Croatia. *Vector Borne and Zoonotic Diseases (Larchmont, N.Y.)*, **13**, 772–774. https://doi.org/10.1089/vbz.2012.1236

Bażanów B, Jansen van Vuren P, Szymański P, Stygar D, Frącka A, Twardoń J, Kozdrowski R, Pawęska JT (2018) A Survey on West Nile and Usutu Viruses in Horses and Birds in Poland. *Viruses*, **10**, 87. https://doi.org/10.3390/v10020087

Beck C, Jimenez-Clavero MA, Leblond A, Durand B, Nowotny N, Leparc-Goffart I, Zientara S, Jourdain E, Lecollinet S (2013) Flaviviruses in Europe: Complex Circulation Patterns and Their Consequences for the Diagnosis and Control of West Nile Disease. *International Journal of Environmental Research and Public Health*, **10**, 6049–6083. https://doi.org/10.3390/ijerph10116049
Becker N, Jöst H, Ziegler U, Eiden M, Höper D, Emmerich P, Fichet-Calvet E, Ehichioya DU, Czajka C, Gabriel M, Hoffmann B, Beer M, Tenner-Racz K, Racz P, Günther S, Wink M, Bosch S, Konrad A, Pfeffer M, Groschup MH, Schmidt-Chanasit J (2012) Epizootic emergence of Usutu virus in wild and captive birds in Germany. *PloS One*, **7**, e32604.
https://doi.org/10.1371/journal.pone.0032604

Becker N, Petrić D, Zgomba M, Boase C, Madon M, Dahl C, Kaiser A (2010) *Mosquitoes and their control*. Springer Science & Business Media, Heidelberg, Dordrecht, New York.
https://doi.org/10.1007/978-3-540-92874-4

Ben Hassine T, De Massis F, Calistri P, Savini G, BelHaj Mohamed B, Ranen A, Di Gennaro A, Sghaier S, Hammami S (2014) First detection of co-circulation of West Nile and Usutu viruses in equids in the south-west of Tunisia. *Transboundary and Emerging Diseases*, **61**, 385–389.
https://doi.org/10.1111/tbed.12259

Bouchez-Zacria M, Calenge C, Villers A, Lecollinet S, Gonzalez G, Quintard B, Leclerc A, Baurier F, Paty M-C, Martinet J-P, Faure É, Eraud C, Decors A (2024) Data used by Bouchez-Zacria et al. in "Analysis of the Usutu episode of summer 2018 in birds in France.". *Zenodo*
https://doi.org/10.5281/zenodo.10992556

Bournez L, Umhang G, Faure E, Boucher J-M, Boué F, Jourdain E, Sarasa M, Llorente F, Jiménez-Clavero MA, Moutailler S, Lacour SA, Lecollinet S, Beck C (2019) Exposure of Wild Ungulates to the Usutu and Tick-Borne Encephalitis Viruses in France in 2009–2014: Evidence of Undetected Flavivirus Circulation a Decade Ago. *Viruses*, **12**, 10.
https://doi.org/10.3390/v12010010

Boutin J, Barbier L, Roux D (2001) Suivi des effectifs nicheurs d'alaudidés, colombidés et turdidés en France: le programme ACT. In: *Actes du 25ème Colloque francophone d'ornithologie, 18 et 19 mars 2000* r Caen. S.E.O.F, pp. 53-61. Alauda, Caen.
https://apisemanticscholar.org/CorpusID:133039289

Boutin J, Roux D, Eraud C (2003) Breeding bird monitoring in France: the ACT survey. *Ornis Hungarica*, **12**, 1-2. http://zeus.nye.hu/~szept/EBCC/25-Boutin_Eraud.pdf

Busquets N, Alba A, Allepuz A, Aranda C, Ignacio Nuñez J (2008) Usutu virus sequences in *Culex pipiens* (Diptera: *Culicidae*), Spain. *Emerging Infectious Diseases*, **14**, 861–863.
https://doi.org/10.3201/eid1405.071577

Cadar D, Becker N, Campos R de M, Börstler J, Jöst H, Schmidt-Chanasit J (2014) Usutu virus in bats, Germany, 2013. *Emerging Infectious Diseases*, **20**, 1771–1773.
https://doi.org/10.3201/eid2010.140909

Cadar D, Maier P, Müller S, Kress J, Chudy M, Bialonski A, Schlaphof A, Jansen S, Jöst H, Tannich E, Runkel S, Hitzler WE, Hutschenreuter G, Wessiepe M, Schmidt-Chanasit J (2017) Blood donor screening for West Nile virus (WNV) revealed acute Usutu virus (USUV) infection, Germany, September 2016. *Euro Surveillance: Bulletin European Sur Les Maladies Transmissibles* = *European Communicable Disease Bulletin*, **22**, 30501.
https://doi.org/10.2807/1560-7917.ES.2017.22.14.30501

Calisher CH, Gould EA (2003) Taxonomy of the virus family *Flaviviridae*. *Advances in Virus Research*, **59**, 1-19. https://doi.org/10.1016/s0065-3527(03)59001-7

Calenge C, Villers A, Bouchez-Zacria M (2024) R package usutuFrance describing the data and analysis of the Usutu episode in France in summer 2018. *Zenodo*
https://doi.org/10.5281/zenodo.10992191

Cavrini F, Gaibani P, Longo G, Pierro AM, Rossini G, Bonilauri P, Gerunda GE, Benedetto FD, Pasetto A, Girardis M, Dottori M, Landini MP, Sambri V (2009) Usutu virus infection in a patient who underwent orthotropic liver transplantation, Italy, August-September 2009. *Eurosurveillance*, **14**, 19448. https://doi.org/10.2807/ese.14.50.19448-en

Chaintoutis SC, Dovas CI, Papanastassopoulou M, Gewehr S, Danis K, Beck C, Lecollinet S, Antalis V, Kalaitzopoulou S, Panagiotopoulos T, Mourelatos S, Zientara S, Papadopoulos O (2014) Evaluation of a West Nile virus surveillance and early warning system in Greece, based on domestic pigeons. *Comparative Immunology, Microbiology and Infectious Diseases*, **37**, 131–141. https://doi.org/10.1016/j.cimid.2014.01.004
Chevalier V, Marsot M, Molia S, Rasamoelina H, Rakotondravao R, Pedrono M, Lowenski S, Durand B, Lecollinet S, Beck C (2020) Serological Evidence of West Nile and Usutu Viruses Circulation in Domestic and Wild Birds in Wetlands of Mali and Madagascar in 2008. *International Journal of Environmental Research and Public Health*, **17**, 1998.
https://doi.org/10.3390/ijerph17061998

Clé M, Barthelemy J, Desmetz C, Foulongne V, Lapeyre L, Bolloré K, Tuaillon E, Erkilic N, Kalatzis V, Lecollinet S, Beck C, Pirot N, Glasson Y, Gosselet F, Alvarez Martinez MT, Van de Perre P, Salinas S, Simonin Y (2020) Study of Usutu virus neuropathogenicity in mice and human cellular models. *PLoS Neglected Tropical Diseases*, **14**, e0008223.
https://doi.org/10.1371/journal.pntd.0008223

Clé M, Beck C, Salinas S, Lecollinet S, Gutierrez S, Van de Perre P, Baldet T, Foulongne V, Simonin Y (2019) Usutu virus: A new threat? *Epidemiology and Infection*, **147**, e232.
https://doi.org/10.1017/S0950268819001213

Constant O, Bollore K, Clé M, Barthelemy J, Foulongne V, Chenet B, Gomis D, Virolle L, Gutierrez S, Desmetz C, Moares RA, Beck C, Lecollinet S, Salinas S, Simonin Y (2020) Evidence of Exposure to USUV and WNV in Zoo Animals in France. *Pathogens*, **9**, 1005.
https://doi.org/10.3390/pathogens9121005

Constant O, Gil P, Barthelemy J, Bolloré K, Foulongne V, Desmetz C, Leblond A, Desjardins I, Pradier S, Joulié A, Sandoz A, Amaral R, Boisseau M, Rakotoarivony I, Baldet T, Marie A, Frances B, Reboul Salze F, Tinto B, Van de Perre P, Salinas S, Beck C, Lecollinet S, Gutierrez S, Simonin Y (2022) One Health surveillance of West Nile and Usutu viruses: a repeated cross-sectional study exploring seroprevalence and endemicity in Southern France, 2016 to 2020. *Euro Surveillance: Bulletin Europeen Sur Les Maladies Transmissibles = European Communicable Disease Bulletin*, **27**, 2200068.
https://doi.org/10.2807/1560-7917.ES.2022.27.25.2200068

Cramp S (1988) *Handbook of the birds of Europe, the Middle East and North Africa: the birds of the Western Palearctic*. Vol.5, Tyrant flycatchers to thrushes.
https://archive.org/details/handbookofbirdso0000unse_b6k8

Csank T, Drzewnioková P, Korytár L, Major P, Gyuranecz M, Pistl J, Bakonyi T (2018) A Serosurvey of Flavivirus Infection in Horses and Birds in Slovakia. *Vector Borne and Zoonotic Diseases (Larchmont, N.Y.)*, **18**, 206–213.
https://doi.org/10.1089/vbz.2017.2216

Decors A, Desvaux S, Payne A, Hivert L, Cardoso O, Brulez M, Gauthier D, Faure E (2022) Le réseau SAGIR, la surveillance au carrefour des enjeux. *Bulletin épidémiologique Santé Animale*, **95**, 1–9.
https://mag.anses.fr/sites/default/files/PER-028_2022-05-10_Decors_Sagir_MaqF.pdf

Diagne MM, Ndione MHD, Di Paola N, Fall G, Bedekelabou AP, Sembène PM, Faye O, Zanotto PM de A, Sall AA (2019) Usutu Virus Isolated from Rodents in Senegal. *Viruses*, **11**, 181.
https://doi.org/10.3390/v11020181

Diggle P (1990) *Time series. A biostatistical introduction*. Oxford University Press, Oxford.
https://doi.org/10.1093/oso/9780198522065.001.0001

Diggle PJ (2013) *Statistical Analysis of Spatial and Spatio-Temporal Point Patterns, Third Edition*. CRC Press.
https://doi.org/10.1201/b15326

van Dijk JGB, Hoye BJ, Verhagen JH, Nolet BA, Fouchier RAM, Klaassen M (2014) Juveniles and migrants as drivers for seasonal epizootics of avian influenza virus. *Journal of Animal Ecology*, **83**, 266–275.
https://doi.org/10.1111/1365-2656.12131

Durand B, Haskouri H, Lowenski S, Vachiery N, Beck C, Lecollinet S (2016) Seroprevalence of West Nile and Usutu viruses in military working horses and dogs, Morocco, 2012: dog as an alternative WNV sentinel species? *Epidemiology and Infection*, **144**, 1857–1864.
https://doi.org/10.1017/S095026881600011X

Eiden M, Gil P, Ziegler U, Rakotoarivony I, Marie A, Frances B, L'Ambert G, Simonin Y, Foulongne V, Groschup MH, Gutierrez S, Eloit M (2018) Emergence of two Usutu virus lineages in *Culex pipiens* mosquitoes in the Camargue, France, 2015. *Infection, Genetics and Evolution: Journal*
of Molecular Epidemiology and Evolutionary Genetics in Infectious Diseases, **61**, 151–154.
https://doi.org/10.1016/j.meegid.2018.03.020

Engel D, Jöst H, Wink M, Börstler J, Bosch S, Garigliany M-M, Jöst A, Czajka C, Lühken R, Ziegler U, Groschup MH, Pfeffer M, Becker N, Cadar D, Schmidt-Chanasit J (2016) Reconstruction of the Evolutionary History and Dispersal of Usutu Virus, a Neglected Emerging Arbovirus in Europe and Africa. *mBio*, **7**, 10.1128/mbio.01938-15. https://doi.org/10.1128/mbio.01938-15

Épidémiosurveillance Santé animale (ESA) (2018) Bilan de la circulation virus Usutu en France au 27 août 2018. https://www.plateforme-esa.fr/sites/default/files/2021-11/USUTU%20bilan%20270818_v4.pdf

Erazo D, Grant L, Ghisbain G, Marini G, Colón-González FJ, Wint W, Rizzoli A, Van Bortel W, Vogels CBF, Grubaugh ND, Mengel M, Frieler K, Thiery W, Dellicour S (2024) Contribution of climate change to the spatial expansion of West Nile virus in Europe. *Nature Communications*, **15**, 1196. https://doi.org/10.1038/s41467-024-45290-3

Escribano-Romero E, Lupulović D, Merino-Ramos T, Blázquez A-B, Lazić G, Lazić S, Saiz J-C, Petrović T (2015) West Nile virus serosurveillance in pigs, wild boars, and roe deer in Serbia. *Veterinary Microbiology*, **176**, 365–369. https://doi.org/10.1016/j.vetmic.2015.02.005

Folly AJ, Lawson B, Lean FZ, McCracken F, Spiro S, John SK, Heaver JP, Seilern-Moy K, Masters N, Hernández-Triana LM, Phipps LP, Nuñez A, Fooks AR, Cunningham AA, Johnson N, McElhinney LM (2020) Detection of Usutu virus infection in wild birds in the United Kingdom, 2020. *Euro Surveillance: Bulletin Europeen Sur Les Maladies Transmissibles = European Communicable Disease Bulletin*, **25**, 2001732. https://doi.org/10.2807/1560-7917.ES.2020.25.41.2001732

Franklinos LHV, Jones KE, Redding DW, Abubakar I (2019) The effect of global change on mosquito-borne disease. *The Lancet Infectious Diseases*, **19**, e302–e312.
https://doi.org/10.1016/S1473-3099(19)30161-6

Fros JJ, Miesen P, Vogels CB, Gaibani P, Sambri V, Martina BE, Koenraadt CJ, van Rij RP, Vlak JM, Takken W, Pijlman GP (2015) Comparative Usutu and West Nile virus transmission potential by local *Culex pipiens* mosquitoes in north-western Europe. *One Health*, **1**, 31–36.
https://doi.org/10.1016/j.onehlt.2015.08.002

Gabriel E, Rowlingson BS, Diggle PJ (2013) stpp: An R Package for Plotting, Simulating and Analyzing Spatio-Temporal Point Patterns. *Journal of Statistical Software*, **53**, 1–29.
https://doi.org/10.18637/jss.v053.i02

Gaibani P, Rossini G (2017) An overview of Usutu virus. *Microbes and Infection*, **19**, 382–387.
https://doi.org/10.1016/j.micinf.2017.05.003

Garigliany M-M, Marlier D, Tenner-Racz K, Eiden M, Cassart D, Gandar F, Beer M, Schmidt-Chanasit J, Desmecht D (2014) Detection of Usutu virus in a bullfinch (*Pyrrhula pyrrhula*) and a great spotted woodpecker (*Dendrocopos major*) in north-west Europe. *Veterinary Journal* (*London, England: 1997*), **199**, 191–193. https://doi.org/10.1016/j.tvjl.2013.10.017

Giglia G, Agliani G, Munnink BBO, Sikkema RS, Mandara MT, Lepri E, Kik M, Ijzer J, Rijks JM, Fast C, Koopmans MPG, Verheije MH, Gröne A, Reusken CBEM, van den Brand JMA (2021) Pathology and Pathogenesis of Eurasian Blackbirds (*Turdus merula*) Naturally Infected with Usutu Virus. *Viruses*, **13**, 1481. https://doi.org/10.3390/v13081481

Haba Y, McBride L (2022) Origin and status of *Culex pipiens* mosquito ecotypes. *Current Biology*, **32**, R237–R246. https://doi.org/10.1016/j.cub.2022.01.062

Hamer GL, Anderson TK, Donovan DJ, Brawn JD, Krebs BL, Gardner AM, Ruiz MO, Brown WM, Kitron UD, Newman CM, Goldberg TL, Walker ED (2014) Dispersal of Adult *Culex* Mosquitoes in an Urban West Nile Virus Hotspot: A Mark-Capture Study Incorporating Stable Isotope Enrichment of Natural Larval Habitats. *PLOS Neglected Tropical Diseases*, **8**, e2768.
https://doi.org/10.1371/journal.pntd.0002768

Hubálek Z, Halouzka J, Juřicová Z, Šikutová S, Rudolf I, Honza M, Janková J, Chytil J, Marec F, Sitko J (2008) Serologic Survey of Birds for West Nile Flavivirus in Southern Moravia (Czech Republic). *Vector-Borne and Zoonotic Diseases*, **8**, 659–666.
https://doi.org/10.1089/vbz.2007.0283
Hubálek Z, Wegner E, Halouzka J, Tryjanowski P, Jerzak L, Šikutová S, Rudolf I, Kruszewicz AG, Jaworski Z, Wlodarczyk R (2008) Serologic Survey of Potential Vertebrate Hosts for West Nile Virus in Poland. *Viral Immunology*, **21**, 247–254. https://doi.org/10.1089/vim.2007.0111

Johnson N, Fernández de Marco M, Giovannini A, Ippoliti C, Danzetta ML, Svartz G, Erster O, Groschup MH, Ziegler U, Mirazimi A, Monteil V, Beck C, Gonzalez G, Lecollinet S, Attoui H, Moutailler S (2018) Emerging Mosquito-Borne Threats and the Response from European and Eastern Mediterranean Countries. *International Journal of Environmental Research and Public Health*, **15**, 2775. https://doi.org/10.3390/ijerph15122775

Jöst H, Bialonski A, Maus D, Sambri V, Eiden M, Groschup MH, Günther S, Becker N, Schmidt-Chanasit J (2011) Isolation of usutu virus in Germany. *The American Journal of Tropical Medicine and Hygiene*, **85**, 551–553. https://doi.org/10.4269/ajtmh.2011.11-0248

Knape J (2016) Decomposing trends in Swedish bird populations using generalized additive mixed models. *Journal of Applied Ecology*, **53**, 1852–1861. https://doi.org/10.1111/1365-2664.12720

Lahournat M, Jiguet F, Villers A, Eraud C, Henry P-Y (2021) Where do thrushes migrating to France come from? Within-France distribution and temporal changes over 70 years. *European Journal of Wildlife Research*, **67**, 95. https://doi.org/10.1007/s10344-021-01525-y

Lecollinet S, Blanchard Y, Manson C, Lowenski S, Laloy E, Quenault H, Touzain F, Lucas P, Eraud C, Bahuon C, Zientara S, Beck C, Decors A (2016) Dual Emergence of Usutu Virus in Common Blackbirds, Eastern France, 2015. *Emerging Infectious Diseases*, **22**, 2225. https://doi.org/10.3201/eid2212.161272

Liégeois F (2024) Usutu virus: Importance of a participatory research network to document the spatio-temporal circulation of USUV virus. *Peer Community in Infections*, 100236. https://doi.org/10.24072/pci.infections.100236

Lühken R, Jöst H, Cadar D, Thomas SM, Bosch S, Tannich E, Becker N, Ziegler U, Lachmann L, Schmidt-Chanasit J (2017) Distribution of Usutu Virus in Germany and Its Effect on Breeding Bird Populations. *Emerging Infectious Diseases*, **23**, 1994–2001. https://doi.org/10.3201/eid2312.171257

Lupulovic D, Martín-Acebes MA, Lazic S, Alonso-Padilla J, Blázquez A-B, Escribano-Romero E, Petrovic T, Saiz J-C (2011) First serological evidence of West Nile virus activity in horses in Serbia. *Vector Borne and Zoonotic Diseases* (Larchmont, N.Y.), **11**, 1303–1305. https://doi.org/10.1089/vbz.2010.0249

Manly BFJ (1991) *Randomization and Monte-Carlo Methods in Biology*. Springer US. https://doi.org/10.1007/978-1-4899-2995-2

Mannasse B, Mendelson E, Orshan L, Mor O, Shalom U, Yeger T, Lustig Y (2017) Usutu Virus RNA in Mosquitoes, Israel, 2014-2015. *Emerging Infectious Diseases*, **23**, 1699–1702. https://doi.org/10.3201/eid2310.171017

Martin II LB (2005) Trade-offs between molt and immune activity in two populations of house sparrows (*Passer domesticus*). *Canadian Journal of Zoology*, **83**, 780–787. https://doi.org/10.1139/z05-062

Martinet J-P, Bohers C, Vazeille M, Ferté H, Mousson L, Mathieu B, Depaquit J, Failloux A-B (2023) Assessing vector competence of mosquitoes from northeastern France to West Nile virus and Usutu virus. *PLOS Neglected Tropical Diseases*, **17**, e0011144. https://doi.org/10.1371/journal.pntd.0011144

Martínez-de la Puente J, Soriguer R, Senar JC, Figuerola J, Bueno-Mari R, Montalvo T (2020) Mosquitoes in an Urban Zoo: Identification of Blood Meals, Flight Distances of Engorged Females, and Avian Malaria Infections. *Frontiers in Veterinary Science*, **7**, 460. https://doi.org/10.3389/fvets.2020.00460

McIntosh B (1985) Usutu (SAAr 1776): nouvel arbovirus du groupe B. *Int Cat Arboviruses*, **3**, 1059–1060. https://scholar.google.com/scholar_lookup?journal=Int%20Catalogue%20Arboviruses&title=Usutu%20(SA%20Ar%201776),%20nouvel%20arbovirus%20du%20groupe%20B&author=BM%20McIntosh&volume=3&publication_year=1985&pages=1059-1060&
M'ghirbi Y, Mousson L, Moutailler S, Lecollinet S, Amaral R, Beck C, Aounallah H, Amara M, Chabchoub A, Rhim A, Failloux A-B, Bouattour A (2023) West Nile, Sindbis and Usutu Viruses: Evidence of Circulation in Mosquitoes and Horses in Tunisia. *Pathogens*, **12**, 360.
https://doi.org/10.3390/pathogens12030360

Millot F, Decors A, Mastain O, Quintaine T, Berny P, Vey D, Lasseur R, Bro E (2017) Field evidence of bird poisonings by imidacloprid-treated seeds: a review of incidents reported by the French SAGIR network from 1995 to 2014. *Environmental Science and Pollution Research*, **24**, 5469–5485. https://doi.org/10.1007/s11356-016-8272-y

Moreno J, Sanz JJ, Arriero E (1999) Reproductive effort and T-lymphocyte cell-mediated immunocompetence in female pied flycatchers *Ficedula hypoleuca*. *Proceedings of the Royal Society of London. Series B: Biological Sciences*, **266**, 1105–1109.
https://doi.org/10.1098/rspb.1999.0750

Moreno-rueda G (2010) Experimental test of a trade-off between moult and immune response in house sparrows *Passer domesticus*. *Journal of Evolutionary Biology*, **23**, 2229–2237.
https://doi.org/10.1111/j.1420-9101.2010.02090.x

Morrison CA, Baillie SR, Clark JA, Johnston A, Leech DI, Robinson RA (2015) Flexibility in the timing of post-breeding moult in passerines in the UK. *Ibis*, **157**, 340–350.
https://doi.org/10.1111/ibi.12234

Mossel EC, Crabtree MB, Mutebi J-P, Lutwama JJ, Borland EM, Powers AM, Miller BR (2017) Arboviruses Isolated From Mosquitoes Collected in Uganda, 2008–2012. *Journal of Medical Entomology*, **54**, 1403–1409. https://doi.org/10.1093/jme/tjx120

Moutailler S, Yousfi L, Mousson L, Devillers E, Vazeille M, Vega-Rúa A, Perrin Y, Jourdain F, Chandre F, Cannet A, Chantilly S, Restrepo J, Guidez A, Dusfour I, Vieira Santos de Abreu F, Pereira dos Santos T, Jiolle D, Visser TM, Koenraadt CJM, Wongsokarijo M, Diallo M, Diallo D, Gaye A, Boyer S, Duong V, Piorkowski G, Paupy C, Lourenco de Oliveira R, de Lamballerie X, Failloux A-B (2019) A New High-Throughput Tool to Screen Mosquito-Borne Viruses in Zika Virus Endemic/Epidemic Areas. *Viruses*, **11**, 904. https://doi.org/10.3390/v11100904

Nikolay B, Diallo M, Boye CSB, Sall AA (2011) Usutu virus in Africa. *Vector Borne and Zoonotic Diseases (Larchmont, N.Y.)*, **11**, 1417–1423. https://doi.org/10.1089/vbz.2011.0631

Ochieng C, Lutomiah J, Makio A, Koka H, Chepkorir E, Yalwala S, Mutisya J, Musila L, Khamadi S, Richardson J, Bast J, Schnabel D, Wurapa E, Sang R (2013) Mosquito-borne arbovirus surveillance at selected sites in diverse ecological zones of Kenya; 2007 – 2012. *Virology Journal*, **10**, 140. https://doi.org/10.1186/1743-422X-10-140

Parmesan C, Yohe G (2003) A globally coherent fingerprint of climate change impacts across natural systems. *Nature*, **421**, 37–42. https://doi.org/10.1038/nature01286

Partecke J, Gwinner E (2007) Increased Sedentariness in European Blackbirds Following Urbanization: A Consequence of Local Adaptation? *Ecology*, **88**, 882–890.
https://doi.org/10.1890/06-1105

Pecorari M, Longo G, Gennari W, Grottola A, Sabbatini A, Tagliazucchi S, Savini G, Monaco F, Simone M, Lelli R, Rumpianesi F (2009) First human case of Usutu virus neuroinvasive infection, Italy, August-September 2009. *Euro Surveillance: Bulletin European Sur Les Maladies Transmissibles = European Communicable Disease Bulletin*, **14**, 19446.
https://doi.org/10.2807/ese.14.50.19446-en

Pedersen EJ, Miller DL, Simpson GL, Ross N (2019) Hierarchical generalized additive models in ecology: an introduction with mgcv. *PeerJ*, **7**, e6876. https://doi.org/10.7717/peerj.6876

Rijks JM, Kik ML, Slaterus R, Foppen R, Stroo A, IJzer J, Stahl J, Gröne A, Koopmans M, van der Jeugd HP, Reusken C (2016) Widespread Usutu virus outbreak in birds in the Netherlands, 2016. *Euro Surveillance: Bulletin European Sur Les Maladies Transmissibles = European Communicable Disease Bulletin*, **21**, 30391. https://doi.org/10.2807/1560-7917.ES.2016.21.45.30391

Rizzoli A, Bolzoni L, Chadwick EA, Capelli G, Montarsi F, Grisenti M, de la Puente JM, Muñoz J, Figuerola J, Soriguer R, Anfora G, Di Luca M, Rosà R (2015) Understanding West Nile virus
ecology in Europe: *Culex pipiens* host feeding preference in a hotspot of virus emergence.
*Parasites & Vectors*, **8**, 213. https://doi.org/10.1186/s13071-015-0831-4

Romeo C, Lecollinet S, Caballero J, Isla J, Luzzago C, Ferrari N, García-Bocanegra I (2018) Are tree squirrels involved in the circulation of flaviviruses in Italy? *Transboundary and Emerging Diseases*, **65**, 1372–1376. https://doi.org/10.1111/tbed.12874

Sanz JJ, Moreno J, Merino S, Tomás G (2004) A trade-off between two resource-demanding functions: post-nuptial moult and immunity during reproduction in male pied flycatchers. *Journal of Animal Ecology*, **73**, 441–447. https://doi.org/10.1111/j.0021-8790.2004.00815.x

Sauvage A (2016) *Merle noir Turdus merula*, in LPO *Les oiseaux de Champagne-Ardenne*. *Nidification, migration, hivernage*. Ouvrage collectif des ornithologues champardennais. Delachaux et Niestlé, Paris, pp: 360-362. https://www.delachauxetniestle.com/livre/les-oiseaux-de-champagne-ardenne/9782603020982

Simonin Y, Sillam O, Carles MJ, Gutierrez S, Gil P, Constant O, Martin MF, Girard G, Van de Perre P, Salinas S, Leparc-Goffart I, Foulongne V (2018) Human Usutu Virus Infection with Atypical Neurologic Presentation, Montpellier, France, 2016. *Emerging Infectious Diseases*, **24**, 875–878. https://doi.org/10.3201/eid2405.171122

Snow DW (1958) The Breeding of the Blackbird *Turdus merula* at Oxford. *Ibis*, **100**, 1-30.
https://doi.org/10.1111/j.1474-919X.1958.tb00362.x

Snow DW (1969) The moult of British thrushes and chats. *Bird Study*, **16**, 115-129.
https://doi.org/10.1080/00063656909476230

Steinmetz HW, Bakonyi T, Weissenböck H, Hatt J-M, Eulenberger U, Robert N, Hoop R, Nowotny N (2011) Emergence and establishment of Usutu virus infection in wild and captive avian species in and around Zurich, Switzerland--genomic and pathologic comparison to other central European outbreaks. *Veterinary Microbiology*, **148**, 207-212.
https://doi.org/10.1016/j.vetmic.2010.09.018

Thieurmel B, Elmarhraoui A (2022) suncalc: Compute Sun Position, Sunlight Phases, Moon Position and Lunar Phase. R package version 0.5.1.
https://doi.org/10.32614/cran/package.suncalc

Toulotte F, Dupuy J, Eraud C, Sallé L (2022) *Merle noir Turdus merula*, in Dupuy J. & Sallé, L. 2022, *Atlas des oiseaux migrateurs de France*. LPO, Rochefort; Biotope editions, Mèze; Muséum national d'histoire naturelle, Paris. pp 897-900.
https://sciencepress.mnhn.fr/fr/collections/inventaires-biodiversite/atlas-des-oiseaux-migrateurs-de-france-0

Tuten HC (2011) Habitat Characteristics of Larval Mosquitoes in Zoos of South Carolina, USA. *Journal of the American Mosquito Control Association*, **27**, 111–119. https://doi.org/10.2987/10-6061.1

Vilibic-Cavlek T, Kaic B, Barbic L, Pem-Novosel I, Slavic-Vrzic V, Lesnikar V, Kurecic-Filipovic S, Babic-Erceg A, Listes E, Stevanovic V, Gjenero-Margan I, Savini G (2014) First evidence of simultaneous occurrence of West Nile virus and Usutu virus neuroinvasive disease in humans in Croatia during the 2013 outbreak. *Infection*, **42**, 689–695. https://doi.org/10.1007/s15010-014-0625-1

Vilibic-Cavlek T, Petrovic T, Savic V, Barbic L, Tabain I, Stevanovic V, Klobucar A, Mrzljak A, Ilic M, Bogdanic M, Benvin I, Santini M, Capak K, Monaco F, Listes E, Savini G (2020) Epidemiology of Usutu Virus: The European Scenario. *Pathogens* (Basel, Switzerland), **9**, 699.
https://doi.org/10.3390/pathogens9090699

Villers A, Prudhomme O, Eraud C (2021) *Suivi des populations nicheuses (1996-2019) et hivernantes (2000-2019)*. Réseau national d'observation "Oiseaux de passage", OFB-FNC-FDC, Rapport interne OFB, version 2021-07-06, 22 pages.

Virgillito C, Manica M, Marini G, Rosà R, Della Torre A, Martini S, Drago A, Baseggio A, Caputo B (2022) Evaluation of *Bacillus thuringiensis* Subsp. *Israelensis* and *Bacillus sphaericus* Combination Against *Culex pipiens* in Highly Vegetated Ditches. *Journal of the American Mosquito Control Association*, **38**, 40–45. https://doi.org/10.2987/21-7024
Vittecoq M, Lecollinet S, Jourdain E, Thomas F, Blanchon T, Arnal A, Lowenski S, Gauthier-Clerc M (2013) Recent circulation of West Nile virus and potentially other closely related flaviviruses in Southern France. *Vector Borne and Zoonotic Diseases* (Larchmont, N.Y.), **13**, 610–613.
https://doi.org/10.1089/vbz.2012.1166

Vogels CBF, Möhlmann TWR, Melsen D, Favia G, Wennergren U, Koenraadt CJM (2016) Latitudinal Diversity of *Culex pipiens* Biotypes and Hybrids in Farm, Peri-Urban, and Wetland Habitats in Europe. *PLOS ONE*, **11**, e0166959. https://doi.org/10.1371/journal.pone.0166959

Walther G-R, Post E, Convey P, Menzel A, Parmesan C, Beebee TJC, Fromentin J-M, Hoegh-Guldberg O, Bairlein F (2002) Ecological responses to recent climate change. *Nature*, **416**, 389–395. https://doi.org/10.1038/416389a

Wand MC, Jones M (1995) *Kernel Smoothing*. Chapman and Hall/CRC, New York.
https://doi.org/10.1201/b14876

Weissenböck H, Bakonyi T, Rossi G, Mani P, Nowotny N (2013) Usutu virus, Italy, 1996. *Emerging Infectious Diseases*, **19**, 274–277. https://doi.org/10.3201/eid1902.121191

Weissenböck H, Kolodziejek J, Url A, Lussy H, Rebel-Bauder B, Nowotny N (2002) Emergence of Usutu virus, an African mosquito-borne flavivirus of the Japanese encephalitis virus group, central Europe. *Emerging Infectious Diseases*, **8**, 652–656.
https://doi.org/10.3201/eid0807.020094

Whelan CJ, Wenny DG, Marquis RJ (2008) Ecosystem Services Provided by Birds. *Annals of the New York Academy of Sciences*, **1134**, 25–60. https://doi.org/10.1196/annals.1439.003

Wood S (2006) *Generalized Additive Models: An Introduction With R*. CRC Press, Boca Raton, FL, USA. https://www.maths.ed.ac.uk/~swood34/igam/

Wood SN (2017) *Generalized Additive Models: An Introduction with R*. Chapman & Hall/CRC, New York. https://doi.org/10.1201/9781315370279

Zaaijer HL, Slot E, Molier M, Reusken CBEM, Koppelman MHGM (2019) Usutu virus infection in Dutch blood donors. *Transfusion*, **59**, 2931–2937. https://doi.org/10.1111/trf.15444
# Appendices to the article: Analysis of the Usutu episode of summer 2018 in birds in France.

Malika Bouchez-Zacria, Clément Calenge, Alexandre Villers, Sylvie Lecollinet, Gaelle Gonzalez, Benoit Quintard, Antoine Leclerc, Florence Baurier, Marie-Claire Paty, Éva Faure, Cyril Éraud, Anouk Decors.

## Contents

<table><tr><td><strong>Appendix A Analysis of the SAGIR dataset</strong></td><td style="text-align: right;"><strong>3</strong></td></tr><tr><td>A.1 The datasets</td><td style="text-align: right;">3</td></tr><tr><td>A.2 Spatial analysis of the bird occurrences</td><td style="text-align: right;">5</td></tr><tr><td>A.2.1 K function and the complete spatial randomness</td><td style="text-align: right;">5</td></tr><tr><td>A.2.2 Fit of a Thomas process to the spatial distribution of occurrence</td><td style="text-align: right;">7</td></tr><tr><td>A.2.3 Effect of environmental variables on the probability of a cluster</td><td style="text-align: right;">13</td></tr><tr><td>A.3 Temporal and Spatio-temporal distribution of bird occurrences</td><td style="text-align: right;">16</td></tr><tr><td>A.4 Synthesis figure</td><td style="text-align: right;">22</td></tr><tr><td><strong>Appendix B Assessment of the population trends as a function of the estimated level of infection</strong></td><td style="text-align: right;"><strong>24</strong></td></tr><tr><td>B.1 Definition of the regions with the three levels of infection</td><td style="text-align: right;">24</td></tr><tr><td>B.2 Formatting the results of ACT survey</td><td style="text-align: right;">27</td></tr><tr><td>B.3 Model fit and prediction</td><td style="text-align: right;">28</td></tr><tr><td>B.4 Calculation of trends</td><td style="text-align: right;">30</td></tr><tr><td>B.5 Summary plot</td><td style="text-align: right;">31</td></tr></table>
Introduction

This vignette is the supplementary material of the article of Bouchez-Zacria et al. (prep). The aim of
this paper is to analyze the space-time dynamics of the Usutu infection in French bird populations during
the summer 2018, and to relate it with two environmental variables, and with the resulting changes in
population size of blackbirds. A companion package named usutuFrance contains the data and R code
used for this paper, and is required to reproduce the calculations in this document. The present document
is also available as a vignette of this package. To install this package, first install the package devtools
and use the function install.github to install usutuFrance:

## If devtools is not yet installed, type
install/packages("devtools")

## Install the package usutuFrance
devtools::install_github("ClementCalenge/usutuFrance", ref="main")

Remark: on Windows, it is required to also install the Rtools (https://cran.r-project.org/bin/windows/Rtools/) on your computer to have a working devtools package (see https://www.r-project.org/nosvn/pandoc/devtools.html).

Throughout this vignette, we suppose that the reader is familiar with the model and simulations devel-
oped in the main paper. We describe the R code used for this analysis, as well as additional checks that
helped to better understand the dataset and underlying infection process.
# Appendix A Analysis of the SAGIR dataset

## A.1 The datasets

In this appendix, we provide the R code used to analyse the infection process using the data collected by the SAGIR network. We first load the package:

```r
library(usutuFrance)
```

```r
## Loading required package: mgcv
```

```r
## Loading required package: nlme
```

```r
## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.
```

We then load the data collected by the SAGIR network:

```r
data(usutu)
```

```r
str(usutu)
```

```r
## Classes 'tbl DF', 'tbl' and 'data.frame': 60 obs. of 5 variables:
```

```r
## $ species: chr "Turdus merula" "Strix nebulosa" "Turdus merula" "Turdus merula" ...
```

```r
## $ date : Date, format: "2018-07-22" ...
```

```r
## $ x : num 507007 599352 520995 579927 531686 ...
```

```r
## $ y : num 2188525 2428094 2144710 2317124 2274347 ...
```

```r
## $ result : chr "Positive" "Positive" "Positive" "Positive" ...
```

This data.frame contains, for each one of the 60 dead birds collected by the network and tested for USUV:

* The latin name of the bird species (13 bird species);

* The date of collection;

* The x and y coordinates (Lambert II projection) of the bird;

* The result of the test (either "Positive" or "Negative").

In the paper, we focused only on the birds with a positive test result. Birds with a negative test result do not bring any useful information in this context: these are not dead birds randomly sampled over France, but rather birds for which local observers had strong suspicion that their death was caused by USUV:

```r
usutup <- usutu[usuturesult=="Positive",]
```

We can show a map of these bird locations over a map of France. We load the dataset sfFrance, which is an object of class sfc containing the coordinates of the contour of mainland France; we then plot the locations on this map:

```r
## Required libraries
library.sf)
```

```r
## Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2()
```

```r
## is TRUE
```

```r
library(ggplot2)
```
```r
library(ggspatial)

data(cfFrance)

theme_set theme_bw())

ggplot()+
  geom_sf(data=sfFrance)+
  geom_point(aes(x=x,y=y),data=usutup,size=3)+
  theme( panel.grid.major = element_line(color = gray(0.5), linetype = "dashed",
                                           linewidth = 0.5),
        panel-background = elementRect( fill = "aliceblue") )+
  theme( axis line=elementDA blank(),
        axis.text.x=elementDA blank(),
        axis.text.y=elementDA blank(),
        axis.ticks=elementDA blank(),
        axis.title.x=elementDA blank(),
        axis.title.y=elementDA blank(),
        legend(position="none")+
  annotation_scale(location = "bl", width_hint = 0.5, plot_unit="m")
```
<figure><img src="image_6.png" /></figure>

We first proceed with the spatial analysis of this point pattern, then we consider the time distribution of these cases, and finally we analyse the space-time distribution of cases.

## A.2 Spatial analysis of the bird occurrences

### A.2.1 K function and the complete spatial randomness

As noted in the paper, we first calculate the L function to identify scales of clustering in the spatial distribution of dead birds (Baddeley et al., 2016, p. 203 and 207). We do not account for the species of the birds in this analysis:

```r
## Distances for which the K function is desired (from 1 to 250 km)
vs <- seq(1000, 250000, length=100)
```

```r
## The boundary of the study area
oo <- sfFrance[[1]][[1]]
```

```r
## Coordinates of the points
```
```r
xy <- as.matrix(unlist([,c("x", "y")]))
```

```r
## Calculation of the K function (divide all coordinates and distance
## by 1000: representation in km is clearer than representation in
## meters)
```

```r
kh <- splancs::khat(xy/1000, oo/1000, vs/1000)
```

```r
## Calculation function L
```

```r
kh <- sqrt(kh/pi)-vs/1000
```

We then simulate the complete spatial randomness (Baddeley et al., 2016, p. 132) 999 times to derive the envelope of the set of K functions expected under this process. WARNING: THIS CALCULATION TAKES A VERY LONG TIME !!! Note that we have included the results of this calculation as a dataset of the package, so that the reader does not need to launch this function to reproduce further calculations:

```r
## Calculation envelope by simulation of complete spatial randomness
## (min and max for different values of distances)
```

```r
simKCSR <- splancs::Kenv.csr(nrow(xy), oo/1000, 999, vs/1000, quiet=TRUE)
```

We load the results of these simulations:

```r
data(simKCSR)
```

And we plot the L-function with le limits of the envelope under complete spatial randomness:

```r
## transformation envelope of K into envelope of L
simKCSR$lower <- sqrt(simKCSR$lower/pi)-vs/1000
simKCSR$upper <- sqrt(simKCSR$upper/pi)-vs/1000
```

```r
## Display the function L with envelope expected under CSR
```

```r
plot(vs/1000, kh, ty="l",
     ylim=range(c(kh, unlist(simKCSR))), xlab="distance t (km)", ylab="L(t)")
polygon(cbind(c(vs/1000, rev vs/1000)),
        c(simKCSR$lower, rev(simKCSR$upper))), col="grey", border = NA)
abline(h=0)
```
<figure><img src="image_7.png" /></figure>

Which illustrates the clustering of the point pattern at all scales, and in particular at 5 km, 125 km and 200 km.

### A.2.2 Fit of a Thomas process to the spatial distribution of occurrence

Given the previous result, we decided to fit a Thomas process to the location data using the minimum contrast method based on the function K (Baddeley et al., 2016, p. 463). We use the functions of the package `spatstat` for this fit:

```r
## Define a point pattern (similarly, divide by 1000: it is clearer to
## work in kilometers).
```

```r
ppo <- spatstat geom::ppp(xy[,1]/1000, xy[,2]/1000,
                          window =
                            spatstat geom::owin(poly=oo[nrow(oo):1,/1000)))
```

```r
## Fit the Thomas process (no trend)
```

```r
ppro <- spatstat model::kppm(ppo, ~1, "Thomas")
```
We then simulate this process 999 times to derive an envelope for the set of K functions expected under this process. WARNING: THIS CALCULATION TAKES A VERY LONG TIME !!!! Note that we have included the results of this calculation as a dataset of the package, so that the reader does not need to launch this function to reproduce further calculations:

```r
simKraw <- sapply(1:999, function(r) {
  cat(r, "\r")
  s <- simulate(ppro, 1)[[1]]
  khs <- splancs::khat(cbind(s$x, s$y), oo/1000, vs/1000)
  khs <- sqrt(khs/pi) - vs/1000
  return(khs)
})
```

We load the results of these simulations:

```r
data(simKraw)
```

And we show the fit based on the L-function:

```r
## And show the fit with:
rk <- t(apply(simKraw, 1, function(x) range(x)))
vsb <- vs/1000

## prepare the plot
plot(vsb, vsb, ty="n", ylim=range(c(c(rk), kh)),
     ylab="Function L", xlab="Distance")

## show the envelope
polygon(c(vsb, rev(vsb)), c(rk[,1], rev(rk[,2])),
        col="grey", border=NA)

## and the function L
lines(vsb, kh)
```
<figure><img src="image_8.png" /></figure>

Note that at very small scales, the L-function shows a clustering not accounted for by the Thomas process. This corresponds to the small peak already identified at 5 km. This peaks corresponds to a tendency of observers to look for other dead birds in the same commune when a dead bird has been identified there by the observer.

*Remark:* One possible reason for this tendency to report several birds at the same place could be that the observers tend to search for other dead birds belonging to different species when they have reported a first dead bird nearby. If this is the case, then one would expect that animals reported at small distance from each other belong to different species more frequently than expected by chance. We carried out a randomization test to test this hypothesis (Manly, 1991). More precisely, we tested whether the pairs of birds reported at a distance lower than 5 km were more often of the same species than expected by chance. There are actually 10 pairs of birds located at less than 5 km from each others, among which 4 are from the same species:

### Removal of all points located at less than 5 km from another point
### Calculation of distance between points
dp <- as.matrix(dist(xy))
### identification, for eachpoint, of the index of the first point
### for which the distance is lower than 5 km (when the distance between
i,j is lower than 5 km, if i<j, then the first point on the i-th row
will be the point i, and the first point on the j-th row will also be
the point i).

firstL5 <- sapply(1:nrow(dp), function(x) min(which(dp[x,]<5000)))
## Remove duplicated points
dup <- duplicated(firstL5)

## Number of pairs:
sum(dup)

[1] 10

## Number of pairs of the same species:
w <- which(dup)
sum(utup$species[w] == utup$species[w-1])

[1] 4

## Which corresponds to two pairs blackbirds and two pairs of Great
## grey owl
utup$species[w][utup$species[w] == utup$species[w-1]]

[1] "Strix nebulosa" "Turdus merula" "Turdus merula"
[4] "Strix nebulosa"

We carried out this test by randomly permuting the species associated to the birds in our dataset and
counting the number of times pairs of the same species were located at less than 5 km. We performed
999 simulations of this process. The P-value of this test is calculated below:

set.seed(777)
ra <- sapply(1:999, function(r) {
  ## Random permutation
  w <- which(c(FALSE, sample(dup[-1])))
  ## Number of pairs of the same species
  sum(utup$species[w] == utup$species[w-1])
})

## The P-value of this test is
mean.ra<=4)

[1] 0.4454454

We therefore fail to identify any significant tendency to report less birds of the same species when they are
collected at the same place. Thus, the observed clustering at small distances is not clearly caused by a ten-
dency to report birds of different species. We decided to continue to ignore the species in further analyses.

To improve the fit of the Thomas process, we decided to thin the point pattern prior to the fit, keeping
only the first found bird when two birds were collected at a distance lower than 5 km:

xy_thin <- xy[!dup,]

We then calculated the K and L functions on this thinned point pattern:
```
thThin <- splancs::khat(xy_thin/1000, oo/1000, vs/1000)
```

## Function L

```
thThin <- sqrt(kh_thin/pi)-vs/1000
```

And we fitted the Thomas process again to this thinned dataset:

```r
## Define a point pattern for the package spatstat (divide by 1000 to work on km)
ppo_thin <- spatstat geom::ppp(xy_thin[,1]/1000, xy_thin[,2]/1000,
                               window =
                                 spatstat geom::owin(poly=oo[nrow(oo):1,]/1000))
```

## Fit the Thomas process

```
ppro_thin <- spatstat model::kppm(ppo_thin, ~1, "Thomas")
```

Again, we can simulate this Thomas process 999 times, to derive an envelope for the set of K functions expected under this process. WARNING: THIS CALCULATION TAKES A VERY LONG TIME !!!!
Note that we have included the results of this calculation as a dataset of the package, so that the reader
does not need to launch this function to reproduce further calculations:

```r
simKthin <- sapply(1:999, function(r) {
  cat(r, "\r")
  s <- simulate(ppro_thin,1)[[1]]
  khs <- splancs::khat(cbind(s$x,s$y), oo/1000, vs/1000)
  khs <- sqrt(khs/pi)-vs/1000
  return(khs)
})
```

We load the results of these simulations:

```
data(simKthin)
```

And we show the L-function together with this calculated envelope:

## And show the fit with:

```
rk <- t(apply(simKthin,1, function(x) range(x)))
vsb <- vs/1000
```

## prepare the plot

```
plot(vsb,vsb, ty="n", ylim=range(c(c(rk),kh_thin)),
    ylab="Function L", xlab="Distance")
```

## show the envelope

```
polygon(c(vsb, rev(vsb)), c(rk[,1], rev(rk[,2])),
        col="grey", border=NA)
## and the function L
lines(vsb, kh_thin)
```
<figure><img src="image_9.png" /></figure>

The fit is now correct. The estimated parameters of the models are therefore:

```R
## The log-density of clusters
co <- coef(summary(ppro_thin))
```

```R
## get the estimate and confidence limits, calculate the exponential
## to estimate the density, and multiply by 10000 to have the mean
## number of cases per 10000 squared km (mean number of cases in
## a square of 100 x 100 km)
exp(co[c(1,3,4)])*10000
```

```R
## Estimate CI95.lo CI95.hi
## (Intercept) 0.7408151 0.3895074 1.408977
```

```R
## All the other parameters of the process
ppro_thin
```

```R
## Stationary cluster point process model
## Fitted to point pattern dataset 'ppo_thin'
```
## Fitted by minimum contrast
## Summary statistic: K-function
##
## Uniform intensity: 7.408151e-05
##
## Cluster model: Thomas process
## Fitted cluster parameters:
## kappa scale
## 1.699252e-05 7.690762e+01
## Mean cluster size: 4.359655 points
##
## Cluster strength: phi = 0.7918
## Sibling probability: psib = 0.4419

Thus, there is in average 0.74 cases in a square of 100×100 km (the parameter **intensity** is the density per squared km). This density of cases is the product of a density of 0.17 clusters in a square of 100 × 100 km (parameter **kappa**) and a mean cluster size of 4.36 cases per cluster. The standard deviation of the distribution of distances between the cases and the centroid of the cluster is of 77 km (parameter **scale**).

Note that this estimated standard deviation of 77 km is consistent with the K function, which identified clusters at a scale of 130 km. Indeed, this standard deviation measures the typical distance between the cases and the centroid of the cluster, whereas the K function is based on the distance between two cases. We can check by simulation that the mean distance between points sampled for a bivariate Gaussian distribution with a standard deviation of 77 km is approximately of 130 km:

## Simulation of 1000 points drawn from a bivariate Gaussian
distribution with standard deviation equal to 77 km, and
calculation of the mean distance between pairs of points.
mean(dist(cbind(rnorm(1000,0,77),rnorm(1000,0,77))))

[1] 135.0801

### A.2.3 Effect of environmental variables on the probability of a cluster

We now load the raster maps of potential wetlands and log-human population density in France (see the help page of the dataset `usutuEnvir` for more details on how these maps were built). We then show these maps, together with a kernel smoothing of the point pattern (smoothing parameter *h*=125 km, Wand and Jones, 1995), to visualize the relationship between the density of cases and the environmental variables:

```r
library(sp)
data(usutuEnvir)

par(mfrow = c(1,2),mar=c(0,0,2,0))
image(usutuEnvir[1], col=viridis::magma(10))
plot.sfFrance, add=TRUE)
kd <- MASS::kde2d(xy[,1], xy[,2], 125000, n=400, lims=c(range(oo[,1]), range(oo[,2])))
contour(kd$x, kd$y, kd$z, add=TRUE, lwd=2, col="lightgrey", nlevels=4,
        levels=seq(0,1.4e-11, length=10)[c(1,4,7,10)], drawlabels = FALSE)
title("Potential Wetlands")
image(usutuEnvir[2], col=viridis::magma(10))
plot.sfFrance, add=TRUE)
contour(kd$x, kd$y, kd$z, add=TRUE, lwd=2, col="lightgrey", nlevels=4,
        levels=seq(0,1.4e-11, length=10)[c(1,4,7,10)], drawlabels = FALSE)
```
```r
levels=seq(0,1.4e-11, length=10)[c(1,4,7,10)],drawlabels = FALSE)
title("Log-Human population density")
```

Potential Wetlands

<figure><img src="image_10.png" /></figure>

Log-Human population density

<figure><img src="image_11.png" /></figure>

There seems to be a relationship between cases and both variables: some clusters are found in high human density areas (Paris and Nantes), whereas others are located close to wetlands (Nantes again, as well as the Brenne region – which corresponds to the main cluster at the south-west of Paris).

We used a randomization approach to test whether the mean value of the environmental variables were greater in the places where the dead birds positive to USUV were found than expected by chance (Manly, 1991). More precisely, we simulated the Thomas process fitted in the previous section 999 times over France, and we calculated for each simulation the mean of the density of humid areas and the mean of the log of human population density. We then compared these observed values with the simulated distribution. We first calculated the vector of means of these variables:

```r
sppo <- sp::SpatialPoints(xy_thin, proj4string=sp::CRS(sp::proj4string(usutuEnvir)))
(obs Means <- colMeans(sp::over(sppo, usutuEnvir)))
```

## Wetlands PopDens
## 4.813901 5.017741

The probability of wetlands is an index varying from 0 (no wetlands) to 20 (wetland certain). The population density here correspond to the mean of the log-transformed population density (transformation log(x+1), where x is expressed in number of inhabitant per squared kilometers).

We then simulate the Thomas process 999 times, and calculate the mean values of the two variables for each simulation. WARNING: THIS CALCULATION TAKES A VERY LONG TIME !!!! Note that we have included the results of this calculation as a dataset of the package, so that the reader does not need to launch this function to reproduce further calculations:

```r
simEnv <- t(sapply(1:999, function(r) {
  cat(r, "\r")
  s <- simulate(ppro_thin, 1)[[1]]
  sppob <- sp::SpatialPoints(cbind(s$x, s$y)*1000,
                             proj4string=sp::CRS(sp::proj4string(usutuEnvir)))
  return(colMeans(sp::over(sppob, usutuEnvir), na.rm=TRUE))
}))
```
We load the result of these simulations:

data(simEnv)

We can compare the observed means for the two variables with the simulated mean and standard error for these variables, derived from the simulations:

```r
## Mean and standard error for the index of probability of wetlands
c(mean(simEnv[,1]), sd(simEnv[,1]))
```

```r
## [1] 2.6915701 0.6507697
```

```r
## Recall the observed mean:
obs Means[1]
```

```r
## Wetlands
## 4.813901
```

```r
## Mean and standard error for the log-human population density
c(mean(simEnv[,2]), sd(simEnv[,2]))
```

```r
## [1] 3.5023693 0.2507175
```

```r
## Recall the observed mean:
obs Means[2]
```

```r
## PopDens
## 5.017741
```

Clearly, the areas where many dead birds positive to USUV are found are characterized by a high probability of wetlands and a high density of human population.

We can display the results of this bivariate test graphically, using the function `biv.test` from the package `adehabitatHS`:

```r
adehabitatHS::biv.test(as.data.frame(simEnv), obs Means)
```

```r
## Registered S3 methods overwritten by 'adehabitatMA':
```

```r
## method             from
```

```r
## print.SpatialPixelsDataFrame sp
```

```r
## print.SpatialPixels sp
```
<figure><img src="image_12.png" /></figure>

The abscissa corresponds to the probability of wetlands, and the ordinate corresponds to the log-human population density. Each small black point corresponds to the vector of means of environmental variables for one simulation of the Thomas process, whereas the grey point corresponds to the observed value for this vector of means. The margins of this graph show the statistical distribution of these means expected by chance, as well as the observed value for each variable. These margins also show the P-values of each univariate test (which are presented in the main text). This test confirms that USUV cases are more frequent in places where both human population and wetlands are higher.

## A.3 Temporal and Spatio-temporal distribution of bird occurrences

We now present how the mean number of cases per day change throughout the study period. We smoothed this curve using moving average approach (Diggle, 1990, p. 22): we present, for each day *d* of the study period, the mean number of cases per day calculated over a week centred on *d*:

```r
library(lubridate)
##
## Attaching package: 'lubridate'
```
```r
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union

seqPeriod <- seq(ymd('2018-07-15'), ymd('2018-08-31'), by = '1 day')
movAvg <- sapply(1:length(seqPeriod), function(i) {
  x <- seqPeriod[i]
  ## days of the week centred on the day i
  d <- seq(x-3, x+3, by="1 day")
  ## number of days of the week in the study period
  ndays <- mean(d>=seqPeriod[1]&d<=seqPeriod[length(seqPeriod)])
  ## number of cases reported in this period
  ncases <- sum(unutup>= (seqPeriod[i]-3)&unutup$\d <= (seqPeriod[i]+3))
  ## number of cases per day:
  return.ncases/ndays
})

da <- data.frame(date=seqPeriod,
                 NumberCases=movAvg)

ggplot(da, aes(x date, y=NumberCases))+geom_line() +xlab("Date") +
  ylab("Mean number of cases reported per day")
```
<figure><img src="image_13.png" /></figure>

A graphical analysis of the spatio-temporal patterns in USUV occurrence starts with a map of the
occurrence where the color indicates the time at which the bird was found:

```r
theme_set themeBW())

labd <- pretty (usutup $ date)

ggplot () +
  geom.sf (data = sfFrance) +
  geom_point (aes (x = x, y = y, col = as.numeric (date), group = as.factor (date)), data = usutup, size = 3) +
  theme (panel.grid.major = element_line (color = gray (0.5),
 linetype = "dashed",
 linewidth = 0.5),
 panel-background = elementRect (fill = "aliceblue")) +
  theme (axis line = elementDA blank (),
 axis.text.x = elementDA blank (),
 axis.text.y = elementDA blank (),
 axis.ticks = elementDA blank (),
 axis.title.x = elementDA blank (),
 axis.title.y = elementDA blank (),
 legend position = "none") +
```
viridis::scale_color_viridis(breaks = as.numeric(labd),
                             labels = labd)+
annotation_scale(location = "bl", width_hint = 0.5, plot_unit="m")

<figure><img src="image_14.png" /></figure>

Blue points correspond to birds found at mid-July, whereas yellow points correspond to birds reported by the end of August. It seems that the main clusters from center and western France appeared approximately at the same time at about mid-July.

We calculated the Space-Time K-function based on this point pattern (Gabriel *et al.*, 2013):

```r
## get define space-time points
po <- stpp::as.3dpoints(usutup$x/1000, usutup$y/1000,
                       as.numeric(usutup hédate= min (usutup$date))))

## The set of distances for which the function K is calculated
## (between 1 to 220 km)
vs <- seq(1000, 220000, length=100)
## The weeks
```
```r
vt <- 1:10

## the dates
dates <- as.numeric(unutup$date-min(unutup$date))

## Space-time function
ch <- stpp::STIKhat po, s.region = oo/1000,
               t.region = c(min(dates-1), max(dates+1)),
               times=vt, dist=vs/1000, infectious=TRUE)
```

We then wanted to know whether birds found at small distance from each other were also found at close dates. We carried out a set of simulations to test whether this space-time K function was greater than expected by chance, by randomizing the collection dates associated with the bird locations. WARNING: THIS CALCULATION TAKES A VERY LONG TIME !!!! Note that we have included the results of this calculation as a dataset of the package, so that the reader does not need to launch this function to reproduce further calculations:

```r
## Randomization test: WARNING! THIS PART IS VERY SLOW
si <- list()

## 999 randomization (slow calculation !!!)
for (i in 1:999) {
  cat(i, "\r")
  pos <- po
  pos <- stpp::as.3dpoints(unutup$x/1000, unutup$y/1000, sample(dates))
  chs <- stpp::STIKhat(pos, s.region = oo/1000,
                     t.region = c(min(dates-1), max(dates+1)),
                     times=vt, dist=vs/1000, infectious=TRUE)
  si[[i]] <- chs$Khat
}

## Calculation of the proportion of simulations greater than observation
KstFun <- chs$Khat
## Initialization
for (i in 1:nrow(KstFun)) {
  cat(i, "\r")
  for (j in 1:ncol(KstFun)) {
    ot <- sapply(1:length(si), function(k) si[[k][i,j])
    KstFun[i,j] <- mean(ot>chs$Khat[i,j])
  }
}
```

We load the dataset containing the resulting matrix:

```r
data(KstFun)
```

We then show, for each distance and duration between occurrence, the probability that the observed space-time K function is greater than the set of functions simulated under the hypothesis of absence of time structure:

```r
## Show the results:
dok <- data.frame(duration=rep(vt,length(vs)),
                  distance=rep(vs/1000,each Service length(vt)),
                  Prob. =asvector(t(KstFun)))

ggplot2::ggplot(dok, ggplot2::aes(x=duration, y=distance))+
```
```r
ggplot2::geom tile(ggplot2::aes (fill=Prob.))+
viridis::scale fill viridis (option="magma", limits=c(0,1))+
ggplot2::geom contour(ggplot2::aes (z=Prob.), breaks = 0.05, col="white",
 alpha=0.5, linewidth=1)+ ggplot2::theme bw ()+
xlab("Duration (days)") + ylab("Distance (km)")
```

<figure><img src="image_15.png" /></figure>

The contour limit identifies the set of distances/duration pairs for which the proportion of simulated K functions greater than the observed K function is lower than 5%.

There is a strong clustering for small distances (as soon as 2 or 3 days after a case) of about 10 to 20 km, expanding to 50-60 km five days after a case. But we also observe a clustering on larger distances (150 km) starting as soon as 2 or 3 days after a case, which increases to reach a large range of distances after a week (90 to 150 km).
A.4 Synthesis figure

We finally draw the synthesis figure that is displayed in the main paper (Fig. 1), and which summarizes
the main results of this analysis:

```r
##png(filename="SynthesisFigure.png", width=1000, height=1000, pointsize=20)
library(gridBase)
library(grid)

par(mfrow=c(2, 2))
plot(vs/1000, kh, ty="l",
     ylim=range(c(kh,unlist(simKCSR))), xlab="Distance t (km)", ylab="L(t)",
     main="(A)")
polygon(cbind(c(vs/1000, rev vs/1000)),
        c(simKCSR$lower, rev(simKCSR$upper))), col="grey", border = NA)
abline(h=0)

opar <- par(mar=c(1,1,4,1))
kd <- MASS::kde2d(xy[,1], xy[,2], 125000, n=400, lims=c(range(oo[,1]), range(oo[,2])))
plot.sfFrance, main="(B)", col=grey(0.8), border=NA)
contour(kd$x, kd$y, kd$z, add=TRUE, lwd=2, col=grey(0.6), nlevels=4,
        levels=seq(0,1.4e-11, length=10)[c(1,4,7,10)], drawlabels = FALSE)
usud <- unclass(usutup$x$date-min usutup$x$date)
points usutup[,c("x","y")], pch=3, col="black", cex=1.5)

da <- data.frame(date=seqPeriod,
                 NumberCases=movAvg)
par(opar)
plot(da$xdate, da$NumberCases, ty="l", xlab = "Date",
     ylab="Mean number of cases reported per day", main="(C)")
plot.new()
vps <- baseViewports()
pushViewport(vps$figure)
vp1 <- plotViewport(c(1.8,1,3,1))
p <- ggplot2::ggplot(dok, ggplot2::aes(x=duration, y=distance))+
    ggplot2::geom tile(ggplot2::aes(fill=Prob.))+
    viridis::scale_fill_viridis option="magma", limits=c(0,1))+
    ggplot2::geom contour(ggplot2::aes(z=Prob.), breaks = 0.05, col="white",
                            alpha=0.5, linewidth=0.8)+ggplot2::theme bw()+
    theme(axis.text=element_text(size=12),
          axis.title=element_text(size=13))+
    xlab("Duration (days)") + ylab("Distance (km)") + title("(D)")
print(p, vp = vp1)
##dev.off()
```
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>
# Appendix B Assessment of the population trends as a function of the estimated level of infection

In this appendix, we provide the R code used to correlate the level of infection by the Usutu virus in a region with the changes in blackbird population size in this region.

## B.1 Definition of the regions with the three levels of infection

We will need the object `kd` created in appendix A, which stores the density of cases over continental France, estimated by the kernel smoothing method. We first convert this object into an object of class `spatRaster` from the package `terra`:

```R
## Load the package terra
library(terra)
```

```R
## terra 1.7.65
```

```R
##
```

```R
## Attaching package: 'terra'
```

```R
## The following object is masked from 'package:grid':
```

```R
##
```

```R
## depth
```

```R
## The following object is masked from 'package:knitr':
```

```R
##
```

```R
## spin
```

```R
## Convert to spatRaster
```

```R
m <- kd
```

```R
resx <- diff(kd$x) [1]
```

```R
resy <- diff(kd$y) [1]
```

```R
xmn <- min(m$x) - 0.5 * resx
```

```R
xmx <- max(m$x) + 0.5 * resx
```

```R
ymn <- min(m$y) - 0.5 * resy
```

```R
ymx <- max(m$y) + 0.5 * resy
```

```R
z <- m$z[,ncol(m$z):1]
```

```R
r1 <- rast(ext(xmn, xmx, ymn, ymx), resolution=c(resx, resy))
```

```R
r1[] <- as.numeric((z))
```

```R
## Conversion of sfFrance to spatVector
sfv <- vect.sfFrance)
```

```R
## Set the map to NA outside the limits of continental France
rr <- rasterize(sfv, r1)
```

```R
r1 <- rr*r1
```

```R
image(r1)
```
<figure><img src="image_20.png" /></figure>

Now, we load the location of the sampling routes and points in the ACT survey, stored in the dataset `resop`:

```
data(resop)
```

```
head(resop)
```

```
## id_point cd route no_point cd route point insee dep
```

```
## 1   38   1015      3   1015_3      22
```

```
## 2   39   1016      3   1016_3      22
```

```
## 3   40   1017      3   1017_3      22
```

```
## 4   41   1018      3   1018_3      56
```

```
## 5   42   1019      3   1019_3      56
```

```
## 6   43   1020      3   1020_3      56
```

```
##       x     y
```

```
## 1 255180.0 2415340
```

```
## 2 255380.0 2403005
```

```
## 3 253523.8 2378385
```

```
## 4 247617.0 2359967
```

```
## 5 252360.0 2341680
```
## 6 250360.0 2322560

plot(resop[,c("x","y")], asp=1)

<figure><img src="image_21.png" /></figure>

We then join these locations to the map of the USUV density, and store the result in the object `resop`:

```
resop$usutu<-terra::extract(r1, as.matrix(resop[,c("x", "y")]))[,1]
```

We then average the USUV density per route (calculation of the mean density over all sampling points of a route, as well as calculation of the centoid of the route):

```
val Agg <- resop |> dplyr::group_by(cd route) |>
dplyr::summarise(usutu\_agg=mean(usutu, na.rm=T),
                  x=mean(x), y=mean(y))
head(val\_agg)
```

## # A tibble: 6 x 4
```
## cd_route usutu Agg x y
## <chr> <db1> <db1> <db1>
## 1 1015 5.73e-13 255204. 2415061.
## 2 1016 1.10e-12 255307. 2402961.
## 3 1017 2.59e-12 253418. 2378671.
## 4 1018 3.30e-12 247564. 2359936.
## 5 1019 3.32e-12 252242. 2341596.
## 6 1020 2.66e-12 250260. 2322750.
```

We then calculate the median and third quartile of the density distribution. These density values will serve to define three levels of infections: low (areas with density lower than the median), medium (areas with density comprised between the median and the third quartile), and high (areas with density greater than the third quartile):

```
## Median
usutu_med <- median(val Agg $usutu Agg)

## Third quartile
usutu_75 <- quantile(val Agg $usutu Agg, PROBS = 0.75)

## new variable in val Agg, defining the three levels of infection
val Agg $usutuF <- ifelse(val Agg $usutu Agg <= usutu_med, "low",
                       ifelse(val Agg $usutu Agg > usutu_75, "high", "medium"))
val Agg $usutuF <- factor(val Agg $usutuF, LEVELS = c("low", "medium", "high"))

head(val Agg)

## # A tibble: 6 x 5
## cd route usutu agg x y usutuF
## <chr> <db1> <db1> <db1> <fct>
## 1 1015 5.73e-13 255204. 2415061. medium
## 2 1016 1.10e-12 255307. 2402961. medium
## 3 1017 2.59e-12 253418. 2378671. high
## 4 1018 3.30e-12 247564. 2359936. high
## 5 1019 3.32e-12 252242. 2341596. high
## 6 1020 2.66e-12 250260. 2322750. high
```

## B.2 Formatting the results of ACT survey

Now, we load the results of the counts of singing male blackbirds for all routes between 1996 and 2022, which are stored in the dataset merula:

```
data(merula)

head(merula)

## tot_MN x y datetime year
## 1 6 357582.0 2315954 1996-04-01 07:22:00 1996
## 2 12 605490.0 1753844 1996-04-01 07:00:00 1996
## 4 7 403730.0 2073529 1996-04-02 07:05:00 1996
## 5 1 849576.0 2115168 1996-04-02 08:00:00 1996
## 6 0 865403.2 1833320 1996-04-02 07:30:00 1996
## 7 5 938164.4 1974468 1996-04-02 07:15:00 1996
## TSR cd routeF yearF jd passageN
## 1 -20.233333 1420 1996 92 First
## 2 -32.883333 2348 1996 92 First
```
```
## 4 -34.216667   1632  1996  93 First
## 5   44.016667   3230  1996  93 First
## 6   12.183333   3244  1996  93 First
## 7    2.116667   3537  1996  93 First
```

The meaning of the columns of this data.frame is explained on the help page of merula. We join this dataset and the raster map r1 storing the density of usutu cases, and we store the results in the dataset merula. We also define a factor containing the three levels of infections (based on the limits calculated before):

```
## Join with the density map
ex2 <- terra::extract(r1, as.matrix(merula[,c("x","y")]))
```

```
## Store the results in the dataset:
merula$usutu Agg <- ex2[,1]
```

```
## median and third quartile limits:
limi <- c(USUTU_med, USUTU_75)
```

```
## Cut the density in three classes and store the result in merula:
cu <- cut(merula$usutuagg, c(-10,limi,10))
cu <- factor(cu, labels=c("low","medium","high"))
merula$usutuF <- cu
```

## B.3 Model fit and prediction

We now fit the Generalized Additive Mixed Model (GAMM) to predict the count outcome as a function of year, julian date, time since sunrise, sampling route, and the level of infection by the Usutu virus (see main text for an explanation). WARNING: THIS CALCULATION TAKES A VERY LONG TIME (2-3 hours on a computer equipped with Intel Core i7) !!!! The results of this fit, which occupy a large space (52 MB), are not included in the package. After a check of the model fit, we used the function compute_distri from our package to simulate 1000 count predictions for each year between 1996 and 2022 on an average road located in regions characterized by varying levels of Usutu virus infection – low, medium, or high. More precisely, 1000 coefficient vectors are simulated from a multivariate Gaussian distribution using GAMM results (with distribution parameters corresponding to the mean coefficient vector and the covariance matrix of coefficients extracted from the model). Each simulated coefficient vector is used to predict the response variable (total number of male blackbirds counted on an average route) for each year in each of the three regions. These 1000 simulations per region allow the assessment of prediction uncertainty. Note that we have included the results of this simulation as a dataset of the package, so that the reader does not need to execute this function to reproduce further calculations:

```
## Load the package mgcv
library(mgcv)
```

```
## Model fit with bam (large dataset)
m_per <- bam(tot_MN ~ usutuF + s(year, k=10) + s(year, by=usutuF, k=10) +
            s(cdRouteF, bs="re") + s(yearF, bs="re") +
            s(jd, k=6) + s(TsSR, k=6), family=tw(),
            data=merula, method="REML")
```

```
## Check the model fit
gratia::appraise(m_per)
```

```
## Nothing worrying there
```

```
## Summary of the model
```
summary(m_per)

## Plot the smooth functions of explanatory variables:

plot(m_per)

## A data.frame for predictions, with unknown levels for the 2 random effects factors,
## to compute changes independently of these effects

df_pred<-data.frame(expand.grid( usutuF=levels( merula$usutuF ), cd_routeF="1", year=1996:2022, jd=mean( merula$jd ), TsSR=mean( merula$TsSR ), yearF=as.factor( 2030 )))

## Load the functions to compute simulations and estimate change

simu_byRegions<-compute_distri(m_per, df_pred, time_var="year", factor_name="usutuF", nreplicates=1000)

We load the resulting dataset:

data(simu_byRegions)

str(simu_byRegions)

## List of 3

## $ low : num [1:1000, 1:27] 6.15 5.74 6.23 6.2 5.93 ...

## ..- attr(*, "dimnames")=List of 2

## ...$ : NULL

## ...$ : chr [1:27] "1996" "1997" "1998" "1999" ...

## $ medium: num [1:1000, 1:27] 7.14 6.7 7.09 7.5 7.47 ...

## ..- attr(*, "dimnames")=List of 2

## ...$ : NULL

## ...$ : chr [1:27] "1996" "1997" "1998" "1999" ...

## $ high : num [1:1000, 1:27] 7.02 6.58 6.85 6.58 6.87 ...

## ..- attr(*, "dimnames")=List of 2

## ...$ : NULL

## ...$ : chr [1:27] "1996" "1997" "1998" "1999" ...

See the help page of this dataset for more details. We transform this list into a data.frame (layout: one simulation per row, corresponding to a unique combination of simulation ID, infection level, and year):

# convert to a data frame from which to average values per relevant time periods

df_4analyses<- do. call("rbind", simu_byRegions) |> as.data.frame() |>

dplyr::mutate( usutuF=rep( levels( merula$usutuF ), each=1000 ),
               id_sim=rep( 1:1000, 3)) |>

tidyr::pivot_longer(!c( usutuF, id_sim ), names_to="year",
                     values_to="abundance") |>

dplyr::mutate_at('year', as.numeric)

head(df_4analyses)

## # A tibble: 6 x 4

## usutuF id_sim year abundance

<chr> <int> <dbl> <dbl>

## 1 low 1 1996 6.15

## 2 low 1 1997 6.17

## 3 low 1 1998 6.19

## 4 low 1 1999 6.20

## 5 low 1 2000 6.22

## 6 low 1 2001 6.24
## B.4 Calculation of trends

We now calculate, for each simulation and each infection level, the average abundance over the years preceding the usutu infection (ranging between 2015 and 2018) and the years following this infection (ranging between 2019-2022). This comparison will enable to assess the differences in mean abundance before and after the infection:

```r
## A data frame where abundance is averaged by simulation,
## by period (before after the 2018 Usutu episode) and
## by areas of Usutu "prevalence" (usutuF)
df_4trends <- df_4analyses |
dplyr::filter(year > 2014) |
dplyr::mutate usutu_BA = factor(ifelse(year < 2019, "2015-2018", "2019-2022"),
  levels = c("2015-2018", "2019-2022")) |
dplyr::group_by(id_sim, usutu_BA, usutuF) |
dplyr::summarise(mean_abund = mean(abundance))

## 'summarise()' has grouped output by 'id_sim', 'usutu_BA'.
## You can override using the '.groups' argument.

## Define the order of the levels of infection
df_4trends $usutuF <- factor(df_4trends $usutuF, levels = c("low", "medium", "high"))

#Compute the trends for the 3 areas according to the previous values
df_trends <- df_4trends |
dplyr::group_by(id_sim, usutu_BA, usutuF) |
tidyr::pivot_wider(names_from = usutu_BA, values_from = mean_abund) |
as.data.frame()

c头(df_trends)

## id_sim usutuF 2015-2018 2019-2022
## 1 1 high 6.707626 5.815715
## 2 1 low 6.540335 6.614478
## 3 1 medium 6.898128 6.519382
## 4 2 high 6.232447 5.316868
## 5 2 low 6.111560 6.195160
## 6 2 medium 6.676841 6.234477
```

We can calculate the rate of variation $100 \times (X_2 - X_1)/X_1$ for each simulation and each level of infection (with $X_1$, $X_2$ the predicted count respectively before and after the infection). Then, we average the 1000 rates of variations for each level of infection and we calculate the 95% confidence interval:

```r
## Low level of infection:
dl <- df_trends[df_trends $usutuF == "low",]
trl <- 100 * (dl[,4] - dl[,3]) / dl[,3]
tr_low <- c(mean(trl), quantile(trl, c(0.025,0.975)))

## Medium level of infection:
dm <- df_trends[df_trends $usutuF == "medium",]
trm <- 100 * (dm[,4] - dm[,3]) / dm[,3]
tr_medium <- c(mean(trm), quantile(trm, c(0.025,0.975)))

## High level of infection:
dh <- df_trends[df_trends $usutuF == "high",]
trh <- 100 * (dh[,4] - dh[,3]) / dh[,3]
tr_high <- c(mean(trh), quantile(trh, c(0.025,0.975)))
```
```r
# the final data.frame which can be turned into a table for the MS
res_trends<-data.frame(UsutuArea=c("low", "medium", "high"),
                       as.data.frame(rbind(tr_low, tr_medium, tr_high)))
colnames(res_trends)<-c("Usutu area", "Mean trend", "2.5%CI", "97.5%CI")
rownames(res_trends)<-NULL
```

### Results:

```
res_trends
```

```
## Usutu area Mean trend      2.5%CI   97.5%CI
## 1      low      0.7176586  -0.3102161  1.748856
## 2      medium  -7.3722991  -11.1112177  -3.669769
## 3      high    -12.7091388  -16.0061863  -9.122092
```

These values correspond to those presented in the results section of the paper.

## B.5 Summary plot

We finally build the summary graph, corresponding to Fig. 2 of the paper:

```r
## format the data.frame df_4trends for an easier handling
df_stat_abund<-df_4trends |>
dplyr::group_by(usutu_BA, usutuF) |>
dplyr::summarise(averag_abund=mean(mean_abund),
                 CI_low=quantile(mean_abund, PROBS=0.025),
                 CI_high=quantile(mean_abund, PROBS=0.975)) |>
dplyr::mutate_at(c('CI_low', 'CI_high'), as.numeric)
```

```r
## 'summarise()' has grouped output by 'usutu_BA'. You can
## override using the '.groups' argument.
```

```r
## correspondance between Usutu areas and colours
## for the following plots (scale_colour_manual)
col_low<-"palegreen4"
col_medium<-"darkorange2"
col_high<-"deeppink4"
val_colours<-c(col_low, col_medium, col_high)
```

```r
## MAIN PLOT
## A plot of each simulated average values of abundance before and
## after, and the mean and 95%CI for each area
p_trends <- ggplot2::ggplot(df_4trends) +
  ggplot2::geom_line(ggplot2::aes(x=usutu_BA, y=mean_abund,
                                 group=id_sim, colour=usutuF), alpha=0.075) +
  ggplot2::geom_point(data=df_stat_abund,
                     ggplot2::aes(x=usutu_BA, y=averag_abund,
                                     colour=usutuF), size=4) +
  ggplot2::geom_errorbar(data=df_stat_abund,
                        ggplot2::aes(x=usutu_BA, ymin=CI_low, ymax=CI_high,
                                          colour=usutuF), width=0.1, linewidth=1) +
  ggplot2::scale_colour_manual(name="Usutu area", values=val_colours) +
  ggplot2::xlab("Period") + ggplot2::ylab("Abundance index") +
  ggplot2::facet-wrap(~usutuF) +
  ggplot2::theme Environment(text = ggplot2::element Marginal)
```
## Histogram of Usutu densities values for the centroid of each route
## and categorisation in 3 groups

usutu_histo <- ggplot2::ggplot(valagg) +
  ggplot2::geom_histogram(ggplot2::aes(x=usutu�, fill=usutuF)) +
  ggplot2::geom_vline(xintercept=c usutu_p, usutu_75),
  colour="black", linewidth=1,
  linetype=c("solid", "dashed")) +
  ggplot2::ylab("Count") + ggplot2::xlab("Density of Usutu cases") +
  ggplot2::scale_fill_manual(name="Usutu area", values=val_colours)

## Spatial plot of ACT and density of Usutu cases
library(tidyterra)

##

## Attaching package: 'tidyterra'

##

## The following object is masked from 'package:stats':

##

## filter

usutu_sp_act<-ggplot2::ggplot() + tidyterra::geom_spatraster(data=r1) +
  ggplot2::geom_point(data=val p, ggplot2::aes(x=x, y=y, colour=usutuF)) +
  ggplot2::scale_colour_manual(name="Usutu area", values=val_colours) +
  tidyterra::scale_fill_whitebox_c(name="Density of Usutu cases", palette="arid") +
  ggplot2::theme(aspect.ratio=1, axis.text.x=elementos, axis.text.y=elementos, axis.ticks=elementos),
  axis.title.x=elementos, axis.title.y=elementos)

# the figure for the blackbird part

#

##ng("Fig2.png", width=3000, height=3000, pointsize=3, res=300)

library(ggpubr)

##

## Attaching package: 'ggpubr'

##

## The following object is masked from 'package:terra':

##

## rotate

ggpubr::ggarrange(ggpubr::ggarrange (usutu_histo,
usutu_sp_act, ncol=2,
labels=c("A", "B"), widths = c(0.6, 1)),
p_trends, nrow=2, labels=c("", "C"))

## 'stat_bin' using 'bins = 30'. Pick better value with

## 'binwidth'.

##dev.off()
<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

Usutu area

low

medium

high

## References

Baddeley, A., Rubak, E., and Turner, R. (2016). *Spatial point patterns: methodology and applications with R*. CRC press.

Bouchez-Zacria, M., Calenge, C., Lecollinet, S., Villers, A., Eraud, C., and Decors, A. (in prep.). Analysis of the Usutu episode of summer 2018 in birds in France. *PCI infections*.

Diggle, P. (1990). *Time series. A biostatistical introduction*. Oxford University Press, Oxford.

Gabriel, E., Rowlingson, B. S., et Diggle, P. J. (2013). stpp: an R package for plotting, simulating and analyzing Spatio-Temporal Point Patterns. *Journal of Statistical Software*, 53: 1–29.

Manly, B. F. J. (1991). *Randomization and Monte Carlo Methods in Biology*. Springer US.

Wand, M. et Jones, M. (1995). *Kernel smoothing*. Chapman & Hall/CRC.