The possible route of introduction of bluetongue virus serotype 3 into Sicily by windborne transportation of infected Culicoides spp.

Summary In October 2017, the first outbreak of bluetongue virus serotype 3 (BTV‐3) began in Italy, specifically in western Sicily. The route of entrance remains unclear, although since 2016 the same strain had been circulating only 150 km away, on the Tunisian peninsula of Cape Bon. The present analysis assessed the feasibility that wind could have carried BTV‐3‐infected Culicoides spp. from Tunisia to Sicily. An advection‐deposition‐survival (ADS) model was used to estimate when and where Culicoides spp. were likely to be introduced prior to the first BTV‐3 report in Italy. Additionally, the Hybrid Single‐Particle Lagrangian Integrated Trajectory (HYSPLIT) model was used to support ADS outputs. The modelling suggests that during September 2017, strong wind currents and suitable climatic conditions could have allowed the transportation of Culicoides spp. from BTV‐3‐infected areas in Tunisia into Sicily. ADS simulations suggest that particles could have reached the province of Trapani in western Sicily on 2 and 12 September. These simulations suggest the feasibility of aerial transportation of infected Culicoides spp. from Tunisia into Sicily. They demonstrate the suitability of the ADS model for retrospective studies of long‐range transportation of insects across large water bodies, which may enhance the early detection of vectorial disease introduction in a region.


In November 2016, a BT outbreak was identified in the Beni
Khalled delegation in the Nabeul Governorate in north-eastern Tunisia (OIE, 2017a) ( Figure 1). The causative agent was identified as a new strain of BT serotype 3 (BTV-3) called BTV-3 TUN2016 . This was the first time BTV-3 had ever been reported in the country. BTV-3 RNA and antibodies were found to be widespread in Tunisia at the end of 2016 and beginning of 2017, and a different strain of BTV-3 was detected farther south (Lorusso et al., 2018). In December 2017, Italian authorities reported to the OIE an outbreak of BTV-3 beginning on 26 October in the Trapani municipality in western Sicily (OIE, 2017b) ( Figure 1).
This serotype was not known to be circulating in Italy or any other European country. One crossbred sheep (Ovis aries) in a flock of 443, presented BT-like clinical signs (fever, head oedema, nasal discharge and depression). Samples from the affected sheep were tested at the Istituto Zooprofilattico Sperimentale of Sicily, and found to be positive for BTV RNA and anti-BTV antibodies . Samples were sent to the Italian National Reference Laboratory for confirmation and strain identification, where a seroneutralization test with all reference BTV serotypes gave a positive result for BTV-3  on 24 November 2017 (OIE, 2017b). Furthermore, preliminary sequencing analysis demonstrated high nucleotide identity with strain BTV-3 TUN2016 circulating in northern Tunisia . Fewer than 150 km separate the western coast of Sicily and the eastern coast of Cape Bon peninsula, where the Nabeul Governorate is located. This proximity suggests the possibility of windborne BTV-3 introduction from Tunisia to Sicily.
To determine the feasibility of windborne introduction, this study applied two atmospheric dispersion models to examine longrange aerial transportation of infected Culicoides spp. Recently, our group developed an advection-deposition-survival (ADS) model to assess the risk of introduction of engorged and potentially infected midges in Spain (Fernández-Carrión et al., 2018). This model is updated hourly, making it a potentially powerful tool for surveillance. In parallel, we analysed data using the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (version 4) of the Air Resources Laboratory at the US National Oceanic and Atmospheric Administration (Draxler & Hess, 1998). Our analyses suggest that the ADS model can be used in retrospective studies of disease vector transportation across large water bodies and, therefore, as a potential early warning tool.

| Analysis of windborne Culicoides trajectories and ground deposition
In order to determine whether windborne advection is a feasible route of introduction of BTV-3 in western Sicily, continuous simulations were carried out with ADS and HYSPLIT models to identify possible days and times of day when introduction occurred.
The ADS model was developed for estimating risk of wind-borne introduction of flying insects into a country. It identifies geographic areas and periods where and when risk of vector-borne diseases incursion is high. Concretely, the model predicts the number density of introduced small flying insects such as Culicoides spp. over space and time based on wind advection, vertical deposition rate and insect survival rate (Fernández-Carrión et al., 2018). The model assumes an initial number density of midges in the source territory (in this case, the Cape Bon peninsula in Tunisia), based on the Culicoides spp. activity function, which correlates with local temperature. Then the model simulates the long-range transportation of midges using Lagrangian particle tracking based on wind direction and intensity until their deposition in the target territory (in this case, Italy and Malta). The model takes into account specific physical and aerodynamic properties of the transported midges (i.e., midge density, Reynolds number and drag coefficient). It also computes a survival rate, which depends on temperature and reflects the probability that an insect may survive transportation at high altitudes, where temperatures can be extreme (Fernández-Carrión et al., 2018). A modification was implemented for the estimation of Culicoides' density (ρ = m/V), for which midge volume (V) was estimated according to the formula: ln(V) = 1.019ln(m dry ) + 1.46, where m dry refers to dry midge weight (Kühsel, Brückner, Schmelzle, Heethoff, & Blüthgen, 2017), which in turn was estimated to be 30% of total weight (Wigglesworth, 1972). Hourly mean temperature (in °C), mean F I G U R E 1 Map showing areas affected by BTV-3 TUN2016 in Tunisia and Sicily. Light grey areas show the administrative areas affected in Tunisia and Sicily. Dark grey areas show the affected Tunisian Delegations where BTV-3 TUN2016 were detected in Nabeul and Monastir Governorates (Lorusso et al., 2018). The Sicilian outbreak of BTV-3 is marked with a black circle horizontal wind direction and wind speed (in m/s) were retrieved from an online climatic repository (www.forec ast.io) for the area and the study period, at temporal resolution of 1 hr (Fernández-Carrión et al., 2018).
The HYSPLIT model (available online at https ://ready.arl.noaa. gov/) is a complex atmospheric system for computing simple air particle trajectories as well as complex transport, dispersion and deposition simulations using archived meteorological data (Stein et al., 2015). The HYSPLIT model is widely used in atmospheric sciences (Stein et al., 2015), and it has even been used to study windborne transportation of Culicoides spp. (Durr, Graham, & van Klinken, 2017;Eagles, Walker, Zalucki, & Durr, 2013;García-Lastra et al., 2012), even though their physical and aerodynamic features differ substantially from those of sand or aerosols. The model relies on meteorological data obtained through the Real-time Environmental

Applications and Display System and archived in the Global Data
Analysis System (GDAS1), with temporal resolution of 3 hr.
In both models, the initial altitude and initial positions of particles can be specified. More detailed information is available elsewhere about the methodology and specifications of the ADS model (Fernández-Carrión et al., 2018) and HYSPLIT model (Draxler & Hess, 1998;Stein et al., 2015).

| Study period and area
BT serological response appears between 7 and 14 days post-infection (OIE, 2014). Since the BTV-3 outbreak in Sicily started on 26 October 2017 (OIE, 2017a), our modelling focused on the period from 1 August to 18 October 2017. To configure the ADS model, we set the extension of the study area to lie between latitudes 34ºN and 39ºN and between longitudes 9.5ºE and 16.5ºE. This area includes north-eastern Tunisia, Sicily, Malta, and the south-western area of mainland Italy ( Figure 2). Within the study area, a grid of 49 x 49 points was created for retrieval of meteorological data. In the ADS and HYSPLIT models, the source of potentially infected midges was set as the Cape Bon peninsula (Tunisia), since BTV-3 was known to be in that area (Lorusso et al., 2018). The area of potential introduction was mainland Italy, the Italian islands and Malta (which fell within the study area).

| Model simulations
A simulation was performed with the ADS model at 1,000 metres above sea level for the study period in order to identify the days and times of highest probability of midge introduction in Sicily from the Cape Bon peninsula. Table 1 shows the major physical and aerodynamic properties of the midges. Once the potential days of introduction were determined, they were studied in greater detail using both the ADS and HYSPLIT models. In the HYSPLIT model, forward  Figure 3 shows the density of particles deposited in Sicily during the study period in the ADS simulation. Slight incursions of particles were predicted to scatter along Sicily, with main deposits on the southern coast and the south-western area of the island. In August, no particle was predicted to land in Sicily. In September, only 4 days of potential arrival of infected midges to target territory were predicted: days 2, 10, 11, and 12. In October, deposition was predicted only on day 6.

| Model outputs
The most probable times of midge incursion and deposition in Sicily on these five dates were predicted using both models (Table 2). No particles were predicted to reach Malta during the study period. Figure 4 shows the density of particles deposited in the target territory on the days when either model predicted that particles were deposited in western Sicily. On 2 September, Trapani province was predicted by both models to accumulate most of the deposition in the study period. On 10 September, the HYSPLIT model predicted deposition of particles in the Trapani region, in contrast to the ADS model. Both models predicted that on 12 September, there was scarce F I G U R E 2 Study domain between latitudes 34ºN and 39ºN, and between longitudes 9.5ºE and 16.5ºE. Dots represent the center of grid cells of 49 x 49 (2,401 in total) used for climatic data interpolation. Grid cell size is 0.1 x 0.14°9 During the dates considered to be at high risk of vector-borne introduction (2, 10, and 12 September), wind directions were predicted to enhance the arrival of midges to the location of the outbreak in Sicily ( Figure 5). On 2 and 12 September, wind speeds at ground level were moderate to low during the entire day, ranging from four to less than 11 m/s on 2 September and, from one to less than 7 m/s on 12 September. On 10 September, winds were stronger and reached speeds of up to 17 m/s during some periods of the day.

The prediction of September as the month when infected
Culicoides spp. were introduced into Trapani correlates with the Culicoides spp. population peak in Tunisia . The average temperature in Sicily in September 2017 was 25 ºC, and adult Obsoletus complex and Culicoides sonorensis midges can survive for up to three months at, respectively, 17-25ºC and 10ºC (Goffredo, Romeo, Monaco, Di Gennaro, & Savini, 2004;Lysyk & Danyk, 2007), while C. sonorensis can live for up to 28 days at 30ºC (Lysyk & Danyk, 2007). Although C. imicola lifespan has been less studied, adults can survive more than 15 days (Nevill, 1971;Paweska, Venter, & Mellor, 2002). Given that the outbreak in Sicily involved BTV-3-naïve animals and did not lead to multiple outbreaks, we hypothesize that a small number of BTV-3-infected midges arrived at one location on the western Sicilian coast and were sufficient to transmit BTV to a susceptible host (Baylis, O'Connell, & Mellor, 2008)  temperature requirements can differ substantially between C. imicola and the Palaearctic species (Wilson & Mellor, 2009), the latter appear to be absent from Tunisia or at least present at very low levels . Thus, the model is suitable for the study in this region.
The ADS and HYSPLIT models showed good agreement (Table 2, Figure 4). They predicted the same prevailing winds and differed primarily in how far into Sicily the particles were deposited (Figure 4).
These differences can be explained in part by the fact that the ADS model takes into account particle size and density (Fernández-Carrión et al., 2018), while the HYSPLIT model does not, since velocity is defined for dry deposition (Durr et al., 2017). In addition, both models rely on different meteorological datasets. The HYSPLIT model uses data from the GDAS1, which includes vertical wind velocity, 1-degree grid resolution and a time window of 3 hr (Stein et al., 2015).
Meteorological data in the ADS model lacks some aspects of vertical wind velocity but features spatial resolution of 1 hr (Fernández-Carrión et al., 2018). We suggest that 1-hr temporal resolution is sufficient to allow precise real-time and retrospective analysis of Culicoides spp. aerial transport without posing excessive computational demands.

| CON CLUS IONS
Two models suggest that BTV-3 could have arrived via windborne transport of infected midges to the Trapani province of Sicily, which may explain the occurrence of BT in a single animal there in October F I G U R E 4 Particle deposition according to ADS and HYSPLIT models for the potential days of introduction. The blue dot shows the location of the Sicilian outbreak of BTV-3. Stars represent the source of midges and particles in the Cape Bon peninsula. Deposition time of the ADS model and the initial and final time of HYSPLIT outputs are gathered in Table 2

CO N FLI C T O F I NTE R E S T S
The authors have no conflict of interest to declare.