Role of eco-climatic factors in the distribution of bluetongue in endemic areas in Tunisia.

Background: Bluetongue (BT) is an important infectious, non-contagious, OIE-listed viral disease of domestic and wild ruminants. The disease is transmitted among susceptible animals by a few species of an insect vector in the genus Culicoides. Recently, during the fall of 2020 (September and October), a Bluetongue virus-4 epizootic marked the epidemiological situation in several delegations of Tunisia with clinical cases recorded in sheep and cattle. Aim: Determine the eco-climatic variables most likely associated with delegations reporting BT cases.
Methods: A logistic regression model (LRM) was used to examine which eco-climatic variables were most likely associated with delegations reporting BT cases.
Results: Based on the LRM, our findings demonstrated that the key factors contributing significantly to BT cases' distribution among delegations in Tunisia included day land surface temperatures (DLST), night land surface temperatures (NLST) and normalized difference vegetation index (NDVI). A positive correlation between sheep distribution and rainfall amounts was demonstrated. Statistical analysis focusing on the most affected delegations during the BT epidemic (the Sahel and the Centre of Tunisia) demonstrated that the epidemic situation seems to be a consequence of the combination of the following environmental parameters: NDVI with values ranging between 0.16 and 0.2, moderate rainfall 2-4-fold above the normal (10-50 mm) and DLST values between 32°C and 34°C in September.
Conclusion: These findings suggest and develop a robust and efficient early warning surveillance program in risk areas based on eco-climatic risk factors.

Introduction

Bluetongue (BT), an emerging and re-emerging vector-borne disease (VBD) transmitted by insect vector species of the genus Culicoides, is caused by different serotypes of non-enveloped double-stranded RNA viruses in the Orbivirus genus of the Reoviridae family. Infection with Bluetongue virus (BTV) is reported in wild and domestic ruminants and can generate serious economic consequences (Talavera ; Ben Hassine ; Hassani ). In addition to the 24 well-known serotypes of BTV, other atypical novels BTV serotypes have been described and reported throughout the world (Schulz et al., 2008; Rupner ). The principal insect vector of BT in Tunisia is Culicoides imicola (C. imicola) (Hammami, 2004; Hammami ; Sghaier , 2014). Other species of Culicoides, such as C. obsoletus, C. newsteadi, C. punctatus, and C. kingi, have been shown to transmit BTV in other parts of the world. Although these species have been identified in Tunisia, their possible role in BT transmission is unclear yet (Sghaier ). Since the first detection of BT serotype 2 in Tunisia in 1999 (Ben Fredj ; Hammami, 2004), many incursions have been reported; BTV-1 in 2006 and BTV-4 in 2009 (Lorusso et al., 2012, Sghaier ). In 2016, a novel BT serotype 3 was reported (Lorusso et al., 2018). In 2018 and 2019, the resurgence of BTV-2 and BTV-4 were reported, respectively, to the World Organization of Animal Health (OIE-Wahis, 2021). Recently, in September and October 2020, the epidemiological situation was marked by the re-emergence of BTV-4. Many clinical cases were reported in sheep and cattle. The most affected areas are localized in the Sahel and the Centre of Tunisia, with over 80% reported cases. BT outbreaks were highly seasonal. They occur during the late summer and autumn, closely related to climatic and ecological factors, which influence the abundance and survival of the vectors.These factors include temperature (Wittman ), vegetation index (Carpenter ), relative humidity (Carpenter ; Kumar ), wind speeds (Ducheyne ; Carpenter ; Baylis ), higher rainfall in the preceding month of the occurrence of the disease (Gloster ; Eticha and Mekonne, 2019), elevation (Hendrickx ; Guis ), livestock distribution (Van Schaik ; Jacquot ) and soil types (Purse ). There are many dynamic climate zones in Tunisia with different agriculture practices and livestock distribution among regions.To design interventions to limit the spread of BT, it is essential to identify these factors and determine their importance in spatial disease spread. Due to the complexity involved in determining BTV and other VBD transmission, scientists often use mathematical modelling to understand the transmission dynamics of these diseases (Gubbins ; Turner ; Brand ). Spatial analytical tools have created an opportunity to evaluate environmental variables generated by remote sensing satellite sensors and analyzed by Geographic Information Systems (Tanser and Le Sueur, 2002; Vinodh Kumar ). These technologies provide useful tools for identifying and quantifying areas at high risk for viruses’ circulation and spread. This study aimed to investigate eco-climatic variables associated with BTV incursions in Tunisia, which may subsequently provide useful information for developing an early warning system and the development of effective prevention and control strategies forthe disease.

Materials and Methods

Study area and data collection

In the first phase of the study, the whole Tunisian territory was considered for the logistic regression analysis model. The administrative structure of Tunisia is composed of governorates, delegations and imadats.The delegation was chosen as the epidemiological unit and classified as an infected or non-infected delegation (Fig. 1). As of 2018, 24 governorates in Tunisia were divided into 264 delegations. Information about infected delegations was obtained from the Veterinary Research Institute of Tunisia and confirmed by the General Directorate of Veterinary Services of Tunisia (DGSV). Normalized difference vegetation index (NDVI) values were extracted in comma-separated values (CSV) format from the Global Agricultural Monitoring System (http://glam1.gsfc.nasa.gov/ last accessed on May 11, 2021). Rainfall data were extracted from the Famine Early Warning Systems Network (https://earlywarning.usgs.gov/. last access on November 10, 2021). Daily and nightly land surface temperature (LST) was extracted from MOD11A2 NASA product (1 km spatial resolution, temporal resolution 8 days) since 2010. All data were taken from the Land Processes Distributed Active Archive Centre service at NASA website (http://e4ftl01.cr.usgs.gov, last accessed on November 10, 2021). Semi-automatics Classification and Zonal statistics plugins were used in QGIS 3.10 to treat and extract data by delegation. Elevation data for Tunisia were downloaded from the Global 30-Arc-Second Elevation Data Set for the World, developed by the United States Geological Survey (http://eros.usgs.gov, last access on November 11, 2021). Sheep and cattle distributions, isused in the logistic regression model (LRM), are downloaded for each delegation from the development offices website in Tunisia (www.odco.nat.tn. last access on November 11, 2021). In a second phase of the study, environmental variables proven to be associated with the BT distribution have been elucidated and compared to the mean in two different bioclimatic floors in Tunisia, the Sahel (zone 1) and the Centre (zone 2), the most affected areas in 2020 epidemic (Fig. 1).

Fig. 1.

Study area, in the left all delegations of Tunisia are classified as infected (1) or non-infected (0) used for LRM, in the right the two zones reporting the maximum of cases (Centre and Sahel of Tunisia).

Statistical analysis

For the study’s first phase, the relationship between the eco-climatic variables and the presence or absence of BT cases in each delegation was determined using LRM, the dependent dichotomous variable. Environmental factors that may influence BT distribution were initially selected. Data were subjected to initial screening, and only parameters with significant differences were then subjected to a multivariate LRM. Analyzed environmental variables included DLST (mean), NLST (mean), NDVI (min, mean, max), elevation and precipitation. NDVI, DLST, and NLST were obtained for the year 2020 and for the month of September 2020. The level of correlation among variables was tested. The first step to identify the presence of multicollinearity was based on the correlation coefficients among the explanatory variables. Tolerance and variance inflation factor (VIF) values were used to check collinearity. This collinearity was confirmed based on condition index, eigen values, and variance proportions. Last, correlations between the remaining variables were rechecked. Minimum and maximum NDVI were highly correlated, and only mean NDVI was used in the final model. Eleven variables were used in the model. All the independent variables were continuous and converted to factors (class factor). The overall Chi-square values were examined using the output and were considered to be significantat p < 0.05. Finally, the odds ratios were determined for each variable in the model.

For the second phase of the study, significant differences between NDVI, LST, and precipitation values observed in September 2020 and the mean value calculated between 2000 and 2018 were determined using a paired comparison through the Wilcoxon Rank Sum test. The analysis is conducted separately in zone 1 and zone 2. R software, version 2.13.1 (R Development Core Team 2011), was used to perform data analyses.

Ethical approval

Ethical approval is not required for this study.

Results

Results of logistic regression for all delegations in Tunisia

Logistic regression analysis results of risk factors associated with BTV infection in sheep and cattle are presented in Tables 1 and 2. Logistic regression analysis identified NDVI and DLST registered in September 2020, and the mean values of these factors recorded in 2020 are the significant risk factors associated with infected delegations. NDVI class, which is significantly associated with infected delegations, has values between 0.15 and 0.2 in September and between 0.15 and 0.3 for the mean recorded during the same year 2020. The values of DLST class, significantly associated with infected delegations, varied between 30°C and 35°C (Table 2).

Table 1.

Univariate analysis of risk factors associated with BT infected delegations.

Risk factors	Min	Max	Mean	p value
Sheep (unite)	0	137,200	18,942	0.011 *
Cattle (unite)	0	25,596	1,753	0.278
Elevation (m)	1	935.7	199.9	0.074
NDVI mean_2020	0.079	0.669	0.280	0.023*
NDVI September 2020	0.001	0.457	0.179	0.037*
NDVI August 2020	0.003	0.597	0.187	0.826
Rainfall September (mm)	0.00	147.42	47.04	0.051.
DLST 2020 (°C)	22.24	34.36	28.60	0.000 ***
DLST Sep-2020 (°C)	27.30	44.07	33.09	0.002 **
NLST 2020 (°C)	10.52	18.29	14.45	0.352
NLST Sep-2020 (°C)	15.39	26.48	20.11	0.527

Signif. codes: 0 “***” 0.001 “**” 0.01 “*” 0.05 “.” 0.1 “ ” 1.

Table 2.

Multivariate logistic regression analysis of risk factors associated with BT infected delegations.

		p value	OR IC 2.5%	OR IC 97,5%
NDVI September	[0.1, 0.15)[0.15, 0.2)[0.2, 0.3)[0.3, 0.6]	ref0.002**0.8680.335	-1.660.310.46	-9.632.657.11
NDVI Mean 2020	[0.1, 0.15)[0.15, 0.2)[0.2, 0.3)[0.3, 0.7]	ref0.046* 0.006**0.167	-1.513.360.80	-158.93309.9280.69
Sheep distribution	[1e + 03, 1e + 04)[1e + 04, 5e + 04)[5e + 04, 1e + 05)[1e + 05, 1.5e + 05]	ref0.8280.9480.123	-0.540.190.62	-2.243.58157.32
DLST Mean 2020 (°C)	[20, 25)[25, 30)[30, 35]	ref0.4580.017*	-0.501.49	-7.2420.27
DLSTMean September2020 (°C)	[25, 30)[30, 35)[35, 40)[40, 45]	ref0.042 *0.2810.984	-1.767.371.41	-13.251.873.93
Precipitation September 2020(mm)	[0, 50)[50, 70)[70, 100)[100, 150]	ref0.7870.4430.840	-0.400.570.23	-1.963.404.39

Signif. codes: 0 “***” 0.001 “**” 0.01 “*” 0.05 “.” 0.1 “ ” 1.

The result of the logistic model shows a positive correlation between the probability of BT occurrence and the mean values of NDVI and DLST recorded in 2020. A negative correlation is observed for NDVI and DLST recorded in September 2020 (Fig. 2).

Fig. 2.

Effect of each variable used in the LRM using the “ggeffects” package in R.

NDVI

NDVI values recorded in September 2020 (0.207 and 0.185, respectively, in zone 1 and zone 2) were significantly higher than the mean NDVI value calculated between 2001 and 2018 for the same period in the two zones with p values of 0.027 and 0.008, respectively, in zone 1 and zone 2 (Fig. 3).

Fig. 3.

Zone 1 Box plots of NDVI values recorded in September 2020, mean value = 0.207 (0.160–0.293), significantly higher than the mean for the same period registered between 2001 and 2018, mean = 0.172 (0.163–0.186), p value = 0.027 (Wilcoxon Rank Sum Test for paired samples p value < 0.05). Zone 2 Box plots of NDVI values recorded in September 2020, mean value = 0.185 (0.163–0.242), significantly higher than the mean for the same period registered between 2001 and 2018, mean = 0.159 (0.153–0.173) p value = 0.008 (Wilcoxon Rank Sum Test for paired samples p value < 0.05).

Rainfall

Only the rainfall values recorded during September in zone 1 are significantly higher than the mean of the same period with 25.85 mm (3.08–54.62) and p value = 0.047 (Fig. 4). The precipitation registered in zone 2 was not significantly higher than the mean, but it can be classified as abnormal with a maximum of 20 mm four times higher than the normal (Fig. 4).

Fig. 4.

Zone 1; Box plots of precipitation values recorded in September 2020, mean value = 25.85 (3.08–54.62), significantly higher than the mean for the same period registered between 2001 and 2018, mean = 8.34 (7.06–10.88), p value = 0.00047 (Wilcoxon Rank Sum Test for paired samples p value < 0.05). Zone 2; Box plots of precipitation value recorded in September 2020, mean value = 10.32 (3.38–20.69), compared to the mean for the same period registered between 2001 and 2018, mean = 8.34 (7.06–10.88), p value = 0.07 (Wilcoxon Rank Sum Test for paired samples p value < 0.05).

Day land surface temperature

In zone 1, the mean DLST registered in September 2020 is 33.81°C (33.53–33.97). It is significantly lower than the mean 36.61°C (32.88–40.57) (p value = 0.047). In zone 2, the DLST value recorded in September 2020, 33.52°C (32.28–34.62), is significantly lower than the mean 37.1°C (33.41–41.23) (Fig. 5).

Fig. 5.

Box plots of DLST during September 2020 and compared to the mean of 2002–2018. In zone 1, the DLST recorded in September 2020, 33.81°C (33.53–33.97) is significantly lower than the mean 36.61°C (32.88–40.57) p value = 0.047. In the zone 2, the DLST value recorded in September 2020, 33.52°C (32.28–34.62) is significantly lower than the mean for the same period recorded between 2002 and 2018, mean = 37.15°C (33.41–41.23) (Wilcoxon Rank Sum Test for paired samples p < 0.05).

Discussion

Historically, BT has been circulating for a long time in tropical regions and has been considered an exotic disease in Tunisia before 1999. The disease became a worldwide threat causing high mortality rates among ruminants and a serious impact on international animal trade (Labadie ; Mayo ). Climate has already been suggested as an important driver of the BTV distribution (Samy and Peterson, 2016) closely related to arthropod lifecycle (Aguilar-Vega , 2020). Several studies have been conducted to contribute to the eradication of the disease, especially those targeting the variability of vector activity and distribution as well as abiotic factors depending on modelling (Acevedo ; Labadie ). Herein, we examined the influence of some environmental factors on BT distribution, which is reported for the first time in Tunisia. Environmental factors that may influence BT occurrence and distribution in Tunisia were initially selected based on previous studies (Talavera ; Samy and Peterson, 2016; Jacquot ; Ma ). LRM is used to obtain the odds ratio of selected environmental factors. The main advantage of LRM is to avoid confounding effects by analysing the association of all variables together (Sperandei, 2014).

Based on the LRM, our findings demonstrated that the key factors that contributed significantly to BT cases distribution among delegations in Tunisia included DLST and NDVI. The distribution of BT cases in delegation is directly associated with a favourable larval and adult habitat of C. imicola. Such observations are quite expected given that vegetation and temperature are among environmental variables related to the high risk of BT outbreak by favoring vector abundance (Diarra ; Ciss ; Cuéllar ; Gahn ). The distribution of BT in the Tunisian delegation seemed to be not significantly related to sheep and cattle distribution. But a positive correlation between the probability of BT occurrence and the sheep distribution was demonstrated. This result is in line with the study conducted by Ma . They demonstrated that areas with sheep density higher than 500 heads/km2 were found to be areas of high risk for BT distribution (Ma ). In 2007, during the epidemic in southern Spain, a significant increase in BTV-1 transmission in certain areas correlated with higher densities of sheep (Napp ). The role of cattle density, despite its possible role in the BTV distribution (Jacquot ), is yet to be clarified in the Tunisian epidemiological context. No significant association with elevation was demonstrated in this study. However, this finding is considered low reliability due to a methodological issue that considers the mean elevation in each delegation (mean = 200 m). It is assumed that C. imicola is highly abundant in lower altitudes (Verhoef ; Cuéllar ).

Environmental variables proven to be associated with BT distribution were better elucidated in two different bioclimatic floors in Tunisia, the Sahel and the Centre, the most affected areas in 2020. NDVI, DLST, and precipitation values recorded in September 2020, considered as the risk period, were downloaded and compared to the mean average values of these factors for the same period. Our findings showed that NDVI values recorded in September 2020, 0.207 (0.160–0.293) and 0.185 (0.163–0.242), respectively, in the Sahel and the Central of Tunisia, were significantly higher than the mean average; 0.172 (0.163–0.186) and 0.159 (0.153–0.173) respectively in zone 1 and zone 2. In line with our results, more sparsely vegetated habitats with values ranging between 0.27 and 0.31 and low forest cover that was exposed to full sunlight in Italy supports the occurrence of adult populations of C. imicola (Conte ). NDVI importantly inform on vector distribution that can be interpreted from the requirements of C. imicola for the humid organically enriched area as breeding sites (Kameke ). Significant high values of NDVI recorded in the Sahel of Tunisia were correlated to significantly high precipitation amounts recorded in the same period (3.08–54.62 mm) in this region. In seasons of above-average precipitation, the species C. imicola can increase its population size 200-fold (Meiswinkel, 1996). The role of precipitation was also reported in the prediction model of BT spread in Spain (Acevedo ). However, the Centre of Tunisia demonstrated no evident correlation between NDVI values and precipitation amounts. Significant high NDVI values registered in this period may be related to the irrigation practice in this region. Sidi Bouzid region, in the Central of Tunisia, is a good example of this situation. Indeed, Sidi Bouzid is known for its agricultural-based economy. The agriculture sector in this region is reinforced by the abundant water resources to irrigate large areas and fertile farmland. Irrigation practice is a critical factor for C. imicola in the risk period. Super abundant C. imicola is especially found in irrigated pastures on clayey, moisture-retentive soils where livestock is present. In countries around the Mediterranean Basin, it was confirmed that C. imicola prefers a semi-moist muddy habitat (Braverman et al., 1974). It breeds in small seepages with little free surface water created by leaks from irrigation pipes and covered by fresh grass growth (Mellor and Pitzolis, 1979). When compared to other species of Culicoides, C. imicola is found to generally prefer drier habitats.

The DLST recorded in September 2020; 33.81°C (33.53–33.97) and 33.52°C (32.28–34.62) respectively in the Sahel and the Centre of Tunisia, were significantly lower than the mean average in the two zones, 36.61°C (32.88–40.57) in zone 1 and 37.15°C (33.41–41.23) in zone 2. For BT, the temperature was shown to be an important factor in determining the magnitude of the basic reproduction number R0 (Pioz ). The distribution of C. imicola is highly dependent on-air temperature, a major limiting factor. However, the precise mechanism of this dependency is not yet known (Van Doninck ). Previously conducted laboratory studies revealed that vector competence varies with temperature and infecting BTV serotype. Vector competence for BTV-4 is lower than that for BTV-1 after incubation at 23°C–25°C (Paslaru ). The optimal temperature of 37°C for the infectious blood during artificial feeding of Culicoides was reported by several studies (Venter ; Goffredo ; Balenghien ). In 2018, it was demonstrated that, at 37°C, species of Culicoides hardly took blood with an increase of their mortality rate, so the authors opted for lower temperatures (Paslaru ).

All these results suggest a close relationship between suitable NDVI (0.16–0.3), suitable temperature (32°C–34°C), moderate rainfall above the normal, and the BT epidemic resulting from this co-variation. Under this combination, vector competence and abundance seem to reach the highest level triggering BT outbreaks in September and October 2020. A similar type of association has been reported in other BT outbreaks (Braverman ; Subramanian and Piramanayagam, 2001).

These findings suggest and develop a robust and efficient surveillance program in the risk areas based on eco-climatic risk factors. In fact, further epidemiological research should be conducted to know the status of the disease considering the specific situation of each agro­ecological zone in the country. Our choice of NDVI, precipitation, and DLST variables fit into the context of ecological niche estimation that could influence the global potential geographic distribution of BTV under present-day circumstances and give rise to more in-depth future analyses (Jones ). Further studies should focus on other environmental drivers such as soil type, wind and humidity since they were proven to be associated with BT distribution (Brand ; Kameke ).