A quantitative analysis of the spatial and temporal evolution patterns of the bluetongue virus outbreak in the island of Lesvos, Greece, in 2014

Bluetongue (BT) is an infectious, vector-borne disease of ruminants caused by the bluetongue virus (BTV), a member of the genus Orbivirus which belongs to the family Reoviridae. Twenty-seven serotypes of BTV have so far been isolated and identified (Rojas, Rodríguez-Martín, Martín, & Sevilla, 2019). BTV is transmitted between hosts by bites of Culicoides midges and replicates in all ruminants, but severe disease is restricted mainly to certain breeds of sheep and some deer (Taylor, 1986). BT can be enzootic in areas where the mammalian reservoirs, the virus and the insect vectors have the opportunity to coexist in Received: 6 April 2019 | Revised: 14 January 2020 | Accepted: 11 March 2020 DOI: 10.1111/tbed.13553

climatic conditions conducive to BTV replication and transmission (Caporale et al., 2014). As a result, BTV was historically present exclusively in tropical and subtropical areas of the world, where suitable conditions exist. However, over the last decades, the global distribution of BT has expanded dramatically, potentially due to a variety of factors, including an increased global travel and commerce, deforestation and recent global climatic changes in combination with the related spread of the main BTV vectors of Culicoides imicola (Caporale et al., 2014;Purse et al., 2005;Randolph & Rogers, 2010) and Obsoletus complex (Federici et al., 2019). In addition, recent evidence from California indicated the possibility that the virus can survive through the winter in long-living C. soronensis female midges, which had been infected during the previous flying activity period (Mayo et al., 2014). The latter can be affected by temperature, light intensity, lunar cycles, relative humidity, wind velocity and other weather conditions, with the mean distance flown by female Culicoides spp. being about 2 km (Mellor, 1990).
In Europe, BT occurred for the first time between 1956and 1960in Spain and Portugal (Manso-Ribeiro et al., 1957, while the first epidemic of the disease in Greece occurred on the island of Lesvos in October 1979 and was caused by serotype 4 (Nomikou, Mangana-Vougiouka, & Panagiotatos, 2004). After 1998, BT spread from Turkey to the Greek mainland and northwards and, by 2001, reached Italy (Calistri et al., 2004;Nomikou et al., 2004). Between 2006 and 2008, an extensive epidemic of BTV-8 occurred in Europe . At the end of May 2014, a case of bluetongue virus serotype 4 was clinically diagnosed in Greece, in south Peloponnese, and by the end of August 2014, it was prevalent in the whole country, including the island of Lesvos (World Organisation for Animal Health, 2015). Cases of the disease caused by  were then (July-December 2014) reported in Albania, Croatia, the Republic of Northern Macedonia, Hungary, Montenegro, Romania, Bulgaria, Serbia and Turkey (Niedbalski, 2015;World Organisation for Animal Health, 2015).
Notwithstanding suggestions that vaccinations could have led to limit the disease at the initial stages in Greece (Vasileiou et al., 2016), the authorities had opted to limit the epidemic by controlling and Veterinary Medicine, regional unit of Lesvos (DREVM). The legal definition of a BT case was described as follows: (a) when a ruminant exhibits clinical signs compatible with BT, (b) when a sentinel bovine animal presents seroconversion, (c) when an animal is serologically BTV-positive or BTV has been detected by molecular techniques and (d) when BTV is isolated from an identified animal (Greek Ministry of Rural Development and Food, 2015).
In the current paper, the authors provide descriptive summaries of the spatial and temporal patterns of the available data regarding the 2014 epidemic in Lesvos Island, Greece. Measures of spatial autocorrelation in order to explore potential spatial patterns in disease transmission are also presented. Additionally, in the present study, a spatio-temporal regression-based Bayesian modelling approach is developed in order to characterize the spread of the BT virus between the infected farms of the island. The transmission was modelled via a spatial kernel-based component describing the spread of disease from infected to non-infected farms. Our model includes an Ornstein-Uhlenbeck (OU) process, a continuous-time model used for absorbing model mis-specification, capturing distributional deviations and the inherent serial correlation of the data.

| Study area
The study area was the Greek island of Lesvos, located in the Northeastern Aegean Sea. It has an area of 1,633 km 2 with 320 km of coastline, making it the third largest island in Greece. It is close to Turkey, separated by a narrow Strait (see Figure 1). A total of 165,831 small ruminants had been reported by the DREVM as susceptible, farmed within 617 holdings, while the reported number of bovine was 6,696 animals. In previous years, serotypes  were present in the island (Kyriakis et al., 2015;Nomikou et al., 2004).
On 17 December 2013, the average temperature in Lesvos was 8°C, with 66% humidity and 13 km/hr wind speed as recorded by Mytilene Airport Weather Station, while on 30 May 2014, the average temperature was 21°C, with 68% humidity and 8 km/hr wind speed. In addition, aggregate farm-level data on the number of sheep and goat flocks for each municipality of the island were recorded and utilized for subsequent analyses, especially concerning the spatial epidemic modelling. We assumed in the subsequent analyses that the unit of population is the farm. The aggregated information on the number of susceptible farms in the island for each municipality provided by the DREVM did not contain the exact spatial location for each farm. Therefore, each farm was randomly allocated in space within each municipality using the Random Points tool available in the QGIS software (QGIS Development Team, 2015). In particular, using the command 'Randomly generate points within the polygon', we placed randomly distributed points within the boundaries of the polygons (in this case the boundaries of each municipality in the island). Next, having predefined the sampling of random point locations based on characteristics such as the frequency, that is the field with the number of points of each polygon and the minimum distance allowed, we have sampled the random spatial distribution of the farms within each municipality.

| Purely spatial measures
We examined the characteristics of the temporal and spatial spread of BT during the 2014 epidemic, by using appropriate spatial measures. Specifically, we estimated the Global Moran's I and the Geary's C statistic in order to test for spatial autocorrelation in disease transmission (Lessler, Salje, Grabowski, & Cummings, 2016).
We later employ suitable Bayesian spatio-temporal models that utilize both spatial and temporal information as well as including a continuous-time stochastic Ornstein-Uhlenbeck part that allows modelling the inherent serial correlation of the data-generating mechanism (Malesios, Demiris, Kalogeropoulos, & Ntzoufras, 2017; Oravecz & Tuerlinckx, 2011).

| Global Moran's I
The global Moran's I test (Li, Calder, & Cressie, 2007;Moran, 1950) is a measure of spatial autocorrelation used to test for spatial clustering of a region throughout the study area without the ability of locating specific clusters sites. The global weight matrix is used to define the spatial relationships so that regions close in space are given a greater weight. Moran's I takes values between −1 and +1 with positive values (near +1) indicating the presence of clustering while a negative I value indicates tendency towards dispersion. For the Global I statistic, the null hypothesis states that the attribute being analysed is randomly distributed among the features in the selected study area. The I test statistic is defined as: where n is the number of observations on a variable x at locations i , j, x is the mean of the x variable, w ij are the elements of the weight matrix, and S 0 is the sum of the elements of the weight matrix: In order to identify the location of potentially significant clusters, suitable heat maps are also presented. All spatio-temporal maps for the visualization of the spread of the epidemic were constructed with the use of the QGIS software.

| Geary's C
Geary's C statistic (Geary, 1954) is based on the deviations in responses of each observation with one another and is defined as: Geary's C ranges from 0 (maximal positive autocorrelation) to a positive value for high negative autocorrelation. Note that its expectation is 1 in the absence of autocorrelation and regardless of the specified weight matrix (Sokal & Oden, 1978). (1)

F I G U R E 1 Island of Lesvos, Greece (in grey colour)
Moran's I is more sensitive to extreme values and is thought of as a global measure, whereas the Geary's C is more sensitive to differences between values in neighbouring areas.

| Spatio-temporal modelling
To estimate the BTV occurrence at the farm-level, a spatio-temporal disease transmission model was developed following the Bayesian paradigm. The total number of farms is subdivided in two sets, the farms which were infected during the epidemic and the remaining set of ultimately susceptible farms.

| The model
We now describe the general class of spatio-temporal epidemic models utilized for our analysis. If y i denotes the number of herds infected with BT virus at time t i where i ∈ 0,1,...,132 is ordered chronologically by day for the complete duration of the epidemic period, then it is assumed that the number of infected herds y i at time t i is described by the following model: where g denotes the assumed distribution of the data and i is the rate at which new infections take place during the disease progress.
Equation (5) is used to model the serial correlation in the data, in the form of an Ornstein-Uhlenbeck process (see e.g. Taylor, Cumberland, & Sy, 1994). Here, B t denotes standard Brownian motion and the mean reversion parameter (i.e. the rate at which the process returns to its mean) of the OU process. The piecewise constant deterministic process μ t is given by: with each μ (i) corresponding to t i ≤ t < t i+1 i = 0,1,...,132 given by: In Equation (7), 0 represents the intercept and d k the spatial kernel for the incorporation of spatial information associated with the rate at which infection passes from an infected farm at times i − j (j = 1,2,3,4) to a susceptible farm k at time i . We assumed a 4-day incubation period (Scott, 2011; World Organisation for Animal Health, 2013). The kernel component of the model captures the fact that BT is more likely to spread between farms located nearby than farms located more distant apart: Here, | | d i | | denotes the cardinality of d i , the term d k denotes some specified function of the distance between the infected and susceptible farms based upon their distance d k and d min , which is fixed a priori, restricts the minimum distance over which infections do not occur (Deardon et al., 2010). For our analysis, we have set The functional forms that were tested for are given in Table 1 Regarding the specification of function g in Equation (3) Fat-tailed (C)).

| Model selection
We additionally present the results of selecting the most-suited model by evaluating their performance on capturing disease evolution. Deviance-based measures have been utilized for the comparison and selection among the P, NB, ZIP and ZINB models due to the well-known equivalence in model selection using cross-validation or AIC (see Stone, 1977). Smaller values indicate better fit (Spiegelhalter, Best, Carlin, & Linde, 2002). (3)

| Implementation
A weakly informative N 0,10 4 prior was used for the intercept 0 and the parameters of the kernel component . All models were estimated using the WinBUGS (Lunn, Thomas, Best, & Spiegelhalter, 2000) software. The posterior results were obtained after discarding the initial 5,000 iterations, using an additional sample of 10,000 iterations (using a thinning lag of 10 for minimizing autocorrelation).  Table A1 in the Appendix.

| Descriptive analysis on the 2014 bluetongue epidemic in Lesvos Island, Greece
The spatial distribution and density of sheep and goat farms in the island of Lesvos in 2014 is shown in Figure 3. The majority of sheep and/or goat farms network lies within the northwest part of the island, mainly in the prefectures of Kaloni and Mantamados.

| Spatio-temporal description of the epidemic
In Figure     Moran's I values revealed that the spatial distribution of the BT cases in the Lesvos Island is spatially non-random. Specifically, the null hypothesis that the BTV cases were randomly distributed was rejected for the majority of the time periods, suggesting spatial autocorrelation for BT incidence in the Lesvos Island (Table 2) Mytilene (see Figure 4, map E).
Further, under the null hypothesis of no global spatial autocorrelation, the results of Geary's C coefficients indicate a positive spatial autocorrelation in our data (Table 3).

| Spatio-temporal regression modelling results
Among the 12 candidate models (Table 4), the ZIP class had the best fit to the data in the sense that its posterior mean deviance is the smallest. Further, the best fitted kernel was the Gaussian kernel B with the exponential kernel A being very close. In contrast to the good fit of the two light-tailed kernels, the worst fit was observed for the fat-tailed kernel C across the different distributional assumptions. In summary, these results indicate that BTV is more likely to spread between farms located nearby as opposed to farms more distant apart.
Next, for the best spatio-temporal model (i.e. ZIP model with a Gaussian kernel B), we present the posterior estimates of model parameters (posterior medians and 95% credible intervals) in Table 5. Note that in all the spatio-temporal models, the spatial kernel parameter was significant, suggesting a positive association between new infections and spatial spread of the disease.
For the best model, the distance affects both the occurrence (infection rate part) and non-occurrence of the disease (excess zero part), since the parameter is statistically significant in both componenets.   (Kyriakis et al., 2015;Vasileiou et al., 2016).

| D ISCUSS I ON
In the current article, we present for the first time the data The Bayesian regression models suggest that BTV is more likely to spread between closely located farms. BTV is a vector-borne disease, transmitted by Culicoides biting midges; therefore, their flying range provides an implicit upper bound to the distance of disease spread (Kluiters, Swales, & Baylis, 2015;Szmaragd et al., 2009). Furthermore, in all the models, the spatial kernel parameter was statistically significant for the occurrence of the disease, suggesting that there is a positive association between new infections and spatial spread of the disease. Our selected ZIP model is based upon the Gaussian kernel part, which indicates that during the epidemic the majority of the incidents in the spread of BT were located in distances smaller than 12 km.
Since our best selected model is based upon the Gaussian kernel part, we have strong evidence to suggest that a circle radius around 12 km at the onset of a similar epidemic in the region could both limit economic damage and prove sufficient in restricting the initial epidemic growth.
Various studies have also attempted to model the spatial distri- intermediary stops (Hendrickx et al., 2008) and with the midges being able to fly a mean distance of 2 km (Mellor, 1990). Results from the model indicated that temperature, precipitation, farm density, land area and proportion of pasture are important covariates for describing the BT dynamics.
Supporting our suggestion of the 12-km distance between the major incidents in the spread of BT, Hendrickx et al. (2008) showed that 50% of cases occur within 5 km of the previous case, while 95% of cases occurred within 31 km of it. This happens even though BTV is non-contagious since presence of infection in neighbouring municipalities implies that a vector with the virus might be present, can be transported to the neighbouring municipality and may bite the animal in that municipality, causing the transmission of BTV.
We did not consider the effect of meteorological variables in this epidemic. However, it has been observed (see for example Malesios et al., 2016) that their effect is materially negligible in the presence of acute short-lived epidemics, of the kind considered in the present paper. Therefore, the absence of such variables is not expected to have any material bearing upon the conclusions of this work.
In addition, we noted that in the absence of exact spatial allocation of the farm density of the island, we have utilized a random spatial distribution of the farms within each municipality (average surface of municipalities is 124.32 km 2 ). Our results are, therefore, conditional upon this allocation process. However, each municipality represents a relatively small area and the actual position of the susceptible farms cannot vary substantially. Therefore, the lack of geographic coordinates is not expected to materially affect the results and the fit of the model provides evidence in this direction.
In summary, we investigated the effect of spatial transmission of BTV from infected to non-infected small ruminants' farms in the island of Lesvos, Greece. Spatial analysis showed that disease outbreaks appear clustered due to spatial dependence during the 2014 BTV serotype 4 incursion in Greece. Furthermore, we described the application of spatial Bayesian regression to predict the occurrence of the disease. A zero-inflated Poisson regression model, incorporating serial correlation through an Ornstein-Uhlenbeck process and spatial dependence through a light-tailed kernel part had adequate predictive ability. The simultaneous consideration of the temporal component of disease evolution did not alter the intuitively reasonable conclusion of the absence of long-range transmission for these data. Information derived by our model can be used to anticipate changes in the spatio-temporal distribution of disease and to assist in surveillance efforts. This analysis showed that the epidemic emerged from one location and then transmitted following a non-random (mostly localized) spatial pattern.

ACK N OWLED G M ENTS
The authors would like to thank the 'Directorate of Rural Economy and Veterinary Medicine, regional unit of Lesvos' for the assistance in the data acquisition. N.D. acknowledges support from the Athens University of Economics and Business' Research Centre Action: 'Original Scientific Publications'.

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

E TH I C A L S TATEM ENT
The authors declare that ethical statement is not applicable because sample collection or questionnaires from animals has been gathered.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.