# Quantitative risk assessment for the introduction of bluetongue virus into mainland Europe by long-distance wind dispersal of Culicoides spp.: A case study from Sardinia

Amandine Bibard, Davide Martinetti, Aymeric Giraud, Albert Picado, Karine Chalvet-Monfray, Thibaud Porphyre

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

## To cite this version:

Amandine Bibard, Davide Martinetti, Aymeric Giraud, Albert Picado, Karine Chalvet-Monfray, et al.. Quantitative risk assessment for the introduction of bluetongue virus into mainland Europe by long-distance wind dispersal of Culicoides spp.: A case study from Sardinia. Risk Analysis, 2025, 45, pp.108-127. 10.1111/risa.14345. hal-04632962v2

HAL Id: hal-04632962

https://hal.science/hal-04632962v2

Submitted on 20 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
<figure><img src="image_3.png" /></figure>

ORIGINAL ARTICLE

# Quantitative risk assessment for the introduction of bluetongue virus into mainland Europe by long-distance wind dispersal of *Culicoides* spp.: A case study from Sardinia

Amandine Bibard<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> ID

Davide Martinetti<fnref n="4" /> ID

Aymeric Giraud<fnref n="4" />

Albert Picado<fnref n="1" /> ID

Karine Chalvet-Monfray<fnref n="3" /> ID

Thibaud Porphyre<fnref n="2" /> ID

Correspondence

Amandine Bibard, Boehringer Ingelheim Animal Health France, Rue de l'Aviation, 69800 Saint-Priest, France.

Email:
amandine.bibard@boehringer-ingelheim.com

Funding information

French National Research Agency and Boehringer Ingelheim Animal Health France, Grant/Award Number: ANR-16-IDEX-0005; Boehringer Ingelheim France; Agence Nationale de la Recherche, Grant/Award Number: Plan de Relance France (ANR-21-PRR-0075-01)

Abstract

Europe faces regular introductions and reintroductions of bluetongue virus (BTV) serotypes, most recently exemplified by the incursion of serotype 3 in the Netherlands. Although the long-distance wind dispersal of the disease vector, *Culicoides* spp., is recognized as a virus introduction pathway, it remains understudied in risk assessments. A Quantitative Risk Assessment framework was developed to estimate the risk of BTV-3 incursion into mainland Europe from Sardinia, where the virus has been present since 2018. We used an atmospheric transport model (HYbrid Single-Particle Lagrangian Integrated Trajectory) to infer the probability of airborne dispersion of the insect vector. Epidemiological disease parameters quantified the virus prevalence in vector population in Sardinia and its potential first transmission after introduction in a new area. When assuming a 24h maximal flight duration, the risk of BTV introduction from Sardinia is limited to the Mediterranean Basin, mainly affecting the southwestern area of the Italian Peninsula, Sicily, Malta, and Corsica. The risk extends to the northern and central parts of Italy, Balearic archipelago, and mainland France and Spain, mostly when maximal flight duration is longer than 24h. Additional knowledge on vector flight conditions and Obsoletus complex-specific parameters could improve the robustness of the model. Providing both spatial and temporal insights into BTV introduction risks, our framework is a key tool to guide global surveillance and preparedness against epizootics.

KEYWORDS

bluetongue virus, HYbrid Single-Particle Lagrangian Integrated Trajectory, midges, orbivirus, risk

1 | INTRODUCTION

Bluetongue disease is recognized as a global veterinary concern due to its economical and animal health impact. It affects both domestic and wild ruminants, and it is mostly transmitted by small flying Diptera (midges) of the genus *Culicoides*. In Europe, since the bluetongue virus (BTV) serotype 10 (BTV-10) first emerged in Portugal and Spain in 1956–1960, 12 additional viral serotypes among the 36 identified nowadays (Ries et al., 2021) have been reported in the continent (BTV-1, -2, -3, -4, -6, -8, -9, -11, -14, -16, -25, -27) (Saminathan et al., 2020). BTV-3 first emerged in

Sicily (2017), then in Sardinia (2018), and very recently in September 2023, in the Netherlands (Holwerda et al., 2023; Wageningen Bioveterinary Research [WBVR], 2023). In Italy, BTV-3 was reported to affect mainly small ruminants with mild clinical signs and remained mostly subclinical in cattle (Cappai et al., 2019). In Sicily, no new clinical cases have been reported in sheep since 2017, but in 2018, 2019, and 2022, one single cattle was diagnosed as BTV-3 positive each year (Bollettino Epidemiologico Nazionale Veterinario (izs.it)). In Sardinia, the virus has been reported circulating every year since 2018 but the infected area has not so far exceeded the endemic zone in the southern part of the island

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2024 The Author(s). *Risk Analysis* published by Wiley Periodicals LLC on behalf of Society for Risk Analysis.

<footnote n="1">¹Global Innovation, Boehringer Ingelheim Animal Health France, Saint-Priest, France</footnote>

<footnote n="2">²Laboratoire de Biométrie et Biologie Évolutives, UMR 5558, Université Claude Bernard Lyon 1, CNRS, VetAgro Sup, Villeurbanne, France</footnote>

<footnote n="3">³Epidémiologie Des Maladies Animales et Zoonotiques, UMR EPIA, Université Clermont Auvergne, INRAE, VetAgro Sup, Saint-Genès-Champanelle, France</footnote>

<footnote n="4">⁴Biostatistique et Processus Spatiaux, UMR 0546, INRAE, Avignon, France</footnote>
(WAHIS (woah.org)). In contrast, the BTV-3 strain circulating in the Netherlands has been associated with severe clinical signs, including deaths, in both sheep and cattle and has spread very rapidly beyond the initial index sites (Holwerda et al., 2023). A pathogenic viral strain of BTV-3 has also previously been reported in Israel (Golender et al., 2020).

Several routes of BTV introduction have been imputed solely or in combination: import of live animals through legal or illegal trade, import of germplasm (semen and/or embryos), movement of wild animals, airborne dispersion of vectors, import of adult vectors by aircraft or use of poorly attenuated modified live vaccines (Gale et al., 2015; Hartley et al., 2013; Hoar et al., 2004; Mintiens et al., 2008; Napp et al., 2011, 2013; Simons et al., 2019). Although the route of introduction into the Netherlands is still under investigation, the suspected route of BTV-3 introduction into Sicily and Sardinia was the long-distance wind dispersal (LDWD) of infected vectors from Tunisia (Aguilar-Vega et al., 2019; Cappai et al., 2019). The BTV-3 circulation in Sardinia, on the doorstep of mainland Europe, poses a real threat to surrounding free countries in Europe. In this study, we investigated which European countries could be at risk of BTV-3 introductions from Sardinia as a direct consequence of LDWD of infected vectors.

Quantitative Risk Assessment (QRA) frameworks aim at describing the sources of risk and quantifying risk levels and their impact. In the last 15 years, their use has significantly increased and has been improved thanks to methodological guidelines (Vos et al., 2011; World Organisation for Animal Health [WOAH], 2010) and computational capacities. Specifically for BTV, introduction risk assessments were developed mostly at country or regional level to enhance national surveillance in high-risk regions or to define high-risk periods of potential introduction (Gale et al., 2015; Gubbins et al., 2010; Hartley et al., 2013; Hoar et al., 2004; Maurella et al., 2019; Napp et al., 2013; Nelson et al., 2022; Roberts et al., 1993; Sagüés et al., 2019; Simons et al., 2019; Sutmoller & Wrathall, 1997). Recently, a generic risk framework at European scale has been developed for BTV considering different pathways of introduction including the long-distance wind dispersion of vectors (Simons et al., 2019). Nonetheless, the source of infection in that study was arbitrary defined as a buffer zone of 300 km around each country without considering variations in topography, wind direction, speed, and vector survival conditions. On the other hand, evidence of long-distance wind dispersion of *Culicoides* spp. has been shown for up to 700 and 500 km over sea (Ducheyne et al., 2011; Eagles et al., 2012, 2014) and over land (García-Lastra et al., 2012) respectively, under suitable conditions. The LDWD of midges depends on meteorological conditions, and these need to be considered in the analyses.

Nowadays, meteorological data are collected worldwide at high resolution, and atmospheric dispersion models exist to simulate trajectories of particles through the atmosphere, mimicking the wind-mediated passive flight of insects. Such meteorological models were previously used to retrospectively identify the potential source of new virus introductions, such as the BTV-2 and BTV-7 introductions into Australia from Oceania islands (Eagles et al., 2014), the BTV-3 introduction into Sicily from Tunisia (Aguilar-Vega et al., 2019), and *Culicoides imicola* migration to mainland France from Corsica (Jacquet et al., 2016). We hypothesize that such an atmospheric dispersion model could suitably infer the probability of long-distance *Culicoides* dispersion from a known source to any destination in Europe and could therefore inform a QRA framework at European scale.

In this study, we have estimated the risk of BTV incursion into mainland Europe via LDWD from Sardinia and identified key parameters to this risk. To do so, we developed a QRA framework considering the complex vector-host-environment interactions of BTV, with multiple host and vector species, and applied this framework to the risk of BTV-3 introduction from the southwestern parts of Sardinia. From our results, high-risk areas and high-risk periods of BTV-3 introduction into European free areas were identified, providing useful guidance for surveillance.

## 2 | MATERIALS AND METHODS

### 2.1 | Framework overview

A QRA framework (Figure 1) was developed based on three evaluation steps:

1. $Inf_{v,i}$: The level of infection in each vector species $v$ at the source $i$, that is, the number of vectors of species $v$ that are infected (and therefore infectious) within the infection source area $i$;

2. $P_{v,ij}$: The probability of LDWD of living vectors of species $v$ from source $i$ to any destination $j$ in Europe;

3. $R_{0j}$: The potential of onward virus transmission at destination $j$, that is, the geometric average rate of secondary infections generated by one primary vector infection at destination $j$.

Our model framework considers two susceptible host species $h$ ("C" for cattle and "S" for small ruminants, the latter comprising sheep and goats), two competent vector species $v$ ("1" for *C. imicola* and "2" for *Obsoletus* complex, the latter comprising *Culicoides obsoletus* and *Culicoides scoticus* species) and the environmental conditions both at the source $i$ and destination $j$.

We estimated the risk of a new BTV introduction into any destination $j$ ($IR_{ij}$) in Europe from source $i$, by computing outputs from the three evaluation steps, as follows:

$$IR_{ij} = \sum_{v \in \{1, 2\}} (Inf_{v,i} \times P_{v,ij}) \times R_{0j}$$

The value $IR_{ij}$ should be interpreted as a daily rate of new infections generated by both vector species $v$ at destination $j$ for each week, and considering the total number of infectious vectors that could be dispersed by the wind from $i$ to
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** Risk Assessment (RA) framework quantifying the risk of bluetongue virus introduction from a source area *i* to a destination area by long-distance wind dispersal (LDWD) of an infected *Culicoides* vector of species *v*. Introduction risk calculation at destination *j* (*IR*<sub>*ij*</sub>) is a combination of three steps.</figcaption></figure>

*j* for each week as the sum of each vector species *v*. Here, *IR*<sub>*ij*</sub> > 1 indicates that a potential introduction of the virus at destination *j* may originate from source *i* and may be established in the local population, whereas *IR*<sub>*ij*</sub> < 1 indicates that no introduction of the virus into destination *j* may occur. The number of weeks with a risk greater than 1 (*W*<sub>*IR*≥1</sub>) and the mean of *IR* over the entire study period (IR̄) were also calculated for each destination *j* and mapped to evidence high-risk regions.

The study period was limited to the 36 weeks from week 11 (mid-March) to week 46 (mid-November), assuming low outdoor *Culicoides* activity during the winter period (Villard et al., 2019). The study grid coordinates were set between longitudes 13°W and 37°E, and latitudes 28°S and 71°N. All datasets were aggregated at the final resolution of 0.5° (grid cell dimensions of approximately 50 × 50 km²) to compute the model. A minimal threshold of 1 vector and 1 host per km² (2500 individuals per population—vector or host—and per each grid cell) was set to compute the model, assuming that BTV transmission is unlikely in areas with too low numbers of host or vector. We decided to use a lower threshold compared to other studies (1 vs. 10 per km²) (Hartemink et al., 2009) because of the high dimension of our grid and in order not to be more conservative.

Table 1 summarizes the model input parameters, their description, and data sources. Details on equations and database computations are provided in Supplement Information S1. Data aggregation and model computation were performed using R Statistical Software v4.2.2 (R Core Team, 2022) and required specific packages, such as “sf” (Pebesma, 2022), “terra” (Hijmans, 2023), “tidyterra” (Hernangómez, 2023), “rnaturalearth” (Massicotte & South, 2023), and “exactextractr” (Baston, 2022).

## 2.1.1 | BTV infection at source (step 1)

The source area *i* was defined as the southwestern part of Sardinia where BTV-3 outbreaks have been reported (Cappai et al., 2019) (Figure 2). The disease was considered at endemic equilibrium, that is, with a bal-

anced number of midges transitioning among the three infection stages—susceptible, latent, and infectious—over time.

The number of susceptible midges was approximated by the daily vector abundance in source *i*. We calculated a daily vector abundance in a given week of species *v* in source *i*, *w* · *N*<sub>*v*,i*</sub>, where *N*<sub>*v*,i*</sub> is the maximal abundance reported within a given year in source *i* and extracted from the VectorNet (Balenghien et al., 2020) data (maximal number per trap per day); and *w* is a probability function of midge presence (ranging between 0 and 1) depending on environmental conditions of the given week (Conte et al., 2009) to account for seasonal variation of the abundance (see Section S1.2).

The number of infected vectors (i.e., the sum of latent and infectious vectors) of species *v* in (*Inf*<sub>*v*,i*</sub>) was then expressed according to the number of susceptible, as follows (see Section S1.1):

$$ Inf_{v,i} = w \cdot N_{v,i} \frac{1}{\mu_v} \sum_{h \in \{C,S\}} \beta_v a_v (\phi_{vh} \text{prev}_h), \quad (1) $$

where *w* · *N*<sub>*v*,i*</sub> is the expected daily number of susceptible vectors in the area *i* in a given week; *μ*<sub>*v*</sub> the mortality rate of vector of species *v*; *β*<sub>*v*</sub> the probability of transmission from a host species *h* to a vector of species *v* given an effective contact; *a*<sub>*v*</sub> the biting rate of vector species *v*; *Φ*<sub>*vh*</sub> the proportion of vectors species *v* attracted by hosts species *h*; and *prev*<sub>*h*</sub> the proportion of infectious host species *h* in area *i* (disease prevalence in host).

The number of infectious hosts in the source *i* was set to 50 cases for cattle, and 200 cases for small ruminant, approximated from the number of BTV-3 events reported by the official surveillance body of Sardinia in December 2022 (Istituto_Zooprofilattico_Sperimentale_della_Sardegna, 2022). The disease prevalence in *i* for each host species *h* (*prev*<sub>*h*</sub>) was calculated by dividing the number of cases by the respective livestock population retrieved in source *i* from the Livestock Grid database (Gilbert et al., 2018) and was considered constant over the study period (Section S1.3 and Figure S3.2). Details on data sources and computation are provided in Section S1.3.
<table><caption>**TABLE 1** Description, formulas, and sources of input parameters used in the risk assessment model of bluetongue virus (BTV) infected *Culicoides* spp. by long-distance wind dispersal (LDWD).</caption><thead><tr><th>Var</th><th>Description</th><th>Comments and formula</th><th>Source/Reference</th></tr></thead><tbody><tr><td>h</td><td>Host population type, C (cattle) or S (small ruminants)</td><td></td><td>—</td></tr><tr><td>v</td><td>Vector population type, 1 (Culicoides imicola) or 2 (Obsoletus complex)</td><td></td><td>—</td></tr><tr><td>N<sub>h</sub></td><td>Host population species h</td><td>Minimal threshold of 1 animal per km<sup>2</sup> was set</td><td>Livestock GLW4 (Gilbert et al., 2018)</td></tr><tr><td>Prev<sub>h</sub></td><td>Disease prevalence in source area in hosts species h</td><td>Prev<sub>h</sub> = number of cases/N<sub>h</sub></td><td>(IZST, 2022)</td></tr><tr><td>N<sub>v</sub></td><td>Maximal vector abundance recorded in a year per trap of species v</td><td>Minimal threshold of 1 vector per km<sup>2</sup> was set</td><td>VectorNet Data Series3 (Balenghien et al., 2020)</td></tr><tr><td>w</td><td>Probability function of vector presence, such as wN<sub>v</sub>, is the expected vector abundance of species v for each week</td><td>Depend on temperature, humidity, and elevation</td><td>(Conte et al., 2020), Supplement Information S1</td></tr><tr><td>Inf<sub>v,i</sub></td><td>Number of infected vector of species v at source i</td><td>Equation detailed in main text</td><td>(Turner, Bowers & Baylis, 2013), Supplement Information S1</td></tr><tr><td>U<sub>i</sub></td><td>Probability to be uplifted in the air mass</td><td>Constant of 0.25</td><td>—</td></tr><tr><td>H<sub>ij</sub></td><td>Probability of wind dispersal from source i to destination j inferred from HYSPLIT simulations</td><td>Mean frequency of deposition in each destination cell over the total forward simulations run by HYSPLIT model from source i, averaged for each week over the 3 years 2020–2022</td><td>(Stein et al., 2015)</td></tr><tr><td>a<sub>v</sub></td><td>Biting rate of vector species v</td><td>a<sub>1</sub>(T) = 0.00014 × T × (T − 3.6966) × (41.8699 − T)<sup>2.7056</sup> for T > 3.7°C<br>a<sub>2</sub>(T) = 0.000171 × T × (T − 3.6966) × (41.8699 − T)<sup>1.27056</sup> for T > 3.7°C</td><td>(Aguilar-Vega, Bosch, Fernandez-Carrion, Lucientes, & Sanchez-Vizcaino, 2020; Mullens, Gerry, Lysyk & Schmidtmann, 2004)</td></tr><tr><td>ω<sub>v</sub></td><td>Virogenesis rate (i.e., rate at which latent vector species v become infectious = 1/extrinsic incubation period in days)</td><td>ω<sub>1</sub>(T) = 0.016 × (T − 12.7) for T > 12.7°C<br>ω<sub>2</sub>(T) = 0.0003 × T × (T − 10.4057) for T > 10.4°C</td><td>ω<sub>1</sub> (Sanders et al., 2011)<br>ω<sub>2</sub> (Mullens et al. 2004; Turner et al., 2013)</td></tr><tr><td>μ<sub>v</sub></td><td>Natural mortality rate of vector species v</td><td>μ<sub>1</sub>(T) = 1 − (1e − 5 × T<sup>3</sup> − 0.001 × T<sup>2</sup> + 0.0187 × T + 0.8924) for T > 10°C<br>μ<sub>2</sub>(T) = 0.015 × exp (0.063 × T) for T > 0°C</td><td>μ<sub>1</sub> (White, Sanders, Shortall, & Purse, 2017)<br>μ<sub>2</sub> (Wittmann, Mellor, & Baylis, 2002)</td></tr><tr><td>β<sub>v</sub></td><td>Probability of transmission from host to vector of species v</td><td>β<sub>1</sub>(T) = 0.0003699 exp (0.1725 × T)<br>β<sub>2</sub>(T) = 0.005465 exp (0.159 × T)</td><td>(Turner et al., 2013)</td></tr><tr><td>m<sub>hv</sub></td><td>Ratio of vectors species v to hosts species h</td><td>m<sub>vh</sub> = wN<sub>v</sub>/N<sub>h</sub></td><td>—</td></tr><tr><td>Φ<sub>hv</sub></td><td>Proportion of vectors species v attracted to hosts species h</td><td>Φ<sub>Cv</sub> = m<sub>Sv</sub>/(m<sub>Sv</sub> + σ<sub>v</sub> × m<sub>Cv</sub>)<br>Φ<sub>Sv</sub> = 1 − Φ<sub>Cv</sub><br>Φ<sub>h1</sub> = 1 − Φ<sub>h2</sub></td><td>—</td></tr><tr><td>σ<sub>v</sub></td><td>Host preference for vector species v</td><td>σ<sub>1</sub> = σ<sub>2</sub> = 0.15 (σ<sub>v</sub> < 1 indicates a preference for cattle)</td><td>(Aguilar-Vega et al., 2020)</td></tr><tr><td>r<sub>h</sub></td><td>Recovery rate (1/ duration of viremia in days) in host species h</td><td>r<sub>C</sub> = 1/29<br>r<sub>S</sub> = 1/20</td><td>(Dorea et al., 2017)</td></tr><tr><td></td><td></td><td></td><td>(Continues)</td></tr></tbody></table>
<table><caption>TABLE 1 (Continued)</caption><thead><tr><th>Var</th><th>Description</th><th>Comments and formula</th><th>Source/Reference</th></tr></thead><tbody><tr><td>d<sub>h</sub></td><td>Disease induced mortality in host h</td><td>d<sub>C</sub> = 0<br>d<sub>S</sub> = 0.0078</td><td>d<sub>C</sub> (Guis et al., 2012)d<sub>S</sub> (Aguilar-Vega et al., 2020)</td></tr><tr><td>b</td><td>Probability of transmission from vector to host</td><td>b = 0.9</td><td>(Aguilar-Vega et al., 2020; Hartemink et al., 2009)</td></tr></tbody></table>

Note: Rates are given per day.

Abbreviation: HYSPLIT, HYbrid Single-Particle Lagrangian Integrated Trajectory.

<figure><img src="image_5.png" /><figcaption>FIGURE 2 Source area *i* was defined in the southwest of Sardinia. Atmospheric simulations in HYbrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) were initiated from 290 locations located on land surface of this area.</figcaption></figure>

## 2.1.2 | Long-distance wind dispersal (step 2)

Step 2 consisted of estimating the probability of LDWD from source to destination ($P_{v,ij}$), which is the combination of the probabilities to be uplifted in the windstream ($U_i$), to be subsequently transported by the wind to a destination area $j$ ($H_{ij}$), and to survive along the aerial transport until $j$ ($S_{v,j}$), as follows:

$$P_{v,ij} = U_i \times H_{ij} \times S_{v,j}, \quad (2)$$

where $U_i$ estimated the proportion of vectors that will be uplifted in the wind stream in average per week and subsequently long-distance transported. In absence of reference in the literature, we estimated this proportion as fixed (0.25)

and corresponding to the active proportion of midges flying at sunset and sunrise (see the following paragraph on conditions for initiation of atmospheric simulations).

An atmospheric dispersion model, HYSPLIT (Stein et al., 2015) (HYbrid Single-Particle Lagrangian Integrated Trajectory), was used to infer $H_{ij}$. Forward simulation trajectories were initiated from 290 starting points comprised within land surface of the source area *i*, which covers approximately 3300 km² (resolution 0.5°) (Figure 2).

Trajectory simulations were run daily over the 36 weeks of the study period. In each day, wind trajectories simulations were initiated from nine timepoints corresponding to periods of flight activity of vectors (Elbers et al., 2016): around dusk (at sunset and 1 h before and after) and at dawn (at sunrise and in the following 5 h). Trajectories were set to
start at a minimal altitude of 50 m above ground level (a.g.l.). Along each simulated trajectory, meteorological and atmospheric variables were recorded hourly, including air mass temperature and thickness of the planetary boundary layer (PBL). Each hourly checkpoint was deemed a viable destination location for the vector if the temperature was greater than 10°C (Braverman & Chechik, 1996; Sellers & Maarouf, 1991) and within the PBL; otherwise, further checkpoints of the trajectory were discarded. Trajectories were stopped based on two separate flight durations (dispersion): 24 and 48 h.

Finally, $H_{ij}$ was estimated by dividing the number of times per week the destination j was reached by the total number of destinations from all simulations initiated from i during this week. Simulations were run from week 11 to week 46 (calendar weeks), independently for 2020, 2021, and 2022, and the final $H_{ij}$ value represents the average over the years for each calendar week.

The duration of passive flight in our scenarios (24 or 48 h) was considered not negligible as compared to the lifespan of the vector, and the probability to be still alive at destination j, $S_{v,j}$, was expressed as follows:

$$S_{v,j} = e^{-\mu_{v,j}} \times \text{time} \quad (3)$$

With time as the two scenarios of dispersion in days (24 or 48 h) and $\mu_{v,j}$ the mortality rate of vector species v calculated with environmental conditions in destination j.

For the sake of clarity in the presentation of results, $P_{v,ij}$ was only presented considering the conditions of Obsoletus complex and in this case was named $P_{ij}$. The expected probability of LDWD ($P_{ij}$) was computed for each given week of the study period and considering environmental conditions constant over the week. Furthermore, we computed the overall probability of LDWD ($\bar{P}_{ij}$) as the mean of $P_{ij}$ over the whole study period. Both the expected probability $P_{ij}$ and the overall probability $\bar{P}_{ij}$ were then binned into four intervals ($[0, 10^{-4}]$; $[10^{-4}, 10^{-3}]$; $[10^{-3}, 10^{-2}]$; $[10^{-2}, 10^{-1}]$) to map four dispersion rings around the source area.

### 2.1.3 | BTV transmission at destination (step 3)

Disease transmissibility in each destination cell j was estimated through the $R_{0j}$ (basic reproduction number), usually defined as the expected number of secondary cases produced, in a completely susceptible population, by a typical infectious individual (Driessche & Watmough, 2002). We used the $R_0$ formula defined for BTV by Turner et al. (2013) considering two-host and two-vector populations v, with $\in\{1, 2\}$ for C. imicola and Obsoletus complex and $h\in\{C, S\}$ for cattle and small ruminant populations, respectively. It combines four $R_{xy}$ components, with x defining the vector species at destination j that get infected due to one primary infectious

vector of species y incoming from source i, such as

$$R_{0j} = \sqrt{\frac{1}{2} \left[ (R_{11} + R_{22}) + \sqrt{(R_{11} + R_{22})^2 - 4(R_{11}R_{22} - R_{12}R_{21})} \right]}, \quad (4)$$

with

$$R_{11} = \left( \frac{\beta_1 ba_1^2}{\mu_1} \right) \left( \frac{\omega_1}{\omega_1 + \mu_1} \right) \left( \frac{m_{c1}\phi_1^2}{r_c + d_c} + \frac{m_{s1}(1 - \phi_1)^2}{r_s + d_s} \right),$$

$$R_{22} = \left( \frac{\beta_2 ba_2^2}{\mu_2} \right) \left( \frac{\omega_2}{\omega_2 + \mu_2} \right) \left( \frac{m_{c2}\phi_2^2}{r_c + d_c} + \frac{m_{s2}(1 - \phi_2)^2}{r_s + d_s} \right),$$

$$R_{12} = \left( \frac{\beta_1 ba_2 a_1}{\mu_2} \right) \left( \frac{\omega_2}{\omega_2 + \mu_2} \right) \left( \frac{m_{c1}\phi_2\phi_1}{r_c + d_c} + \frac{m_{s1}(1 - \phi_1)(1 - \phi_2)}{r_s + d_s} \right),$$

$$R_{21} = \left( \frac{\beta_2 ba_2 a_1}{\mu_1} \right) \left( \frac{\omega_1}{\omega_1 + \mu_1} \right) \left( \frac{m_{c2}\phi_2\phi_1}{r_c + d_c} + \frac{m_{s2}(1 - \phi_1)(1 - \phi_2)}{r_s + d_s} \right),$$

where $R_{11}$ is the average number of infectious C. imicola produced at destination by one infectious C. imicola incoming from source; $R_{22}$ the average number of infectious Obsoletus complex midges produced at destination by one infectious Obsoletus complex midge incoming from source; $R_{12}$ is the average number of infectious C. imicola produced at destination by one infectious Obsoletus complex incoming from source; and $R_{21}$ is the average number of infectious Obsoletus complex produced at destination by one infectious C. imicola incoming from source. For simplicity, the index j of each $R_{xy}$ of $R_{0j}$ was not expressed but is implicit.

It is worth noting that $R_{0j}$ is actually the square root of the “Type Reproduction Number” $R_T$ defined for vector-borne disease (Yakob & Clements, 2013). The $R_{0j}$ is equivalent to the average number (geometric mean) of secondary infections at destination produced by one single infection. As such, an estimate of $R_{0j} < 1$ implies that disease transmission would not be sufficient for BTV to sustain in these conditions, whereas BTV would spread where $R_{0j} > 1$. The terms $\sqrt{R_{11}}$ and $\sqrt{R_{22}}$ represent the case where BTV would be transmitted by only one-vector species, C. imicola or the Obsoletus complex, respectively. All infected vectors originating from i and surviving movement to j were considered infectious in j, and all local hosts in j were considered susceptible to the disease. Spatial distributions of host populations in Europe as considered in our study were extracted from the Gridded
Livestock of the World (GLW3) (Gilbert et al., 2018) and are shown in Figure S3.2.

Input parameters involved in Equation (4) are described in Table 1. Parameters specific to the vector species v, like the virogenesis rate ($\omega_v$), the mortality rate ($\mu_v$), or the probability of virus transmission from a host to vector ($\beta_v$), are dependent on the temperature, and the weekly mean in each destination cell j was used. The probability of virus transmission from a vector to a host (b) was considered constant, irrespective of the host and vector species involved (Aguilar-Vega et al., 2020; Hartemink et al., 2009). The host recovery rates ($r_h$), defined as the inverse of the viremic period in host, were retrieved from a meta-analysis performed by EFSA (Dórea et al., 2017). The median end day of viremia detection (by PCR or virus isolation) was taken as reference for each host species (29 days for cattle and 20 days for sheep and goats; assuming equal parameters for small ruminants). Disease mortality rates in hosts ($\mu_h$) are highly dependent of the species susceptibility and the viral serotype but were assumed to be low (0.0078) in small ruminants and even null in cattle (asymptomatic) (Maclachlan, 2011; Saminathan et al., 2020).

## 2.2 | Sensitivity analysis

A global sensitivity analysis was performed to identify the most influential parameters of the model. We used the Sobol' Indices method (Saltelli et al., 2010) which is a variance-based statistical method and does not require a priori assumptions (linearity or additivity) of model behavior. First, maximal range domains for each of the 18 parameters tested were defined according to literature, and uniform distributions were applied for all factors (Supplement Information S2). Maximal flight duration and host prevalence were fixed at 24 h and 0.01, respectively. Second, the Sobol indices were estimated using the revised method of Jansen and Sobol (Saltelli et al., 2010), following 4-step procedure: (1) two Latin hypercube samplings of each factor with 500,000 iterations, (2) computation of sampled values leading to $10^7$ sets of parameters combinations, (3) model run for each set of parameters combinations, and (4) calculation of the Sobol sensitivity indices with 100 bootstrap samples. First-order indices represent the impact of each individual parameter change on model output, whereas the total-order indices account for potential interactions among parameters. Sobol indices of each variable were depicted in a tornado chart.

For the parameters identified as most influential, we ran additional model iterations to characterize their individual impact on risk output, considering independent sampling one parameter at a time and a constant environment for the others. Specific R packages “mc2d” (Pouillot, 2022), “lhs” (Carnell, 2022), and “sensitivity” (Iooss et al., 2022) were used for the analysis.

## 2.3 | Uncertainty assessment

Once the most influential input parameters had been identified, we quantitatively estimated how the uncertainty around their values could affect the accuracy of our model outputs, in both time and space. The general approach was to introduce variation in each influencing parameters according to their level of uncertainty, iterate model runs (500 iterations), and quantify statistical quantities among iterations of risk outputs. The methodology used to introduce variation (random sampling on probabilistic distributions) depended on whether influencing parameters were either published as a static point estimate, derived from an environment-dependent function, or derived from a database (see details in Supplement Information S.2). In the absence of uncertainty indications in literature, noise of 10% (Yakob & Clements, 2013) was arbitrary introduced. For each destination j and for each week of the study period, we calculated $P(IR > 1)$, the proportion of iteration runs where $IR_{ij} > 1$.

# 3 | RESULTS

## 3.1 | BTV infection at source (step 1)

Over the 3300 km² source area, the maximal abundance of Obsoletus complex midges was 1.8 in $log_{10}$ higher than for *C. imicola* midges ($1.6 \times 10^6$ vs. $2.4 \times 10^4$ individuals) (Figure 3). Disease prevalence set at 0.06% in small ruminants and 0.66% in cattle led to a median number of 2.81 in $log_{10}$ (equivalent to 654 individuals) infected Obsoletus complex midges per week and -0.04 in $log_{10}$ (equivalent to 0.9 individuals) infected *C. imicola* midges per week, but with a significant variation range over the study period (Figure 3). Infection in the vector populations started to increase since mid of April (week 18) and peaked during summer (weeks 26-34—end of June to end of August) with a number of 2480 Obsoletus complex and 2 *C. imicola* infected vectors per week. The number of infected vectors started to decrease smoothly from September to mid-November. Maximal infection rates, observed in July, were $1.6 \times 10^{-3}$ and $8.2 \times 10^{-5}$ for Obsoletus complex and *C. imicola*, respectively.

## 3.2 | Probability of long-distance wind dispersal (step 2)

Figure 4 depicts the four dispersion rings of the overall probability of LDWD $\bar{P}_{ij}$ from southwestern Sardinia. The dispersion ring (1) of highest probability $[10^{-2}-10^{-1}]$ was restricted to southern part of Sardinia. The dispersion ring (2) of probability between $10^{-3}$ and $10^{-2}$ covers half of Sardinia but never reached any other island. The dispersion ring (3) of probability between $10^{-3}$ and $10^{-4}$ covered western part of Sicily, the southwestern point of Italy, Malta, the Balearic
<figure><img src="image_6.png" /><figcaption>**FIGURE 3** Weekly estimates of vector abundance (A) and infection (B) in the source area *i* (Sardinia), for the two-vector species, in log10 scale. On the x-axis, month name (in brackets) for each associated week number is indicative of the 2020 calendar.</figcaption></figure>

Islands, and the northern and eastern borders of Tunisia, but not Corsica. The lowest probability dispersion ring (4) reached extreme destinations such as Paris in France, south-eastern border of Spain, or eastern countries like Kosovo or Bosnia and Herzegovina.

When considering longer time of passive flight (Figure S3.3), the lowest probability dispersion ring (4) was much more broadly extended. After 48 h of dispersion, destinations as far north as southern Sweden (~1800 km to Malmö) or the border between Ukraine and Belarus (~2000 km) could be reached, but with very low probability. No western continental country was reached with a probability $\bar{P}_{ij} > 10^{-3}$. Surprisingly, the overall probability of LDWD to Corsica was <10⁻⁴, whatever the scenario, despite being less than 260 km away from southwestern Sardinia.

However, the spatial distribution of the overall expected probability LDWD $\bar{P}_{ij}$ from southwestern Sardinia hides marked temporal variations throughout the study period (Figure 5). Unlike in Figure 4, dispersion rings of $P_{ij}$ pointed to the western direction at the beginning of the study period, mainly from mid-March to mid-May (weeks 11–19) or by October (weeks 38 and 40). The dispersion ring (3) of $P_{ij}$ between 10⁻³ and 10⁻⁴ occasionally reached the Mediterranean coast of continental France (weeks 19, 24, and 37) and made a significant incursion over the southwestern part of the country (provinces of “Pyrénées Orientales,” “Aude,” and “Hérault”) at weeks 42 and 43.

When considering 48 h of wind dispersal (Figure S3.4), the incursions over mainland France by the dispersion ring (3) were more frequent (observed 13 times over 36 weeks) and were more marked at weeks 24–25 (June) and weeks 42–43 (third week of October).

## 3.3 | BTV transmission at destination (step 3)

Destinations with $R_{0j}$ greater than 1 constitute areas where BTV is likely to be transmitted to local populations, considering two vector and host populations. Figure 6 shows the spatial distribution of $\bar{R}_{0j}$, the mean of $R_{0j}$ over the study period. In Italy, the whole peninsula was dotted with high-risk areas of transmission, particularly on the coastal line and rarely on the northern region. In France, Corsica and three areas in the mainland country exhibited a high risk of transmission: southeastern Mediterranean coastal line, southwestern area, and northern-central region surrounding the south of Paris. In Spain, the northeastern Mediterranean coast and the northwestern area were shown to be more at risk than the rest of the country. High values of $\bar{R}_{0j}$ were also estimated on the eastern coastal line of the Adriatic Sea and the southwestern border of the Black Sea. Unsurprisingly, null or very low values were estimated in mountainous areas like in the Alps, the Pyrénées, central Spain, Scandinavian, and Balkan Mountains.

Maps of $R_{0j}$ computed for each week of the study period (Figure S3.5) highlight seasonal variations of $R_{0j}$ in line with variations of vector abundance. The disease transmission was unlikely ($R_{0j} < 1$) in most of European destinations before week 17 (mid of April) and after week 44 (end of October). Although $R_{0j}$ integrate the role of both vector species in the spread of BTV, it appears to be mostly driven by Obsoletus complex midges (Figure S3.6), probably due to its higher abundance and widespread distribution in comparison with that of *C. imicola* (Figure S3.1).
24-hours dispersion scenario

<figure><img src="image_7.png" /><figcaption>**FIGURE 4** Spatial distribution of the overall mean probability of long-distance wind dispersal from southwestern Sardinia to any destination in Europe (step2), considering 24-h maximal flight duration and the survival conditions of Obsoletus complex ($\bar{P}_{ij}$ equivalent to $\bar{P}_{2, ij}$ here). Here, overall mean probability $\bar{P}_{ij}$ is computed over the considered 36-week study period.</figcaption></figure>

## 3.4 Introduction risk outputs ($IR_{ij}$)

Figure 7A shows the spatial distribution of $\overline{IR}_{ij}$, the mean of $IR_{ij}$ over the study period. Unsurprisingly, the highest estimates of risk of BTV introduction, that is, $\overline{IR}_{ij} > 10$, were restricted to the southern part of Sardinia, whereas dispersion to the northern part of the island was characterized by a $\overline{IR}_{ij}$ between 1 and 5. Estimates of $\overline{IR}_{ij} > 1$ were found in Sicily, Malta, the southwestern coast of Corsica, and on the southern tip of Italy. When considering longer dispersion duration (Figure S3.10), destination cells with $\overline{IR}_{ij} > 1$ extended west-

ward to Minorca (Balearic archipelago), northward to the Mediterranean coast of continental France, and eastward to sporadic destinations within the central and southern parts of the Italian peninsula. No risk above 1 was reported for any European countries/regions above latitude 45°, whatever the scenario of dispersion. Interestingly, several areas along the northern coast of Tunisia and Algeria showed a risk of BTV introduction $\overline{IR}_{ij} > 1$.

Temporal dynamics of the risk of BTV introduction over the study period and for the different scenarios of maximal flight duration are shown in Figure S3.8 and S3.9. Considering a 24-h maximal flight duration, some areas sporadically
<figure><img src="image_8.png" /><figcaption>**FIGURE 5** Spatial distribution of the mean probability of long-distance wind dispersal from southwestern Sardinia, considering 24-h maximal flight duration and the survival conditions of Obsoletus complex (*P*̄*ij* equivalent to *P*̄2, *ij* here). Here, mean probability *P*̄*ij* is computed for each week of the study period, considering constant environmental conditions.</figcaption></figure>

experienced a *IR*<sub>ij</sub> >1 for a few weeks (Figure 7B). In mainland France, the Mediterranean coast and the southwestern region were exposed at week 24 (mid-June), week 37 (mid-September), and weeks 42–43 (end of October). In Spain, Minorca and, to a lesser extent, Mallorca showed a risk above 1 during 1–10 weeks within the June to September period. In contrast, continental Spain showed only limited risk of BTV incursion from Sardinia, with its extreme eastern point showing a *IR*<sub>ij</sub> >1 only once or twice over the study period. The southern and central parts of the Italian peninsula were sporadically exposed to risk greater than 1 at least once during

the study period, but the central part was exposed later in the year (in October, from weeks 39 to 44).

As dispersal duration increased (scenarios of 48-h maximal flight duration—Figure S3.9), estimates of *IR*<sub>ij</sub> between 1 and 5 extended to the whole south of Corsica, the Mediterranean border of mainland France, and sporadically on the Italian peninsula.

Considering *C. imicola* alone as a competent vector to transmit the disease, overall risk of further incursion *IR*<sub>ij</sub> remains below 1, whatever the scenario of dispersal duration considered (Figure S3.7).
<figure><img src="image_9.png" /><figcaption>**FIGURE 6** Spatial distribution of the overall transmissibility $\bar{R}_{0j}$ of bluetongue virus (BTV) in Europe considering two-host types (cattle and small ruminants) and two-vector populations (*Culicoides imicola* and Obsoletus complex). Here, the values of $\bar{R}_{0j}$ in each grid cell *j* are computed as the mean of $R_{0j}$ over the study period and considering constant weekly environmental conditions.</figcaption></figure>

**3.5 | Sensitivity analysis**

As expected, uncertainties in the probability to be uplifted in the windstream ($U_i$) and the HYSPLIT probability to be airborne to a new destination ($H_{ij}$) have the strongest impact on the model outputs, both independently (first effect) and in interaction with the other variables (total effect) (Figure 8). Uncertainties in parameters specific to the dynamics of Obsoletus complex midges and their role in disease transmission, such as biting rate ($a_2$), abundance ($wN_2$), transmission rate

from host to vector ($\beta_2$), and mortality rate ($\mu_2$), have also a significant impact on the model outputs, mostly when considering the interactions with other parameters. Their individual impacts on the risk outputs (function shapes and directions) are presented in Supplement Information S2.2.

Uncertainty on the *C. imicola*-specific parameters, mostly the biting and mortality rates ($a_1, \mu_1$), as well as the transmission rate from vector to host (*b*), showed a limited impact on the model. The vector preference for a host species ($\sigma_v$), the host-specific parameters (mortality and recovery rates of both
<figure><img src="image_10.png" /><figcaption>**FIGURE 7** Spatial distribution of the overall risk of introduction from southwestern Sardinia to any destination in Europe, considering a 24-h maximum flight duration of vectors. Here, the risk of introduction was expressed as (A) the overall risk of introduction (*IR*) in Europe and (B) the number of weeks in which the mean risk of introduction *IR* is greater than 1 (*W*<sub>*IR*≥1</sub>). Here, the mean risk of introduction *IR* was computed considering constant weekly environmental conditions and a 24-h maximum flight duration of vectors.</figcaption></figure>

host species), and the virogenesis rates of both species of vectors ($\omega_1$, $\omega_2$) do not have any impact on the model variance.

**3.6 | Uncertainty assessment around risk**

Figure 9 highlights the destination cells that have a probability greater or equal to 75% of being at risk of introduction $P(IR > 1) \ge 0.75$, when accounting for uncertainty in influential parameters (full estimates in Figure S3.11). In this analysis, no destination could be identified at risk ($IR > 1$), with at least 75% confidence, before week 22 (end May/beginning June) and after week 44 (end October), which is 5 weeks less than for the baseline period when uncertainty is not considered (from week 19—beginning May—to week 46—mid-November). Sicily and the southern tip of the Italian peninsula remained the main risk areas, although the central and northern parts of the peninsula were exposed sporadically for only 1 or 2 weeks. Interestingly, Northern coast of Tyrrhenian Sea and Corsica were simultaneously at risk on week 30 (end of July). In addition, southwest of Corsica is also exposed later in the summer.

Week 24 (mid-June) was the only week during which mainland France was exposed to a risk greater than 1 with a high confidence level. For Minorca and the southeastern Mediterranean border of mainland Spain, the riskiest week was week 36 (September).

**4 | DISCUSSION**

In this study, we quantified the risk of BTV-3 introduction as a mean daily rate of new infections expected at destination

per week, proportionally to the number of vectors that were infected in the southwest of Sardinia, successfully brought by the wind and did infect local populations at destination. Variations in time and space of this risk estimate were used to predict areas and periods at risk throughout Europe.

First, our results confirmed a high risk of northward expansion of BTV within Sardinia Island, which was historically evidenced for other BTV serotypes (Calistri et al., 2004) but surprisingly not for BTV-3 since 2018. One hypothesis could be that the Obsoletus complex midges have a reduced competence (lower transmission rate from host to vector) against this viral serotype than the one used in our model. No specific vector competence study was done specifically for Obsoletus complex for this strain, and recent entomological trappings in the North of the source area failed to identify Obsoletus complex midges infected by BTV-3, contrary to *C. imicola* (Quaglia et al., 2023). Under the hypothesis that Obsoletus complex would not be competent to replicate the Sardinian BTV-3 strain, the unique *C. imicola* population seems not sufficient to constitute a risk of BTV spread by the wind (Supplement Information S3.7). However, the recent BTV-3 outbreak in the Netherlands (WBVR, 2023) supports the fact that Obsoletus complex midges might be able to transmit BTV-3 virus strains, at least the one circulating in the Netherlands. Further studies on the vector competence of the Obsoletus complex with respect to both BTV-3 strains would be of particular interest to fine-tune our model.

Out of Sardinia, the riskiest destinations for further dispersion are Sicily, Malta, and the southern tip of Italy, due to a major windstream blowing in that direction. This result was observed at both temporal scales (weekly and annually) and remained robust despite uncertainty of the input parameters. This southeastern spread is in line with history of BTV-2
<figure><img src="image_11.png" /><figcaption>**FIGURE 8** Results of sensitivity analysis ranking most influential parameters on model variability.</figcaption></figure>

which emerged in August 2000 in Sardinia and was success-
sively introduced in Sicily and in the Calabria region of the
Italian peninsula (Calistri et al., 2004). Corsica and Minorca
were also shown to be at risk of introduction but more spo-
radically during specific weeks of the study period. Here
again, it appears consistent with historical cases of BTV-4,
which was first detected in Sardinia (August 2003) and sub-
sequently reported on Corsica and Minorca (October 2003)
(Gómez-Tejedor, 2004). Regarding Balearic archipelago, our
results suggest a lower risk for Mallorca than for Minorca,
whereas historical cases of BTV-2 were first reported on Mal-
lorca and subsequently on Minorca in 2000 (Calistri et al.,
2004; Diego et al., 2014). However, it should be noted that
BTV-2 was concomitantly circulating in North Africa in
2000, suggesting that these countries may have been the
original or additional sources of BTV-2 introduction to the
Balearics (Burgin et al., 2013). Additionally, our results sup-
port a direct introduction from Sardinia to Tuscany region

in Italy as was suggested by Purse et al. (2008) for BTV-
2. The risk of introduction into mainland France and Spain
was low, sporadic, and uncertain considering a 24-h max-
imal flight duration but increases when the maximal flight
duration was considered longer. Although these regions have
never reported cases of BTV-3, they should be considered
surveillance hotspots. The first cases of BTV-3 in the Nether-
lands were diagnosed in the region of Utrecht (WBVR, 2023),
in early September 2023. As this location was not identified
by our model as a high-risk destination for the introduction
of BTV-3 by long-distance dispersal of vectors from Sar-
dinia, whatever the scenario of maximal flight duration, it
can be suggested to investigate other routes of introduction
such as trade of live animal and germplasm or aircraft move-
ments. Finally, despite the lack of a true validation process,
our framework was able to predict variations of risk in time
and space, consistent with previous BTV incursion to guide
surveillance at large scale.
<figure><img src="image_12.png" /><figcaption>**FIGURE 9** Areas in Europe showing a high estimated mean risk of introduction from southwestern Sardinia with high level of certainty. Uncertainty is measured as the probability $P_{IR>1}$ for which the estimate of IR is greater than 1, considering uncertainty around key parameters. We considered that values of $P_{IR>1} \ge 0.75$ show a high level of certainty. Here, estimates of IR were computed for each week of the study period and considering constant weekly environmental conditions and a 24-h maximum flight duration of vectors.</figcaption></figure>

Although our model shows consistent results with historical empirical evidences, we acknowledge several limitations at each step of our model framework.

First, vector abundance data were estimated from VectorNet database (Balenghien et al., 2020) which has the advantages to be open access and gathering light trappings data from a huge entomologist surveillance network covering 28 countries of continental Europe and surrounding regions, from 2015 to 2018 (Balenghien et al., 2020). We assumed that the VectorNet metric (maximum number of midges per trap per day) could be equivalent to an average “aggressive” density (Balenghien et al., 2020) per km² (or a mean density per farm in a pixel of 1 km²), with a constant renewal of the vector population. Knowing that trappings data do not represent the absolute *Culicoides* abundance (Meiswinkel & Elbers, 2016), and the attraction range of a light trap is lower than 1 km² (Kirkeby et al., 2013) and varies between species (Viennet et al., 2013), vector abundance estimates used in the model might be underestimated. Moreover, the maximal abundance retrieved for *C. obsoletus/scoticus* in the source area was about 63 times higher than for *C. imicola*,

whereas Italian field trappings normally collect higher *C. imicola* than *C. obsoletus* in this low-altitude area of Sardinia (Foxi & Delrio, 2010; Foxi et al., 2016). This discrepancy may suggest that the VectorNet's algorithm, despite its overall good performance, tends to overestimate abundance for vectors that are ubiquitous and widely distributed such as Obsoletus complex, and for which, the lack of data reporting an actual absence compromises the predictive performance of the model. Weekly abundances were calculated applying a presence probability function *w* to the maximal vector abundance, and this function was applied for both vector species at European scale, whereas it was initially established for *C. imicola* in Italy (Conte et al., 2003). However, the large activity period from spring to late autumn was consistent with the one usually reported. But we could not reproduce, specifically for Obsoletus complex in southern areas, the reduced peak of activity in summer when the conditions are dry and hot (Versteirt et al., 2017). For all these reasons, our estimates of weekly abundance might be biased, particularly for Obsoletus complex, and should be improved by continuous trapping collections and better understanding around its ecology.
Second, disease prevalence in host in source *i* (0.06% in small ruminants and 0.66% in cattle) was considered constant over the study period in absence of official reports of BTV-3 cases per month. Although these estimates can be classified as “low” and “moderate,” respectively, according to the MintRisk scoring for vector-borne disease (Vos et al., 2021), they might overestimate the infection pressure when the vector abundance is low, particularly in the beginning of the period. Despite this conservative assumption, the final risk outputs remained low.

Third, a higher infection rate in Sardinia was calculated for Obsoletus complex than for *C. imicola* (maximum rate of 0.16% and 0.008%, respectively). Although this is in line with more favorable transmission parameters suggested for this first species (Foxi et al., 2016), comparison with field estimations of infection rate derived from vector capture suggests that our estimates appear lower than those reported during the BTV-1/BTV-4 epidemy in Sardinia (Goffredo et al., 2015) (3.6% for Obsoletus complex and 1.7% for *C. imicola* specifically) and little higher than estimates for BTV-9 from the Italian peninsula for Obsoletus complex (0.005%) (Braks et al., 2017; Savini et al., 2005). Vector infection rates reported in literature are heterogeneous due to diverse calculation methodologies (Braks et al., 2017) (pool vs. individuals) or disease context (epidemic vs. endemic situation) and differ a lot between virus serotypes. Specific BTV-3 laboratory studies assessing infection, dissemination, and transmission rates would be valuable for our model parametrization.

One of the greatest strengths in our methodology is to forecast vector airborne dispersion from atmospheric simulations (step 2) based on real meteorological data and considering their variations at a fine spatial (grid cells of 2500 km²) and temporal (weekly) scales. This enables us to precisely identify wind pathways that would be overlooked at higher resolution or timescales. Indeed, HYSPLIT simulations predicted some sporadic wind incursions to the northwest direction toward continental France in June and October, despite a global wind pattern pointing in the southwestern direction. These episodic events may be linked to the northward Sirocco wind flow coming from Sahara Desert, known to result in hot and dusty climatic conditions in European countries such as Spain (Sousa et al., 2019), Portugal, and France. Global climate models have identified the Mediterranean area as a hot spot for climate change with winter rain decline and increasing wind and temperature anomalies (Tuel & Eltahir, 2020). To integrate these future meteorological changes into our model and to guarantee reliable estimations, regular updates of HYSPLIT simulations will need to be performed. Importantly, the sensitivity analysis confirmed that the probability of wind dispersion is one of the most influencing parameters on the risk output (Figure 8). This result highlights the importance of the methodological filters applied to raw simulations to ensure that final trajectories best mimic the flight of a living and active midge. First, the impact of the maximal flight duration was assessed through two scenarios considering a maximal flight duration lasting

for 24 or 48 h, unlike most of the previous studies that defined only one duration, between 12 and 72 h (Alba et al., 2004; Burgin et al., 2013; Eagles et al., 2012, 2013; García-Lastra et al., 2012; Jacquet et al., 2016; Klausner et al., 2018; McGrath et al., 2018; Rajko-Nenow et al., 2020). The 12-h flight duration limit has been supported by some unpublished data reporting that *Culicoides* midges tethered in a blowing tunnel were able to fly for 10 consecutive hours (Burgin et al., 2013). For most of the other studies, no clear justification of the threshold was provided. In absence of robust experimental evidence, we originally judged plausible the two scenarios (and even a third one considering 72 h of dispersion but not presented here) in accordance with vector lifespan, ranging between 14 days and 6 weeks. These long durations of flight were also the opportunity to test extreme situations, the most conservative in terms of risk. However, when considering results from 48-h maximal flight duration scenario (Supplement Information S3), the insects could be air-borne from Sardinia over an extensive and maximal distance up to 2000 km at a probability lower than $10^{-4}$. Although this probability was very low, such a long distance did not appear consistent with past BTV introductions in Europe. Moreover, no marked difference was observed between 48- and 72-h maximal flight duration scenarios, and the risk of introduction, most of the time, did not exceed the Mediterranean basin (results of the 72-h scenario not presented in this article). In this context, we believe that the 24-h scenario could be the most realistic for future studies as the predictions were consistent with historical BTV incursions from Sardinia.

Another filtering condition that was used to infer the probability of air-borne dispersion was the temperature of the air-mass transporting the insect. A cut-off was set at 10°C in accordance with the minimal isotherm for *Culicoides* established between 9.7 and 12.5°C (Conte et al., 2003; Sellers & Mellor, 1993). Such strict limit could be challenged as each *Culicoides* species might tolerate different conditions (Murray, 1995; Venter et al., 2019). To our knowledge, no temperature threshold was previously considered when assessing the risk of vector dispersal due to wind. Some other authors used a survival function depending on temperatures (ADS model) (Aguilar-Vega et al., 2019; Fernández-Carrión et al., 2018) or limited simulations height to 1000 m (Eagles et al., 2013). These temperature thresholds, together with the PBL, are important elements to account for the topography. As expected, our results estimated a very low or null dispersion probability in mountainous areas. This may partly explain the very low dispersal probability observed in Corsica despite its proximity to Sardinia (Figure 4).

Regarding the HYSPLIT simulations, we favored the use of the trajectory mode rather than the concentration mode (prediction of a particle density deposited at destination knowing a fix number released at the source), unlike other applications that used it alone or in combination of both (Aguilar-Vega et al., 2019; Durr et al., 2017; Eagles et al., 2013; McGrath et al., 2018). We considered that flying vectors could not behave like gases or inert particles, characterized by their physical parameters (such as density or drag
coefficients) that determine their likelihood to be deposited but could maintain themselves in the air, even during passive flight with minimal energy. In addition, the number of particles released at the source, which is fixed in the concentration mode, was found to be highly variable in literature (2500–100,000) (Eagles et al., 2013; McGrath et al., 2018) and not supported by empirical data. Finally, our approach aimed to limit complexity and sources of uncertainty according to the principle of parsimony.

The probability for the *Culicoides* spp. to be uplifted in a windstream was identified as a critical but highly uncertain parameter (Figure 8). In absence of any experimental study in this field from our knowledge, we took a constant and conservative value of 25% of the weekly infected vectors that could enter an ascendant windstream up to 50 m a.g.l. (minimal altitude for HYSPLIT to efficiently simulate atmospheric trajectories) and be subsequently windborne transported. However, this estimate is much higher than the range $10^{-5}$–$10^{-3}$ per day (7.00·10$^{-5}$–6.97·10$^{-3}$ per week) estimated by Hall et al. (2022), irrespective of the insect type among midges, mosquitoes, or heavy flier like tabanids. In contrast to this study, we restricted the simulation start time window to dusk and dawn, which are more suitable times for take-off with ascendant airstream than turbulent daytime (Reynolds et al., 2006). *Culicoides* midges, due to its very small size (1–3 mm), are more likely to be windborne transported and found at high altitude (2 km) (Chapman et al., 2004; Reynolds et al., 2006; Sanders et al., 2011) than mosquitoes or big flies. We judged then our estimate plausible but uncertain and included a large range of variation in the uncertainty assessment. Despite this conservative estimate of $U_i$, the predicted at-risk locations remain limited.

When calculating the BTV basic reproduction number in such large geographical scale as Europe (step 3), discrepancies of spatial distributions and abundances between vector species need to be considered. In this case, the two-host, two-vector $R_0$ formula used in our model (Turner et al., 2013) is believed to provide a better estimate of disease transmission potential than the two-host one-vector formula (Gubbins et al., 2008). However, in our model, *C. imicola* contributed very little to final $R_0$, which might explain lower $R_0$ values in southern European regions than usually calculated (Guis et al., 2012). In Spain, our results showed high risk of BTV transmission not only at the northern border, similarly illustrated by Aguilar-Vega et al. (2020), but also at the northeastern Mediterranean border close to France. In the Netherlands, our $R_0$ values were lower and more homogeneous than illustrated by Hartemink et al. (2009) for Obsoletus complex, probably because we used a lower resolution grid. In France, surprisingly, areas of high BTV transmission were evidenced around Paris and Bordeaux cities where past BTV outbreaks were rare. The areas were characterized by low densities of ruminant hosts but high densities of Obsoletus complex midges leading to high value of vector/host ratios. It is important to note that this $R_0$ approach considers the host population as fully susceptible,

which might not be the case in countries where a vaccination program is implemented or previous outbreaks have occurred. Another improvement of the $R_0$ could be to consider additional vector species in the formula as long as they are shown competent to transmit BTV, such as *Culicoides newsteadi* (Foxi et al., 2019; Kundlacz et al., 2019).

Seasonal variations of $R_0$ estimates were observed during the study period (Supplement Information S3), principally due to the weekly vector abundance as discussed previously, but also due to the vector transmission parameters (biting rate, transmission rate from host to vector, mortality rates...), which are all functions of ambient temperature (Haider et al., 2019; Mullens et al., 2004; Turner et al., 2013; White et al., 2017), considering the poikilothermic nature of the arthropod. These functions were mostly retrieved from studies involving the American species *Culicoides sonorensis* and not Obsoletus complex, due to the difficulty of the latter to be bred in laboratory conditions. Nevertheless, functions used in the model led to higher biting rate, higher vector-to-host transmission, and lower mortality rate for Obsoletus complex midges as compared to *C. imicola*, in line with literature reporting a better vectorial capacity, higher Minimum Infection Rate (Foxi et al., 2016), increased longevity (Meiswinkel et al., 2008), and lower isotherm (Versteirt et al., 2017). Further experiments investigating specifically Obsoletus complex midges would be beneficial to reduce uncertainty around the function's parameters.

Globally, our model computes a high number of parameters that can have a significant impact on the final risk outputs (Figure 8, Supplement Information S2.2). Among the most important is the probability to be uplifted in the air-mass $U_i$ and the parameters specific to the Obsoletus complex (abundance, competence for the BTV-3 strain, and activity rate), which were probably overestimated with respect to the assumptions made, the literature available, and the data used in the model. Even under these conservative conditions, our results suggest that the risk of BTV-3 introduction outside Sardinia remains limited to the southern Mediterranean Basin. The destinations identified at-risk remain consistent when accounting for uncertainty in most influential parameters, but the number of weeks a destination was exposed to a risk >1 was reduced, as was the total risk period between the end of May and the end of October (Figure 9).

As mentioned above, most of the parameters of our model were not serotype specific. Therefore, the location of the source determines the specificity of the model application, and here the southwestern part of Sardinia is the main factor that makes this assessment specific to BTV-3. On the other hand, this model constitutes a generic framework that can be easily reused for other BTV serotypes present in different locations, provided that HYSPLIT simulations are available from that source. Adaptation of the model to other *Culicoides*-borne diseases, such as Epizootic Hemorrhagic Disease Virus, can also be envisioned, provided that the transmission pattern and competent vectors can be considered similar to BTV.
5 | CONCLUSION

In this study, we assessed the risk of BTV introduction by long-distance dispersal of infected vectors considering meteorological data and vector biological limits during flight. The at-risk destinations identified by the model applied for BTV-3 are consistent with historical data for other serotypes and provide valuable information that may facilitate preparedness and target surveillance. Further experimental investigations on *Culicoides* flight conditions and Obsoletus complex-specific parameters are needed to improve the robustness of our predictions and continue efforts in preventing arbovirus emergence.

ACKNOWLEDGMENTS

AB gratefully acknowledges the support received from the French Government through the “Plan de Relance” initiative (ANR-21-PRR-0075-01). TP would like to thank the French National Research Agency and Boehringer Ingelheim Animal Health France for support through the IDEXLYON project (ANR-16-IDEX-0005) and the Industrial Chair in Veterinary Public Health, as part of the VPH Hub in Lyon. HYSPLIT data were produced by DM within a mutual agreement between VetAgro Sup and INRAE. KCM declares no commercial or financial relationships that could be construed as a potential conflict of interest. We also thank Thomas Balenghien for his inputs in the development of our model framework.

DATA AVAILABILITY STATEMENT

Raw data regarding HYSPLIT atmospheric simulations are openly available in data INRAE repository within the “Experimental—Observation—Simulation” dataverse, with the following reference link: https://doi.org/10.57745/IS51BX. Other data can be available upon request to the corresponding author.

ORCID

Amandine Bibard https://orcid.org/0009-0006-3103-9110

Davide Martinetti https://orcid.org/0000-0003-2047-1793

Albert Picado https://orcid.org/0000-0001-7344-628X

Karine Chalvet-Monfray https://orcid.org/0000-0003-0904-4539

Thibaud Porphyre https://orcid.org/0000-0002-8552-9899

REFERENCES

Aguilar-Vega, C., Bosch, J., Fernández-Carrión, E., Lucientes, J., & Sánchez-Vizcaíno, J. M. (2020). Identifying Spanish areas at more risk of monthly BTV transmission with a basic reproduction number approach. *Viruses*, *12*(10), 1158. https://doi.org/10.3390/v12101158

Aguilar-Vega, C., Fernández-Carrión, E., & Sánchez-Vizcaíno, J. M. (2019). The possible route of introduction of bluetongue virus serotype 3 into Sicily by windborne transportation of infected *Culicoides* spp. *Trans-boundary and Emerging Diseases*, *66*(4), 1665–1673. https://doi.org/10.1111/tbed.13201

Alba, A., Casal, J., & Domingo, M. (2004). Possible introduction of bluetongue into the Balearic Islands, Spain, in 2000, via air streams. *Veterinary Record*, *155*(15), 460–461. https://doi.org/10.1136/vr.155.15.460

Balenghien, T., Alexander, N., Arnpórsódtír, A. L., Bisia, M., Blackwell, A., Bodker, R., Bourquia, M., Boutsini, S., Carpenter, S., Colenutt, C., Culverwell, L., Cvetkovikj, A., Dascălu, L., Regge, N. D., Dhollander, S., Elbers, A., England, M., Filatov, S., Garros, C., ..., Wint, W. G. R., (2020). VectorNet Data Series 3: Culicoides Abundance Distribution Models for Europe and Surrounding Regions. Open Heal Data 7:2, 7. https://doi.org/10.5334/ohd.33

Baston, D. (2022). exactextractr: Fast extraction from raster datasets using *polygons*. The Comprehensive R Archive Network. https://CRAN.R-project.org/package=exactextractr

Braks, M., Mancini, G., & Goffredo, M. (2017). Risk of vector-borne diseases for the EU: Entomological aspects—Part 1. *EFSA Supporting Publications*, *14*(2), 1173E. 10.2903/sp.efsa.2017.en-1173

Braverman, Y., & Chechik, F. (1996). Air streams and the introduction of animal diseases borne on *Culicoides* (*Diptera, Ceratopogonidae*) into Israel: -EN- -FR- -ES-. *Revue Scientifique et Technique de l OIE*, *15*(3), 1037-1052. 10.20506/rst.15.3.968

Burgin, L. E., Gloster, J., Sanders, C., Mellor, P. S., Gubbins, S., & Carpenter, S. (2013). Investigating incursions of bluetongue virus using a model of long-distance *Culicoides* biting midge dispersal. *Transboundary and Emerging Diseases*, *60*(3), 263-272. https://doi.org/10.1111/j.1865-1682.2012.01345.x

Calistri, P., Giovannini, A., Conte, A., Nannini, D., Santucci, U., Patta, C., Rolesu, S., & Caporale, V. (2004). Bluetongue in Italy: Part I. *Veterinaria Italiana*, *40*(3), 243-251.

Cappai, S., Rolesu, S., Loi, F., Liciardi, M., Leone, A., Marcacci, M., Teodori, L., Mangone, I., Sghaier, S., Portanti, O., Savini, G., & Lorusso, A. (2019). Western bluetongue virus serotype 3 in Sardinia, diagnosis and characterization. *Transboundary and Emerging Diseases*, *66*(3), 1426-1431. https://doi.org/10.1111/tbed.13156

Carnell, R. (2022). *Ihs: Latin hypercube samples*. R package version 1.1.6. https://CRAN.R-project.org/package=Ihs

Chapman, J. W., Reynolds, D. R., Smith, A. D., Smith, E. T., & Woiwod, I. P. (2004). An aerial netting study of insects migrating at high altitude over England. *Bulletin of Entomological Research*, *94*(2), 123-136. https://doi.org/10.1079/ber2004287

Conte, A., Gilbert, M., & Goffredo, M. (2009). Eight years of entomological surveillance in Italy show no evidence of *Culicoides imicola* geographical range expansion. *Journal of Applied Ecology*, *46*(6), 1332-1339. https://doi.org/10.1111/j.1365-2664.2009.01723.x

Conte, A., Giovannini, A., Savini, L., Goffredo, M., Calistri, P., & Meiswinkel, R. (2003). The effect of climate on the presence of *Culicoides imicola* in Italy. *Journal of Veterinary Medicine, Series B*, *50*(3), 139-147. https://doi.org/10.1046/j.1439-0450.2003.00632.x

Diego, A. C. P., Sánchez-Cordón, P. J., & Sánchez-Vizcaíno, J. M. (2014). Bluetongue in Spain: From the first outbreak to 2012. *Transboundary and Emerging Diseases*, *61*(6), e1-e11. https://doi.org/10.1111/tbed.12068

Dórea, F. C., Swanenburg, M., Roermund, H., Horigan, V., Vos, C. J., Gale, P., Lilja, T., Comin, A., Bahuon, C., Zientara, S., Young, B., Vial, F., Kosmider, R., & Lindberg, A. (2017). Data collection for risk assessments on animal health (acronym: DACRAH): Final report. *EFSA Supporting Publications*, *14*(1), 1171E. https://doi.org/10.2903/sp.efsa.2017.en-1171

Driessche, P. V. D., & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. *Mathematical Biosciences*, *180*(1-2), 29-48. https://doi.org/10.1016/s0025-5564(02)00108-6

Ducheyne, E., Lange, M., Stede, Y. V. D., Meroc, E., Durand, B., & Hendrickx, G. (2011). A stochastic predictive model for the natural spread of bluetongue. *Preventive Veterinary Medicine*, *99*(1), 48-59. https://doi.org/10.1016/j.prevetmed.2011.01.003

Durr, P. A., Graham, K., & Klinken, R. D. V. (2017). Sellers’ revisited: A big data reassessment of historical outbreaks of bluetongue and African horse sickness due to the long-distance wind dispersion of *Culicoides* midges. *Frontiers in Veterinary Science*, *4*, 98. https://doi.org/10.3389/fvets.2017.00098

Eagles, D., Deveson, T., Walker, P. J., Zalucki, M. P., & Durr, P. (2012). Evaluation of long-distance dispersal of *Culicoides* midges into northern
Australia using a migration model. *Medical and Veterinary Entomology*, 26(3), 334–340. https://doi.org/10.1111/j.1365-2915.2011.01005.x

Eagles, D., Melville, L., Weir, R., Davis, S., Bellis, G., Zalucki, M. P., Walker, P. J., & Durr, P. A. (2014). Long-distance aerial dispersal modelling of *Culicoides* biting midges: Case studies of incursions into Australia. *BMC Veterinary Research, 10*, 135. https://doi.org/10.1186/1746-6148-10-135

Eagles, D., Walker, P. J., Zalucki, M. P., & Durr, P. A. (2013). Modelling spatio-temporal patterns of long-distance *Culicoides* dispersal into northern Australia. *Preventive Veterinary Medicine*, 110(3–4), 312–322. https://doi.org/10.1016/j.prevetmed.2013.02.022

Elbers, A. R. W., Heuvel, S. V. D., & Meiswinkel, R. (2016). Diel activity and preferred landing sites in *Culicoides* biting midges attacking fjord horses. *Entomologia Experimentalis et Applicata*, 160(3), 272–280. https://doi.org/10.1111/eea.12481

Fernández-Carrión, E., Ivorra, B., Ramos, Á. M., Martínez-López, B., Aguilar-Vega, C., & Sánchez-Vizcaíno, J. M. (2018). An advection-deposition-survival model to assess the risk of introduction of vector-borne diseases through the wind: Application to bluetongue outbreaks in Spain. *PLoS ONE*, 13(3), e0194573. https://doi.org/10.1371/journal.pone.0194573

Foxi, C., & Delrio, G. (2010). Larval habitats and seasonal abundance of *Culicoides* biting midges found in association with sheep in northern Sardinia. *Medical and Veterinary Entomology*, 24(2), 199–209. https://doi.org/10.1111/j.1365-2915.2010.00861.x

Foxi, C., Delrio, G., Falchi, G., Marche, M. G., Satta, G., & Ruiu, L. (2016). Role of different *Culicoides* vectors (*Diptera: Ceratopogonidae*) in bluetongue virus transmission and overwintering in Sardinia (Italy). *Parasites & Vectors*, 9(1), 440. https://doi.org/10.1186/s13071-016-1733-9

Foxi, C., Meloni, G., Puggioni, G., Manunta, D., Rocchigiani, A., Vento, L., Cabras, P., & Satta, G. (2019). Bluetongue virus detection in new *Culicoides* species in Sardinia, Italy. *Veterinary Record*, 184(20), 621. https://doi.org/10.1136/vr.105118

Gale, P., Kelly, L., & Snary, E. L. (2015). Pathways for entry of live-stock arboviruses into Great Britain: Assessing the strength of evidence. *Transboundary and Emerging Diseases*, 62(2), 115–123. https://doi.org/10.1111/tbed.12317

García-Lastra, R., Leginagoikoa, I., Plazaola, J. M., Ocabo, B., Aduriz, G., Nunes, T., & Juste, R. A. (2012). Bluetongue virus serotype 1 outbreak in the Basque Country (Northern Spain) 2007–2008. Data support a primary vector windborne transport. *PLoS ONE*, 7(3), e34421. https://doi.org/10.1371/journal.pone.0034421

Gilbert, M., Nicolas, G., Cinardi, G., Boeckel, T. P. V., Vanwambeke, S. O., Wint, G. R. W., & Robinson, T. P. (2018). Global distribution data for cattle, buffaloes, horses, sheep, goats, pigs, chickens and ducks in 2010. *Scientific Data*, 5(1), 180227. https://doi.org/10.1038/sdata.2018.227

Goffredo, M., Catalani, M., Federici, V., Portanti, O., Marini, V., Mancini, G., Quaglia, M., Santilli, A., Teodori, L., & Savini, G. (2015). Vector species of *Culicoides* midges implicated in the 2012–2014 bluetongue epidemics in Italy. *Veterinaria Italiana*, 51(2), 131–138. https://doi.org/10.12834/vetit.771.3854.1

Golender, N., Bumbarov, V., Eldar, A., Lorusso, A., Kenigswald, G., Varsano, J. S., David, D., Schainin, S., Dagoni, I., Gur, I., Kaplan, A., Gorohov, A., Koren, O., Oron, E., Khinich, Y., Sclamovich, I., Meir, A., & Savini, G. (2020). Bluetongue serotype 3 in Israel 2013–2018: Clinical manifestations of the disease and molecular characterization of Israeli strains. *Frontiers in Veterinary Science*, 7, 112. https://doi.org/10.3389/fvets.2020.00112

Gómez-Tejedor, C. (2004). Brief overview of the bluetongue situation in Mediterranean Europe, 1998–2004. *Veterinaria Italiana*, 40(3), 57–60.

Gubbins, S., Carpenter, S., Baylis, M., Wood, J. L. N., & Mellor, P. S. (2008). Assessing the risk of bluetongue to UK livestock: Uncertainty and sensitivity analyses of a temperature-dependent model for the basic reproduction number. *Journal of the Royal Society Interface*, 5(20), 363–371. https://doi.org/10.1098/rsif.2007.1110

Gubbins, S., Szmaragd, C., Burgin, L., Wilson, A., Volkova, V., Gloster, J., & Gunn, G. J. (2010). Assessing the consequences of an incursion of a vector-borne disease I. Identifying feasible incursion scenarios for bluetongue in Scotland. *Epidemics*, 2(3), 148–154. https://doi.org/10.1016/j.epidem.2010.05.001

Guis, H., Caminade, C., Calvete, C., Morse, A. P., Tran, A., & Baylis, M. (2012). Modelling the effects of past and future climate on the risk of bluetongue emergence in Europe. *Journal of the Royal Society Interface*, 9(67), 339–350. https://doi.org/10.1098/rsif.2011.0255

Haider, N., Kjær, L. J., Skovgård, H., Nielsen, S. A., & Bødker, R. (2019). Quantifying the potential for bluetongue virus transmission in Danish cattle farms. *Scientific Reports*, 9(1), 13466. https://doi.org/10.1038/s41598-019-49866-8

Hall, D. R., Torpy, J., Nye, R., Emma, & Cowled, B. (2022). *Quantitative risk assessment for the introduction of lumpy skin disease virus into Australia via non-regulated pathways*. CABI Compendium, https://doi.org/10.1079/cabicompendium.76779

Hartemink, N. A., Purse, B. V., Meiswinkel, R., Brown, H. E., Koeijer, A. D., Elbers, A. R. W., Boender, G. J., Rogers, D. J., & Heesterbeek, J. A. P. (2009). Mapping the basic reproduction number (R0) for vector-borne diseases: A case study on bluetongue virus. *Epidemics*, 1(3), 153–161. https://doi.org/10.1016/j.epidem.2009.05.004

Hartley, M., Voller, F., Murray, T., & Roberts, H. (2013). Qualitative veterinary risk assessment of the role of wild deer in the likelihood of incursion and the impact on effective disease control of selected exotic notifiable diseases in England. *European Journal of Wildlife Research*, 59(2), 257–270. https://doi.org/10.1007/s10344-012-0674-7

Hernangómez, D. (2023). *tidyterra: Tidyverse methods and ggplot2 helpers for terra objects*. The Comprehensive R Archive Network. https://CRAN.R-project.org/package=tidyterra

Hijmans, R. (2023). *terra: Spatial data analysis*. R package version 1.7-3. https://CRAN.R-project.org/package=terra

Hoar, B. R., Carpenter, T. E., Singer, R. S., & Gardner, I. A. (2004). Probability of introduction of exotic strains of bluetongue virus into the US and into California through importation of infected cattle. *Preventive Veterinary Medicine*, 66(1–4), 79–91. https://doi.org/10.1016/j.prevetmed.2004.08.006

Holwerda, M., Santman-Berends, I. M. G. A., Harders, F., Engelsma, M. Y., Vloet, R. P. M., Dijkstra, E., van Gennip, R. G. P., Mars, J., Spierenburg, M., Roos, L., van denBrom, R., & van Rijn, P. A. (2023). Emergence of bluetongue virus serotype 3 in the Netherlands in September 2023. *bioRxiv*, https://doi.org/10.1101/2023.09.29.560138

Iooss, B., Veiga, S. D., Janon, A., Pujol, G., Broto, B., Boumhaout, K., Delage, T., Amri, R. E., Fruth, J., Gilquin, L., Guillaume, J., Herin, M., Idrissi, M. I., Gratiet, L. L., Lemaitre, P., Marrel, A., Meynaoui, A., Nelson, B. L., Monari, F., ... Weber, F. (2022). sensitivity: Global Sensitivity Analysis of Model Outputs. https://CRAN.R-project.org/package=sensitivity

Istituto Zooprofilattico Sperimentale della Sardegna. (2022). *Osservatorio epidemiologico veterinario regionale*. Istituto Zooprofilattico Sperimentale della Sardegna. https://www.izs-sardegna.it/CdSE-OEVR-blue-tongue.cfm

Jacquet, S., Huber, K., Pagès, N., Talavera, S., Burgin, L. E., Carpenter, S., Sanders, C., Dicko, A. H., Djerbal, M., Goffredo, M., Lhor, Y., Lucientes, J., Miranda-Chueca, M. A., Fonseca, I. P. D., Ramilo, D. W., Setier-Rio, M.-L., Bouyer, J., Chevillon, C., Balenghien, T., ... Garros, C. (2016). Range expansion of the Bluetongue vector, *Culicoides* imicola, in continental France likely due to rare wind-transport events. *Scientific Reports*, 6(1), 27247. https://doi.org/10.1038/srep27247

Kirkeby, C., Græsbøll, K., Stockmarr, A., Christiansen, L. E., & Bødker, R. (2013). The range of attraction for light traps catching *Culicoides* biting midges (*Diptera: Ceratopogonidae*). *Parasites & Vectors*, 6(1), 67. https://doi.org/10.1186/1756-3305-6-67

Klausner, Z., Klement, E., & Fattal, E. (2018). Source–receptor probability of atmospheric long-distance dispersal of viruses to Israel from the eastern Mediterranean area. *Transboundary and Emerging Diseases*, 65(1), 205–212. https://doi.org/10.1111/tbed.12649

Kundlacz, C., Caignard, G., Sailleau, C., Viarouge, C., Postic, L., Vitour, D., Zientara, S., & Breard, E. (2019). Bluetongue virus in France: An illustration of the European and Mediterranean context since the 2000s. *Viruses*, 11(7), 672. https://doi.org/10.3390/v11070672
Maclachlan, N. J. (2011). *Bluetongue: History, global epidemiology, and pathogenesis*. *Preventive Veterinary Medicine*, **102**(2), 107–111. https://doi.org/10.1016/j.prevetmed.2011.04.005

Massicotte, P., & South, A. (2023). *rnaturalearth: World map data from natural earth*. rOpensci. https://docs.ropensci.org/rnaturalearth/

Maurella, C., Mastrantonio, G., Bertolini, S., Crescio, M. I., Ingravalle, F., Adkin, A., Simons, R. R. L., De Nardi, M., Peña, A. E., Horigan, V., & Ru, G. (2019). Social network analysis and risk assessment: An example of introducing an exotic animal disease in Italy. *Microbial Risk Analysis*, **13**, 100074. https://doi.org/10.1016/j.mran.2019.04.001

McGrath, G., More, S. J., & O'Neill, R. (2018). Hypothetical route of the introduction of Schmallenberg virus into Ireland using two complementary analyses. *Veterinary Record*, **182**(8), 226. https://doi.org/10.1136/vr.104302

Meiswinkel, R., Baldet, T., Deken, R. D., Takken, W., Delécolle, J.-C., & Mellor, P. S. (2008). The 2006 outbreak of bluetongue in northern Europe—The entomological perspective. *Preventive Veterinary Medicine*, **87**(1–2), 55–63. https://doi.org/10.1016/j.prevetmed.2008.06.005

Meiswinkel, R., & Elbers, A. R. W. (2016). The dying of the light: Crepuscular activity in *Culicoides* and impact on light trap efficacy at temperate latitudes. *Medical and Veterinary Entomology*, **30**(1), 53–63. https://doi.org/10.1111/mve.12150

Mintiens, K., Méroc, E., Mellor, P. S., Staubach, C., Gerbier, G., Elbers, A. R. W., Hendrickx, G., & Clercq, K. D. (2008). Possible routes of introduction of bluetongue virus serotype 8 into the epicentre of the 2006 epidemic in north-western Europe. *Preventive Veterinary Medicine*, **87**(1–2), 131–144. https://doi.org/10.1016/j.prevetmed.2008.06.011

Mullens, B. A., Gerry, A. C., Lysyk, T. J., & Schmidtmann, E. T. (2004). Environmental effects on vector competence and virogenesis of bluetongue virus in *Culicoides*: Interpreting laboratory data in a field context. *Veterinaria Italiana*, **40**(3), 160–166.

Murray, M. D. (1995). Influences of vector biology on transmission of arboviruses and outbreaks of disease: The *Culicoides* brevittarsis model. *Veterinary Microbiology*, **46**(1–3), 91–99. https://doi.org/10.1016/0378-1135(95)00075-1

Napp, S., Allepuz, A., García-Bocanegra, I., Alba, A., Vilar, M. J., & Casal, J. (2011). Quantitative assessment of the probability of bluetongue virus transmission by bovine semen and effectiveness of preventive measures. *Theriogenology*, **75**(5), 920–932. https://doi.org/10.1016/j.theriogenology.2010.10.035

Napp, S., García-Bocanegra, I., Pagès, N., Allepuz, A., Alba, A., & Casal, J. (2013). Assessment of the risk of a bluetongue outbreak in Europe caused by *Culicoides* midges introduced through intracontinental transport and trade networks. *Medical and Veterinary Entomology*, **27**(1), 19–28. https://doi.org/10.1111/j.1365-2915.2012.01016.x

Nelson, E., Thurston, W., Pearce-Kelly, P., Jenkins, H., Cameron, M., Carpenter, S., Guthrie, A., & England, M. (2022). A qualitative risk assessment for bluetongue disease and African horse sickness: The risk of entry and exposure at a UK zoo. *Viruses*, **14**(3), 502. https://doi.org/10.3390/v14030502

Pebesma, E. (2022). *sf: Simple features for R*. The Comprehensive R Archive Network. https://CRAN.R-project.org/package=sf

Pouillot, R. (2022). *mc2d: Tools for two-dimensional Monte-Carlo simulations*. The Comprehensive R Archive Network. https://CRAN.R-project.org/package=mc2d

Purse, B. V., Brown, H. E., Harrup, L., Mertens, P. P. C., & Rogers, D. J. (2008). Invasion of bluetongue and other orbivirus infections into Europe: The role of biological and climatic processes. *Revue Scientifique et Technique (International Office of Epizootics)*, **27**(2), 427–442.

Quaglia, M., Foxi, C., Satta, G., Puggioni, G., Bechere, R., Ascentis, M. D., D’Alessio, S. G., Spedicato, M., Leone, A., Pisciella, M., Portanti, O., Teodori, L., Di Gialleonardo, L., Cammà, C., Savini, G., & Goffredo, M. (2023). *Culicoides* species responsible for the transmission of epizootic haemorrhagic disease virus (EHDV) serotype 8 in Italy. *Veterinaria Italiana*, **59**, 83–89. https://doi.org/10.12834/vetit.3347.22208.1

R Core Team. (2022). *R: A language and environment for statistical computing (version v4.2.2)*. R Foundation for Statistical Computing. https://www.R-project.org/

Rajko-Nenow, P., Christodoulou, V., Thurston, W., Ropiak, H. M., Savva, S., Brown, H., Qureshi, M., Alvanitopoulos, K., Gubbins, S., Flannery, J., & Batten, C. (2020). Origin of bluetongue virus serotype 8 outbreak in Cyprus, September 2016. *Viruses*, **12**(1), 96. https://doi.org/10.3390/v12010096

Reynolds, D. R., Chapman, J. W., & Harrington, R. (2006). The migration of insect vectors of plant and animal viruses. *Advances in Virus Research*, **67**, 453–517. https://doi.org/10.1016/s0065-3527(06)67012-7

Ries, C., Vögtlin, A., Hüssy, D., Jandt, T., Gobet, H., Hilbe, M., Burgener, C., Schweizer, L., Häfliger-Speiser, S., Beer, M., & Hoffmann, B. (2021). Putative novel atypical BTV serotype ‘36’ identified in small ruminants in Switzerland. *Viruses*, **13**(5), 721. https://doi.org/10.3390/v13050721

Roberts, D. H., Lucas, M. H., & Bell, R. A. (1993). Animal and animal product importation and the assessment of risk from bluetongue and other ruminant orbiviruses. *British Veterinary Journal*, **149**(1), 87–99. https://doi.org/10.1016/s0007-1935(05)80212-x

Sagués, E. M., Fernández-Carrión, E., & Sánchez-Vizcaíno, J. M. (2019). Risk of introduction of infectious animal diseases for Europe based on the health situation of North Africa and the Arabian Peninsula. *Frontiers in Veterinary Science*, **6**, 293. https://doi.org/10.3389/fvets.2019.00293

Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., & Tarantola, S. (2010). Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. *Computer Physics Communications*, **181**(2), 259–270. https://doi.org/10.1016/j.cpc.2009.09.018

Saminathan, M., Singh, K. P., Khorajiya, J. H., Dinesh, M., Gupta, V. K., S., Maity, M., Rahman, A. F., Misri, J., Malik, Y. S., Gupta, V. K., Singh, R. K., & Dhama, K. (2020). An updated review on bluetongue virus: Epidemiology, pathobiology, and advances in diagnosis and control with special reference to India. *Veterinary Quarterly*, **40**(1), 258–321. https://doi.org/10.1080/01652176.2020.1831708

Sanders, C. J., Selby, R., Carpenter, S., & Reynolds, D. R. (2011). High-altitude flight of *Culicoides* biting midges. *Veterinary Record*, **169**(8), 208. https://doi.org/10.1136/vr.d4245

Savini, G., Goffredo, M., Monaco, F., Gennaro, A. D., Cafiero, M. A., Baldi, L., de Santis, P., Meiswinkel, R., & Caporale, V. (2005). Bluetongue virus isolations from midges belonging to the Obsoletus complex (*Culicoides Diptera: Ceratopogonidae*) in Italy. *Veterinary Record*, **157**(5), 133. https://doi.org/10.1136/vr.157.5.133

Sellers, R. F., & Maarouf, A. R. (1991). Possible introduction of epizootic hemorrhagic disease of deer virus (serotype 2) and bluetongue virus (serotype 11) into British Columbia in 1987 and 1988 by infected *Culicoides* carried on the wind. *Canadian Journal of Veterinary Research* = *Revue Canadienne De Recherche Veterinaire*, **55**(4), 367–370.

Sellers, R. F., & Mellor, P. S. (1993). Temperature and the persistence of viruses in *Culicoides* spp. During adverse conditions: -EN- -FR- -ES- *Revue Scientifique et Technique de l OIE*, **12**(3), 733–755. https://doi.org/10.20506/rst.12.3.722

Simons, R. R. I., Horigan, V., Ip, S., Taylor, R. A., Crescio, M. I., Maurella, C., Mastrantonio, G., Bertolini, S., Ru, G., Cook, C., & Adkin, A. (2019). A spatial risk assessment model framework for incursion of exotic animal disease into the European Union Member States. *Microbial Risk Analysis*, **13**, 100075. https://doi.org/10.1016/j.mran.2019.05.001

Sousa, P. M., Barriopedro, D., Ramos, A. M., García-Herrera, R., Espírito-Santo, F., & Trigo, R. M. (2019). Saharan air intrusions as a relevant mechanism for Iberian heatwaves: The record breaking events of August 2018 and June 2019. *Weather and Climate Extremes*, **26**, 100224. https://doi.org/10.1016/j.wace.2019.100224

Stein, A. F., Draxler, R. R., Rolph, G. D., Stunder, B. J. B., Cohen, M. D., & Ngan, F. (2015). NOAA’s HYSPLIT atmospheric transport and dispersion modeling system. *Bulletin of the American Meteorological Society*, **96**(12), 2059–2077. https://doi.org/10.1175/bams-d-14-00110.1

Sutmoller, P., & Wrathall, A. E. (1997). A quantitative assessment of the risk of transmission of foot-and-mouth disease, bluetongue and vesicular stomatitis by embryo transfer in cattle. *Preventive Veterinary Medicine*, **32**(1–2), 111–132. https://doi.org/10.1016/s0167-5877(96)01118-x
Tuel, A., & Eltahir, E. A. B. (2020). Why is the mediterranean a climate change hot spot? *Journal of Climate*, 33(14), 5829–5843. https://doi.org/10.1175/jcli-d-19-0910.1

Turner, J., Bowers, R. G., & Baylis, M. (2013). Two-host, two-vector basic reproduction ratio (R 0) for bluetongue. *PLoS ONE*, 8(1), e53128. https://doi.org/10.1371/journal.pone.0053128

Venter, G. J., Boikanyo, S. N. B., & Beer, C. J. D. (2019). The influence of temperature and humidity on the flight activity of *Culicoides imicola* both under laboratory and field conditions. *Parasites & Vectors*, 12(1), 4. https://doi.org/10.1186/s13071-018-3272-z

Versteirt, V., Balenghien, T., Tack, W., & Wint, W. (2017). A first estimation of *Culicoides imicola* and *Culicoides obsoletus/Culicoides scoticus* seasonality and abundance in Europe. *EFSA Supporting Publications*, 14(2), 1182E. https://doi.org/10.2903/sp.efsa.2017.en-1182

Viennet, E., Garros, C., Gardès, L., Rakotoarivony, I., Allène, X., Lancelot, R., Crochet, D., Moulia, C., Baldet, T., & Balenghien, T. (2013). Host preferences of palaearctic *Culicoides* biting midges: Implications for transmission of orbiviruses. *Medical and Veterinary Entomology*, 27(3), 255–266. https://doi.org/10.1111/j.1365-2915.2012.01042.x

Villard, P., Muñoz, F., Balenghien, T., Baldet, T., Lancelot, R., & Hénaux, V. (2019). Modeling *Culicoides* abundance in mainland France: Implications for surveillance. *Parasites & Vectors*, 12(1), 391. https://doi.org/10.1186/s13071-019-3642-1

Vos, C. D., Hoek, M., Fischer, E., Koeijer, A. D., & Bremmer, J. (2011). Risk assessment framework for emerging vector-borne livestock diseases. Report 11-CVI0168. Lelystad, the Netherlands: Central Veterinary Institute, part of Wageningen UR. Available from: https://edepot.wur.nl/198115

Vos, C. J. D., Hennen, W. H. G. J., Roermund, H. J. W. V., Dhollander, S., Fischer, E. A. J., & Koeijer, A. A. D. (2021). Assessing the introduction risk of vector-borne animal diseases for the Netherlands using MINTRISK: A model for integrated risk assessment. *PLoS ONE*, 16(11), e0259466. https://doi.org/10.1371/journal.pone.0259466

Wageningen Bioveterinary Research (WBVR). (2023). *Bluetongue in ruminants, updates 2023*. Wageningen Bioveterinary Research. https://www.wur.nl/nl/onderzoek-resultaten/onderzoeksinstituten/bioveterinary-research/show-bvr/blauwtong-bij-herkauwers-updates-2023.htm

White, S. M., Sanders, C. J., Shortall, C. R., & Purse, B. V. (2017). Mechanical model for predicting the seasonal abundance of *Culicoides* biting midges and the impacts of insecticide control. *Parasites & Vectors*, 10(1), 162. https://doi.org/10.1186/s13071-017-2097-5

World Organisation for Animal Health (WOAH). (2010). *Handbook on import risk analysis for animals and animal products* (2nd ed., Vol. 1). WOAH.

Yakob, L., & Clements, A. C. A. (2013). A mathematical model of chikungunya dynamics and control: The major epidemic on Réunion Island. *PLoS ONE*, 8(3), e57448. https://doi.org/10.1371/journal.pone.0057448

## SUPPORTING INFORMATION

Additional supporting information can be found online in the Supporting Information section at the end of this article.

**How to cite this article:** Bibard, A., Martinetti, D., Giraud, A., Picado, A., Chalvet-Monfray, K., & Porphyre, T. (2025). Quantitative risk assessment for the introduction of bluetongue virus into mainland Europe by long-distance wind dispersal of *Culicoides* spp.: A case study from Sardinia. *Risk Analysis*, 45, 108–127. https://doi.org/10.1111/risa.14345