Epidemiology and spatial distribution of bluetongue virus in Xinjiang, China

Bluetongue (BT) is a non-contagious disease affecting domestic and wild ruminants. Outbreaks of BT can cause serious economic losses. To investigate the distribution characteristics of bluetongue virus (BTV), two large-scale censuses of BTV prevalence in Xinjiang, China were collected. Spatial autocorrelation analysis, including global spatial autocorrelation and local spatial autocorrelation, was performed. Risk areas for BTV occurrence in Xinjiang were detected using the presence-only maximum entropy model. The global spatial autocorrelation of BTV distribution in Xinjiang in 2012 showed a random pattern. In contrast, the spatial distribution of BTV from 2014 to 2015 was significantly clustered. The hotspot areas for BTV infection included Balikun County (p < 0.05), Yiwu County (p < 0.05) and Hami City (p < 0.05) in 2012. These three regions were also hotspot areas during 2014 and 2015. Sheep distribution (25.6% contribution), precipitation seasonality (22.1% contribution) and mean diurnal range (16.2% contribution) were identified as the most important predictors for BTV occurrence in Xinjiang. This study demonstrated the presence of high-risk areas for BTV infection in Xinjiang, which can serve as a tool to aid in the development of preventative countermeasures of BT outbreaks.


INTRODUCTION
Bluetongue (BT) is a non-contagious disease affecting ruminant and camelid species (Hofmann et al., 2008). BT, caused by Bluetongue virus (BTV), is a vectorborne disease transmitted between ruminant hosts by blood-feeding midges of the Culicoides spp.
Hosts of BTV infection are domestic and wild ruminants, including sheep, goats, cattle and deer. Among these, goats and cattle are often considered as asymptomatic reservoir hosts (Maclachlan, 1994) or sub-clinically affected (Maclachlan et al., 2009). Severe clinical signs are often seen in certain breeds of sheep, European fine wool and mutton breeds for example (Maclachlan et al., 2009). The most commonly observed clinical signs include fever, hyperemia in nasal and oral mucosa, edema in the lip, ulcers of the oral mucosa, cyanosis of the tongue, and skeletal muscle deformation. Cyanotic tongues are the most obvious characteristic that aid in differentiation from other diseases.
Late in the 18th century, BT was first reported officially in Cape of Good Hope, South Africa. After a systematic clinical study, Spreull named the disease "BT" for the first time in 1905, with reference to the characteristic cyanotic tongues of the infected sheep (Spreull, 1905). In 1943, an outbreak of BT occurred in Cyprus, which is believed to be the first occurring outside of Africa (Gambles, 1949). Since that time, BT has subsequently occurred in many regions of the world, with Antarctica being the only continent free of BTV infection for now (Maclachlan et al., 2009). In China, BTV has become widely distributed throughout the mainland since it was first reported in Yunnan Province in 1979. At present, 11 serotypes of BTV have been isolated in China (Yang et al., 2016;Zhang et al., 1999), with BTV-1 and BTV-16 being the most commonly isolated serotypes (Lee et al., 2011;Zhang et al., 1999).
Outbreaks of BT can create serious economic consequences. Extensive measures are required to control the spread of the virus among infected and cohoused livestock. It is estimated that economic losses resulting from a BT outbreak in 1996 totaled more than $3 billion USD worldwide (Tabachnick, 1996). Thus, BT was included on the World Organisation of Animal Health (OIE) list of notifiable diseases in the mid 1960s, and is classified as a Class A disease of concern in China as well. At present, there are no effective treatments for the disease. Thus, measures to prevent and control BT outbreaks are of critical importance.
In this study, data from two recent large-scale sampling investigations (2012 and 2014-2015) in Xinjiang Province, China were further analyzed. To investigate the distribution characteristics of BTV, a global and local spatial autocorrelation analysis was performed. Risk areas for BTV occurrence in Xinjiang were detected using the presence-only maximum entropy (MaxEnt) model. This work presents the risk zones of BT in Xinjiang Province, which may subsequently provide useful information for the development of effective strategies for the prevention and control of BT outbreaks.

Study area
The Xinjiang Uygur Autonomous Region is located along the northwestern border of China ( Fig. 1), extending from latitude 34 22′ to 49 10′N and longitude 73 40′ to 96 23′E. The region has a temperate continental climate. According to the data released by the National Bureau of Statistics of the People's Republic of China (http://www.stats.gov.cn/), Xinjiang is the second largest province in terms of sheep and goat production in China.

Environmental data collection
To characterize the environmental requirements for BTV to be present, 21 environmental factors that may influence BT occurrence were initially selected (García-Bocanegra et al., 2018;Remya, Ramachandran & Jayakumar, 2015). Environmental factors included 19 bioclimatic variables (Bio 1-Bio 19) and gridded sheep and goat densities (SD and GD). The 19 bioclimatic variables obtained from WorldClim (http://www.worldclim.org/), represent annual trends, seasonality and limiting or extreme environmental factors. Two animal distribution variables were obtained from Livestock Geo-Wiki (https://livestock.geo-wiki.org/home-2/), and represent sheep and goat densities. All the environmental parameters were converted to ASCII raster grids and preprocessed to a spatial resolution of 30 arc-seconds (ca. one km 2 at ground level).

Spatial autocorrelation analysis
The distribution of a phenomenon is presented as either a clustered, dispersed, or random pattern within a given space. Spatial autocorrelation analysis was used to investigate geographic patterns of BTV distribution in the Xinjiang Province. BTV positive rate was taken as the attribute value. Both local and global spatial autocorrelations were used to analyze the datasets.

Global spatial autocorrelation analysis
A global spatial autocorrelation was applied to analyze the distribution of BTV where all counties were seen as a whole. Global Moran's I measured the spatial autocorrelation of county locations and the BTV positive rate. Global Moran's I ranges from -1 to 1, which correspond to highly dispersed and highly clustered distributions, respectively. These parameters were calculated as follows (Cliff & Ord, 1973): where: X i = the BTV positivity rate in the ith county; X = the mean of the BTV positivity rate in all counties within Xinjiang Province; X j = the BTV positivity rate in the jth county; W ij = a weight parameter for the pair of counties i and j that represents proximity; n = the number of counties in Xinjiang Province.

Local spatial autocorrelation analysis
Local spatial autocorrelation was applied to explore the BTV distribution mode within a particular county. The local Getis-Ord G Ã i statistic (Hinman, Blackburn & Curtis, 2006) and its Z-value were calculated to test for statistical significance of BTV local autocorrelation values. If G Ã i > 0 and Z > 1.96, the county would be considered as a hotspot area, indicating that BTV distribution within this province were spatially clustered with a significance level of 99% (p < 0.01). Getis-Ord G Ã i was calculated as follows (Getis & Ord, 2010): where: X i = the BTV positivity rate in the jth county; W ij = a weight parameter for the pair of counties i and j that represents proximity; X = the mean of the BTV positivity rate in all counties within Xinjiang; n = the number of counties in Xinjiang; S = the standard deviation.

Maximum entropy modeling
Risk areas for BTV infection in Xinjiang were detected using the presence-only MaxEnt ecological niche model (Phillips, Anderson & Schapire, 2006). BTV positive livestock samples were used as proxies for the vector-borne BTV to model the virus presence. Correlation among environmental factors was assessed using Spearman's rank correlation coefficient to avoid variable multicollinearity that can result in model over-fitting (Graham, 2003). A cross-correlation value less than 0.75 was used as the cut-off threshold to exclude highly correlated variables (Zhang et al., 2012). Variables considered co-linear were excluded and nine variables were selected as evaluator variables. The MaxEnt model used BTV positive locations in 2012, 2014 and 2015 as presence data and 10,000 randomly chosen background points as "Pseudo-Absence" data. To reduce environmental bias resulting from sampling bias introduced from spatially clustered occurrences, all the locality data of 164 BTV positive samples were rarefied at one km 2 . A total of 112 occurrence locality data remained for model development after rarefying the database using the Spatially Rarefy Occurrence Data for SDMs tool (Brown, Bennett & French, 2017). A total of 75% of the BTV occurrence locations were used as the training set for model calibration, and 25% were used as the testing set for model evaluation. A regularization value of one was used to avoid over fitting of the test data. Area under the curve (AUC) of the receiver operating characteristic plot was calculated to evaluate the produced model (Phillips, Anderson & Schapire, 2006). AUC ranges from zero to one, with one indicating perfect discrimination (Fielding & Bell, 1997). The importance of each variable was assessed using the Jackknife test and percent contribution. The potential species distribution map had a series of values from zero to one which indicated low potential to high potential. These values were regrouped into four classes of potential habitats with high potential (>0.421), moderate potential (0.296-0.421) and low potential (<0.296) based on 10th percentile presence threshold (0.296) (Hannah & Edwin, 2012)

Ethics statement
Ethics Committee approval was obtained from the Laboratory Animal Ethics Committee of Northeast Agricultural University to the commencement of the study.

C-ELISA for BTV
During 2012, a total of 1,441 sheep and goat blood samples were collected to assess the prevalence of BT. A total of 19 samples were determined to be positive for BTV (1.32%). The rate of BTV positive samples was the highest in Yuli County, at 3.33%. Qapqal County was found to be free of BTV in 2012. The detailed data are presented in Table 1.
In 2014 and 2015, A total of 2,135 sheep and goat blood samples were collected to determine the prevalence of BT. There were 145 samples found to be BTV positive (6.79%). The rate of BTV positive animals in Hejing County was the highest at 15.33%. Turpan Prefecture was observed to be free of BTV. The detail information is shown in Table 2.

Spatial autocorrelation analysis
Global spatial autocorrelation analysis The results of the global spatial autocorrelation analysis are presented in Table 3. The 2012 distribution of BTV in Xinjiang followed a random pattern, but was significantly clustered in 2014-2015.

Local spatial autocorrelation analysis
The results of the 2012 BTV hotspot analysis in Xinjiang Province are presented in Fig. 1 and Table S1. As is demonstrated in the figure, Balikun County (Z = 3.1256,  Figure 2 showed the result of the hotspot analysis of BTV distribution in Xinjiang from 2014 to 2015. As shown in the figure, during 2014 and 2015, Balikun County (Z = 4.0818, p = 0.0018), Yiwu County (Z = 5.6000, p < 0.001) and Hami City (Z = 4.9151, p < 0.001) were considered to be BTV infection hotspot areas. Table S1 presented G i Ã Z scores and p-values of the hotspot areas.

MaxEnt modeling
Area under the curve score for the training data was 0.880, indicating that the approach fit the training data fairly well. AUC score for the test data was 0.876 (SD = 0.029), also indicating that the model performed well. As is shown in Table 4, five variables, including sheep distribution (25.6% contribution), precipitation seasonality (22.1% contribution), mean diurnal range (16.2% contribution) and isothermality (16.9% contribution), provided over 80% of model contribution. Sheep distribution, precipitation seasonality and mean diurnal range were identified as the most important predictors for BTV occurrence in Xinjiang. Results of Jackknife test and response curves were shown in Figs. S1 and S2. Figure 3 was a representation of the MaxEnt model for BTV. Warmer colors showed areas with better predicted conditions for BTV occurrence, and these areas were identified as "high-risk" areas.

DISCUSSION
The first case of BT in China was identified in Yunnan Province in 1979, and BTV was isolated after the outbreak. Soon thereafter, cases of BT in Hubei, Anhui, Sichuan, Gansu and Shanxi Provinces were reported. Concurrently, BTV seropositive animals were  found in 29 provinces throughout China, including Guangdong, Guangxi, Jiangsu, Xinjiang, and others, indicating a rapid spread of BT throughout the country. A previously published study covering 27 provinces throughout China between 1987 and 1989 indicated a nationwide BTV seroprevalence rate in sheep and goats of 4.73%. In a similar study, the highest seroprevalence was observed in Guangxi Province in 2001, with 31.7% of sheep and goats testing positive for BTV infection. In Inner Mongolia, the BTV seropositivity rates were 11.75% in 2014 and 11.27% in 2015. It has been demonstrated that BT has become widely established in China for several decades. In Xinjiang Province, China, three large-scale epidemiological surveys were conducted since the initial cases of BTV were confirmed. The first epidemiological survey analyzed a total of 160,671 blood samples collected from sheep, goats, cattle, yaks and deer from 1988 to 1989 to determine rates of seroconversion to BTV. The result confirmed for the first time that BT was widespread in Xinjiang Province, with goats exhibiting the highest rate of seroconversion to BTV. The second epidemiological survey was conducted in 2012. In total, 1,441 blood samples from sheep and goats, as well as 701 from cattle were collected. The observed rate of BTV positive samples was 1.32% (19 in 1,441) in sheep and goats and 0% (zero in 701) in cattle. However, it has been demonstrated that the rate of BTV positive sera collected in Southern Xinjiang is higher than was observed in Northern Xinjiang. The third and most recent epidemiological survey covered the period from 2014 to 2015. During these two years, 2,135 blood samples collected from sheep and goats were tested, and the rate of BTV positive sera was 6.79% (145 in 2,135). The investigation indicated that the epidemic status of BT in Xinjiang remains a significant concern. Spatial epidemiology plays an important role in the study of infectious diseases in the field of public health. One such method, spatial autocorrelation, has been used widely in epidemic studies (Al-Ahmadi & Al-Zahrani, 2013;Ma et al., 2017;Ratovonirina et al., 2017). According to our research, in 2012, the distribution of BTV was observed to follow a random pattern, considering Xinjiang Province as a whole. However, local spatial autocorrelation analysis showed that Balikun County (p = 0.0018), Yiwu County (p < 0.001) and Hami City (p < 0.001) were hotspots of BTV infection. It has been proposed by others that case clusters which occur randomly also have an effect on the spread of an infectious disease (Jeefoo, Tripathi & Souris, 2011). During 2014 and 2015, both global and local spatial autocorrelation analyses demonstrated a significantly clustered distribution of BTV in Xinjiang, and that Balikun County (p < 0.001), Yiwu County (p < 0.001) and Hami City (p < 0.001) were also hotspot areas.
In 2011, 48 Culicoides species (Diptera: Ceratopogonidae) were recorded in Xinjiang, China (Tian, Liu & Jia, 2011). Adult female haematophagous midges of the Culicoides spp. are the only known vectors and the only known mode of transmission through which BTV can spread between susceptible ruminant hosts (Toit, 1944). As such, the epidemic distribution of BT is closely related to the activity of the midge vector (Tabachnick, 2004), which is generally distributed between 40 N and 35 S latitude (Gibbs & Greiner, 1994). However, the hotspot areas for BTV infection observed during the study period lie outside of this geographical range, with the northernmost hotspot being located at 45 N latitude. Additionally, the MaxEnt model identified BTV high-risk areas north to 47 N (Fig. 3) in Xinjiang Province. It is possible that this observation is the result of the expansion of the habitat range of the midge vector due to climate change. It has been reported that the distribution of Culicoides are predicted to move northwards up to 53 N latitude with changing climatic and environmental conditions (Zuliani et al., 2015).
Maximum entropy calculates the relationship between the presence data and some environmental predictors which were known to be related to the disease (García-Bocanegra et al., 2018). This study successfully built a presence-only MaxEnt model relying on climatic and environmental data. Sheep distribution, precipitation seasonality and mean diurnal range were identified as the most important predictors for BTV occurrence. Areas with sheep density higher than 500 heads/km 2 were found to be areas of high risk for BTV (Fig. S2). It has been reported that sheep were most susceptible to BTV infection among all the small ruminants (Coetzee et al., 2012). Moreover, a decrease in the variation coefficient of precipitation seasonality resulted in lower risk of BTV presence. Seasonality of BTV occurrence has been investigated, and the infection regularity in different seasons was observed (Ward, 1996). In this study, the relationship between precipitation seasonality and BTV occurrence may be influenced by the extremely arid climate throughout a whole year in Xinjiang Province. The response curve of Bio 2 showed that areas with mean diurnal range of 14 C were found to be areas of high risk for BTV. The effect of temperature and precipitation on BTV infection has been investigated (Brand & Keeling, 2017), as well as the effects of seasonal and meteorological parameters on the Culicoides existence (Ander, Meiswinkel & Chirico, 2012;Racloz et al., 2008;Sanders et al., 2011). However, influence of annual mean temperature and precipitation on BTV occurrence was not observed in our research. As we know, both biological factors (vegetation, human and animal activities etc.) and natural environmental factors (light, temperature, atmospheric gases etc.) have major impacts on the spread of Culicoides spp, which is the major transmission mode (Blanda et al., 2018). Thus, these factors can have a significant influence on the geographical distribution of BTV. For example, an association between BT disease diffusion and some landscape features has been reported (Guis et al., 2007).
As is stated above, data of BTV positive blood samples were used as the input data for the MaxEnt model. Using livestock samples to represent vector-borne BTV may cause the result of niche models and cluster maps to not agree. Future studies should be conducted considering more predictors, including vector, vegetation and other environmental and climatic variables unavailable in this study. Furthermore, factors affecting the distribution of haematophagous midges require much consideration.
Although we believe our research is comparably reliable, limitations still exist. The apparent prevalence may differ from the true prevalence caused by serological tests. A low number of positives with a low number of tests will over inflate the estimates of BTV prevalence, although Empirical Bayes smoothing was performed. The sensitivity and specificity of the test should be performed in the future study. Furthermore, the result of MaxEnt and the hotspot maps do not agree perfectly. Maybe it is because BTV positive cases were used as a proxy for the vector when developing niche models. And the sampling biases, including animal movement and low sample numbers, may also lead to different results.

CONCLUSIONS
The global spatial autocorrelation data on the distribution of BTV in Xinjiang in 2012 exhibited a random pattern, which became markedly clustered in 2014-2015. The hotspot areas for BTV infection included Balikun County, Yiwu County and Hami City in 2012. These three regions were also hotspot areas during 2014 and 2015. A BTV suitability map was generated to show the high-risk areas for BTV occurrence in Xinjiang. This study can serve as a tool for the development of preventative countermeasures for future BT outbreaks.