Culex pipiens distribution in Tunisia: Identification of suitable areas through Random Forest and MaxEnt approaches

Abstract

Background

Tunisia has experienced several West Nile virus (WNV) outbreaks since 1997. Yet, there is limited information on the spatial distribution of the main WNV mosquito vector Culex pipiens suitability at the national level.

Objectives

In the present study, our aim was to predict and evaluate the potential and current distribution of Cx. pipiens in Tunisia.

Methods

To this end, two species distribution models were used, i.e. MaxEnt and Random Forest. Occurrence records for Cx. pipiens were obtained from adult and larvae sampled in Tunisia from 2014 to 2017. Climatic and human factors were used as predictors to model the Cx. pipiens geographical distribution. Mean decrease accuracy and mean decrease Gini indices were calculated to evaluate the importance of the impact of different environmental and human variables on the probability distribution of Cx. pipiens.

Results

Suitable habitats were mainly distributed next to oases, in the north and eastern part of the country. The most important predictor was the population density in both models. The study found out that the governorates of Monastir, Nabeul, Manouba, Ariana, Bizerte, Gabes, Medenine and Kairouan are at highest epidemic risk.

Conclusions

The potential distribution of Cx. pipiens coincides geographically with the observed distribution of the disease in humans in Tunisia. Our study has the potential for driving control effort in the fight against West Nile vector in Tunisia.

INTRODUCTION

West Nile virus (WNV; family Flaviviridae; genus Flavivirus) is an important zoonosis, worldwide distributed, responsible for mild fever to severe, lethal neuroinvasive disease in humans, birds, horses and other wildlife species (Ulbert, 2019). Since its discovery, WNV has caused multiple human and animal disease outbreaks. Infections are associated with economic costs, mainly due to the price of treatment of infected patients, control programmes, and loss of animals and animal products (Barber et al., 2010; Habarugira et al., 2020). WNV is primarily maintained in a bird host–Culex mosquito vector transmission cycle (Colpitts et al., 2012). Occasionally, WNV is transmitted to other vertebrate hosts, including horses and humans, and can cause disease. In humans, most WNV infections produce mild symptoms, but occasionally infections may be severe, causing neurological impairment and even death (Murray et al., 2010). Species belonging to the Cx. pipiens complex are the most widely distributed mosquito species in the world and two biotypes of Cx.pipiens (Linnaeus, 1758) are present in Tunisia, Cx. pipiens pipiens and Cx. pipiens molestus (Krida et al., 2015). Three major West Nile epidemics have affected humans in Tunisia in 2003, 2007 and 2012 (Hammami et al., 2017).

The spatial distribution of vector-borne zoonotic pathogens heavily depends on environmental features and of course, on the presence of both host and vector required for their maintenance in natural foci (Sage et al., 2017). The presence of endemic species like Cx. pipiens, is of public health relevance. Though, albeit the transmission mode of the virus via the mosquito vector is known, the capacity of the public health authorities and scientific community to forecast and fathom the risks of transmission remains still insufficient. Actually, it remains yet to be fully understood whether the mechanisms and conditions that caused the extraordinary increase in West Nile virus cases in Tunisia during 2007 and triggered the WNV outbreak in 2012 were singular or potentially forecastable. Among other variables, the spatial distribution and seasonality of vectors, as well as their interconnections with hosts, are key factors that can significantly impact the risk of dissemination (Gould and Higgs, 2009; Candeloro et al., 2020). Hence, it is crucial to understand how Cx. pipiens is distributed across the environment and what natural and anthropogenic determinants influence its spatial distribution (Manica et al., 2020).

Notwithstanding the progresses in WNV control measures in Tunisia recently (WHO, 2018) with the achievement of an agreement between the Ministries of Health and Agriculture, WNV continues to expand with the death of a person in Msaken in September 2018 and another one in Gabes in October 2018. A total number of 49 human positive cases were confirmed nationwide in several governorates (Sousse, Monastir, Mahdia, Kairouan, Sfax, Sidi Bouzid, Bizerte, Gabès, Tunis, Béjà, Zaghouan, Gafsa) (ONMNE, 2018).

We developed two species distribution models (SDMs) to predict the habitat suitability of Cx. pipiens (Linnaeus, 1758), the main WNV vector mosquitoes species in Tunisia, presenting a country-wide mapping of areas potentially occupied by this species and assessing the importance of different environmental variables in modelling their ecological niches.

Maximum entropy (MaxEnt) is a presence only modelling technique that outperforms most of the other algorithms in terms of both predictive success and accuracy of predictions (Hernandez et al., 2006; Phillips et al., 2006; Ortega-Huerta and Peterson, 2008). Moreover, MaxEnt predictions tend to be relatively conservative and do not overpredict the area where vectors may be present, a common drawback of other algorithms (Moffett et al., 2007; Larson et al., 2010).

The Random Forest (RF) methodology consists of an ensemble of classification and regression trees built on a random subset of both the available samples and the attributes recorded for each data point (Breiman, 2001). RF performs well in the prediction of species distribution with few samples in large undersampled areas (Walter et al., 2018; Arora et al., 2022).

The aim is to predict the most suitable areas for Cx. pipiens mosquitoes using RF and MaxEnt modelling techniques.

MATERIALS AND METHODS Study area

Tunisia is located on the northern coast of Africa. The north and eastern parts of the country face the Mediterranean Sea, with a typical Mediterranean landscape. Tunisia borders with Algeria to the west and Libya to the southeast. Geographically, the southern part of Tunisia is occupied by the Sahara desert; mountains are present in the northwestern part (Atlas Mountains) while the coastal and northern-eastern parts of the country are cultivated. Tunisia occupies an area of 163,610 km², 8250 of which are represented by waters. Tunisia has 38 wetland zones of international importance especially as waterfowl habitat (Ramsar) (Abid, 2013). The southern part of the country has not been considered in all the analyses as it is represented mainly by desert.

Vector presence data

A field campaign of mosquitoes and larvae collection was carried out in Tunisia from 2014 to 2017. To select the sites where to put the traps, first bibliographic research was conducted to identify the areas where the West Nile virus was previously detected or where suitable areas for WNV circulation were described based on human cases reported and climatic and environmental conditions (Ben Hassine et al., 2017). Mosquito traps were placed in urban (houses and gardens) and rural (riding stables and livestock farms) environments. CDC miniature light traps (BioQuip Products, Inc., Rancho Dominguez, USA) baited with CO₂ (1 kg of dry ice) were used for the adult mosquito collections. Traps were located in outdoor and indoor habitats at 6 pm and grouped the next morning. Mosquito larvae were collected from a variety of breeding sites by the classical dipping method using plastic dippers. Adult specimens were fixed and 4th instar larval specimens were mounted on microscope slides and identified morphologically with local and regional identification keys (Brunhes et al., 2000).

In our study, we extracted occurrences of Cx. pipiens from a nationwide survey reporting Culex species presence/absence in Tunisia, carried out from 2014 to 2017. The dataset obtained from the field campaign carried out (in 234 sites), was integrated with a literature dataset of mosquito occurrence in 22 sites across Tunisia (Beji et al., 2017). A total number of 256 places have been investigated focusing the modelling on Cx. pipiens environmental suitability (Krida et al., 2015; Gangoso et al., 2020) across Tunisia (Figure 1).

FIGURE 1. Geographical location of the sampling sites investigated in this study (Cx. pipiens occurrence in red circles, absence in green and the locations from Beji et al., 2017 in orange triangles)

Climatic, environmental and human predictors

Several factors known to be relevant for the presence of Cx. pipiens were considered in the analyses: they refer to temperature, rainfall, vegetation, topography, presence of water, human density (Burkett-Cadena et al., 2013; Foster and Walker; 2019; Paz and Semenza, 2013; Anyamba et al., 2012; Conte et al., 2015; Conley et al., 2014).

We considered historical climatic conditions of the country to derive the main abiotic factors, i.e. the physical features of different Tunisian territories. The idea is to catch the known main differences across the country, in terms of temperatures, precipitations, vegetation abundance and measured as climate metrics (annual trends, seasonality and extreme or limiting environmental factors). Tunisia's climate varies due to the country's diverse geography, which can be divided into three regions: the northern mountainous region has a Mediterranean climate with mild, rainy winters and hot, dry summers; the south has a hot, dry, and semi-arid climate as it enters the Sahara Desert; the eastern coastal border has an arid steppe climate (Ghoneim et al., 2017; Baban et al., 1999).

Considering that at local scale, climatic and environmental changes can occur compared to the historical conditions (e.g. urbanisation areas, desertification areas), to take into account also the current environmental conditions and to characterise the recent vegetation dynamics, other vegetation indices and land surface temperature were also included for the years 2016 and 2017. These determinants are relevant to provide suitable climate conditions for the proliferation of mosquitoes and as resting sites for birds/hosts species: both aspects have influence on the transmission cycles (Tran et al., 2014; Paz and Semenza, 2013).

Two water-related predictors were also considered, reporting the percentage of stagnant, clear or polluted waters due to natural wetlands or due to irrigated fields. Cx. pipiens mosquitoes lay their eggs in water, and larval stages are aquatic, so aquatic habitats are a prerequisite for mosquito populations, and the maintenance of suitable humid conditions for a while, is a condition affecting the abundance of mosquitoes (Rosà et al., 2014).

Cx. pipiens is able to breed and feed in rural ecosystems as well as in human-altered landscapes. Moreover, it was found in human densely populated area, so we also included human population density as a predictor (Paz and Semenza, 2013; Gangoso et al., 2020).

We collected in total thirty-one layers, whose description and source are reported in Table S1. Here follows the description of data download and manipulation to get the data ready for statistical analyses.

The historical climatic conditions are referred to the years 1970–2000: we downloaded 19 ‘bioclimatic’ variables (Fick and Hijmans, 2017) from the WorldClim database (http://www.worldclim.org/) at a spatial resolution of 30 arc seconds (approximately 1 km²). The complete list and a brief description for these environmental predictors are reported in Table S1.

Normalised Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI) were extracted from MOD13Q1 MODIS-NASA product (250 m spatial resolution and a temporal resolution of 16 days) and night time land surface temperature (LSTN) and daytime land surface temperature (LSTD) were extracted from MOD11A2 NASA product (1 km spatial resolution, temporal resolution of 8 days). Data were downloaded from the Land Processes Distributed Active Archive Center (LP DAAC) service at NASA website, (http://e4ftl01.cr.usgs.gov/MOLT, accessed on the 7 October 2020). Data were first submitted to a pre-processing phase of missing data interpolation to fill in pixels with no data due to cloud cover. All the data were then aggregated and mean, maximum value and standard deviation were calculated per pixel for the 2-year period.

Two water-related predictors were derived from the agricultural map of Tunisia (National technical report, confidential data), rasterised and calculated as the percentage of pixel covered by waters. The first layer (IRRIGATED_LAND) take into account waters in irrigated fields, that is the percentage of pixel equipped for irrigation, which corresponds to human-influenced habitat. Irrigated areas are constituted of ‘private perimeters’ and collective or ‘public perimeters’: approximately 40% of the perimeters are private and irrigated from shallow or deep individual wells and shallow drillings, and 60% of the perimeters are supplied for public irrigation from deep wells and great and hillside dams (https://water.fanack.com/tunisia/water-use-tunisia/ last accessed on 7 October 2020). The second layer (HYDROZONES) reports the percentage of pixel occupied by natural surface waters: lakes, rivers, oued (streambeds that remain dry except during the rainy season), chott (shallow saline lakes), sebkha (a smooth flat often saline plain sometimes occupied after a rain by a shallow lake) and garaa (bodies of fresh water in the interior of the country either temporary and permanent).

The human distribution was extracted from the global rural urban mapping project at https://sedac.ciesin.columbia.edu/data/set/grump-v1-population-density/data-download (accessed on the 7 October 2020) that provides the population density with a spatial resolution of 5 km.

Altitude was derived from the digital elevation model (DEM) of the Global 30 arc second elevation data (GTOPO) from https://doi.org/10.5066/f7df6pqs (accessed on the 7 October 2020), with a spatial resolution of approximately 1 km².

All layers were converted into WGS84 Longitude-Latitude geographic projection (when needed), clipped on Tunisia extent and resampled in raster format at 5000 m spatial resolution.

The software ESRI^® ArcGIS 10.6 was used for geographical data management and manipulation, and for the map production. Furthermore, the southern part of the country has been masked in all the analyses as it is represented mainly by desert.

Statistical analyses and modelling

Considering the high number of variables, we calculated the pairwise correlation matrix among predictors, and selected those variables with a Pearson's correlation coefficient less than 0.8 to reduce multi-collinearity problems. This resulted in 13 layers (out of the initial 31 layers) including two vegetation layers, one topographic layer, five bioclimatic, two hydrological layers, two temperature layers and one population layer (Table 1).

TABLE 1 List of climatic, environmental and human predictors selected for Cx. pipiens habitat suitability models

Abbreviation	Climate and environmental parameters	Min value	Max value	Unit of measure
LSTD	Land surface temperature daytime – average of the years 2016–2017	18.55	38.28	°C
LSTN	Land surface temperature night-time – average of the years 2016–2017	9.82	21.69	°C
NDVI MEAN	Average annual NDVI (2016–2017)	−0.07	0.8	Index value
NDVI STD	Standard deviation of NDVI (2016–2017)	0.00	0.28	Index value
DEM	Elevation	−29	1407	Metres
Bio 15	Precipitation seasonality	26.34	104.19	Percentage
Bio 3	Isothermality	34.03	46.23	Percentage
Bio 7	Temperature annual range	22.17	38.18	°C
Bio 8	Mean temperature of the wettest quarter	5.35	22.17	°C
Bio 9	Mean temperature of the driest quarter	22.01	31.70	°C
HYDROZONES	Percentage covered by water	0	100	Percentage
IRRIGATED_LAND	Percentage covered by water	0	100	Percentage
POPULATION	Density of population	0	25,742	Number

We determined the probability of Cx. pipiens occurrence by estimating the relationship between species records at sites and the environmental and/or spatial characteristics of those sites.

Two different machine learning models, MaxEnt and RF, were implemented to predict the occurrence of Cx. pipiens.

MaxEnt, a machine learning algorithm, is designed for use with presence-only data (Elith and Graham, 2009) and has been widely used for modelling species distributions and regularly outperforms similar models (Elith et al., 2011). It generates relative scores estimated from presence data and environmental predictors and minimises the relative entropy between them (Elith et al., 2011). The importance of variable was assessed through the per cent contribution metric (Phillips et al., 2006).

The RF model is based on classification trees and additional bootstrapping (Breiman, 2001; Cutler et al., 2007; Evans et al., 2011). To estimate the optimal number of predictors which were randomly selected as candidate variables at each bootstrapping step, the model was tuned using a number of 1000 trees. The importance of each variable was assessed as the total decrease in node impurities from splitting on the variable, averaged over all trees. For classification, the node impurity is measured by the Gini index. The Gini index represents the probability of a particular variable being wrongly classified when it is randomly chosen, so its values range from 0 to 1. How much the Gini index decreases for a feature at each split define the variable importance during RF training. The more it decreases, the greater the features importance based on this index. The two models differ for type of input data: MaxEnt is a presence-only model that uses absence data for the evaluation of the model performance, while RF requires both presence and absence data to predict the species distribution.

Because only few absence data were available for the analyses, additional pseudo-absence points were randomly generated inside the extent of the region of Cx. pipiens occurrence (drawing a minimum convex polygon) and added to true negative points in order to obtain comparable numbers in both presence and absence records.

The two models were run using the same dataset made of presence data (n = 200 points), absence data (n = 56 true negatives points) and pseudo-absence points (n = 200 points).

To compare the model performances, the area under the Receiving Operator Curve (AUC) was evaluated for MaxEnt and RF. The value of AUC ranges from 0 to 1. An AUC value of 0.50 indicates that the model did not perform better than random, whereas a value of 1.0 indicates perfect discrimination (Swets, 1988). However, the AUC evaluated for the final models (fitted using all available data and then used for prediction), may be overestimated, because the training data are the same used to evaluate the model (Hastie et al., 2004). To measure AUC in a less optimistic way, so to have a more precise idea of the real predictive capability, we applied a k-fold approach (k = 5) for both the models. A probability threshold of 0.5 was adopted for positive/negative classification and for creating the confusion matrix so to calculate classification metrics for the final models.

All the analysis were run in R statistical environment version 3.5.3 (R Core Team, 2020) using the main following packages: dismo (Hijmans et al., 2020), randomForest (Liaw and Wiener, 2002), rasterVis (Lamigueiro and Hijmans, 2021), rgdal (Bivand et al., 2021) and ENMeval (Kass et al., 2021). Maps were elaborated using ESRI^® ArcGIS 10.6.

RESULTS Cx. pipiens sampling campaign

The map in Figure 1 reports the sampling locations across Tunisia. During the study period (2014–2017), specimens of Cx. pipiens were collected; the number of presence/absence of Cx. pipiens species by bioclimatic zone is shown in Table 2. Cx. pipiens was distributed in the humid, sub-humid, semi-arid, arid bioclimatic zones and became numerous in the higher semi-arid and lower arid bioclimatic zones.

TABLE 2 Description of Cx. pipiens presence/absence sites across Tunisia

Bioclimatic zone	Governorate/delegations	Number of presence/absence per delegation	Month of trapping
Higher arid temperate	Sidi Bouzid /Jelma,Ouled Haffouz, Regueb Sidi Bouzid/Meknassi,Sidi Ali Ben Aoun Sidi Bouzid/Sidi Bouzid Est, Sidi Bouzid Ouest Kairouan/Sbikha	3P 4A 5P+1A 5P+6A	September 2015 September 2015 September 2015 November 2015
Higher saharan	Tozeur/Tozeur Tozeur/Degueche Kebili/Kebili Nord Kebili/Souk El Ahed Kebili/Douz	3P 3P+2A 8P+1A 6P 3P+4A	December 2015 December 2015 December 2015 December 2015 December 2015
Higher semi-arid	Tunis/Bab Bhar, El Ouardia, La Marsa, Sidi El Bachir, Sijoumi Nabeul/Beni khiar Nabeul/Bou Argoub Nabeul/Grombalia Nabeul/Korba Nabeul/Menzel Bouzelfa Nabeul/Soliman, Takelsa Ariana/Soukra, Kalaat El Andalous, Raoued, Sebket Ariana Siliana/Bargou, El krib, Makthar, Siliana nord Le Kef/Kef Est, Kef Ouest	6P 1P 2P 3P 4P+1A 1P 1P+1A 4P 6P+9A 3P+1A	December 2014 November 2014 November 2014 November 2014 November 2014 November 2014 November 2014 November 2014 September 2016 November 2015
Lower semi-arid	Sousse/Msaken, Bouficha, Kalaa Sghira, Sidi Bou Ali Monastir/Moknine, Monastir Monastir/Bembla, Jammel, Ksar Hellal Monastir/Sahline, Sebkhat El Moknine, Zeramdine Mahdia/Souassi Mahdia/Sidi Alouane Mahdia/Ksour Essef, Mahdia, El Jem	10P 12P+2A 6P 3P 4P 1P 4P	October 2015 October 2015 October 2015 December 2014 September 2015 October 2015 November 2015
Lower humid	Beja/Nefza Beja/Testour Beja/Teboursouk Beja/Beja sud Bizerte/Sejnene Bizerte/Mateur Bizerte/Utique, Tinja Bizerte/Menzel Bourguiba, Bizerte sud Jendouba/Tabarka Jendouba/Jendouba Nord, Jendouba Sud, Bousalem Jendouba/Fernana, Ghardimaoui	2P 3P 3P 4P+1A 2P 9P 2P 4P 6P+1A 6P+2A 2P	September 2015 December 2014 November 2014 November 2014 December 2015, January 2017 December 2015, January 2017 November 2014 November 2014, January 2017 September 2015 November 2014 May 2015
Middle semi-arid	Zaghouan/Zaghouan, Fahs, Zriba, Bir Mchergua	4P+6A	May 2016
Higher arid	Kasserine/Kasserine Nord, Ezzouhour Kasserine/Thala, Foussana	2P+2A 2P	May 2016 December 2015
Lower arid	Sfax/Sfax Ouest, Sfax Est Sfax/Jebeniana, Sfax Medina, Agareb Medenine/Medenine Sud,Medenine Nord Gabes/Gabes Medina, Gabes Ouest, Gabes sud, Mareth Gafsa/Ksar	3P 3P+2A 8P+1A 28P+8A 1A	December 2015 December 2015 December 2016 December 2016 December 2016

P, presence; A, absence.

Modelling

The environmental suitability for Cx. pipiens was investigated in Tunisia using two statistical models, i.e., RF and MaxEnt. Through these modelling techniques, occurrence data collected in 256 locations were linked to environmental factors. The k-fold approach run for RF and MaxEnt produced average AUCs values of 0.91 (standard deviation = 0.037) and 0.86 (standard deviation = 0.059), respectively, indicating better performances of the RF model. The Figure 2 shows the ROC curve for the final models and the AUC for the final model (lines), along with the AUC values obtained during k-fold (boxplot of the five measures). The AUC for final models derived using all the data available resulted in 0.99 (RF) and 0.91 for (MaxEnt). The misclassifications for the final models are reported in Tables 3 and 4 assuming a threshold of 0.5 for MaxEnt and RF respectively. RF is the most specific and sensitive method (specificity = 0.97, sensitivity = 0.95, Table 4).

FIGURE 2. AUC for the final models is represented by lines, AUC from k-fold approach is summarised by boxplot

TABLE 3 The confusion matrix of the MaxEnt prediction

		True status
		Negative	Positive	Total
Predicted values	Negative	165	26	191
	Positive	33	157	190
Total		198	183	381

TABLE 4 Confusion matrix of the RF prediction

		True status
		Negative	Positive	Total
Predicted values	Negative	192	10	202
	Positive	6	173	179
Total		198	183	381

The probability distribution values found across the territory are reported through violin plots (Figure 3) where the shape represents the density estimate of the probability: the more data points in a specific range, the larger the violin is for that range.

FIGURE 3. Violin plots depicting the distribution of probabilities of Cx. pipiens occurrence for MaxEnt and RF

The thirteen predictors and their contribution for MaxEnt and RF models are reported in Tables 5 and 6. The fitted models were then used to predict the probability of Cx. pipiens occurrence for the whole territory of Tunisia, except the desert area in the south of the country (Figures 4 and 5). Coastal areas of the Sahel in the northeast, the governorates located in the north of the country, oases of the south, and along the eastern Tunisian coast from the Gulf of Hammamet along the Gulf of Gabes are the areas with the highest probability of occurrence. The central governorates represent the regions with the lowest probability of presence of Cx. pipiens. The areas with lower probability are represented by the southwest of Kasserine, the southwest of Kef and the south of Siliana. Most areas shown in the map with high risk are also indicated as areas where WNV outbreaks have been recorded (Wasfi et al., 2016; Beji et al., 2017; Hammami et al., 2017); human WND sporadic cases were recorded in the northwestern (Jendouba governorate), in the southern (Tataouine and Kebili governorates) areas of the country and in the eastern coastal districts of Tunisia (Governorates of Sfax, Mahdia and Monastir) (Ben Hassine et al., 2017; Conte et al., 2015; Hammami et al., 2017).

TABLE 5 Predictors used in the study (13 out of 31 after the correlation) and the two metrics used for importance assessment

Predictor	Per cent contribution MaxEnt	Gini-based importance
POPULATION	26.17	24.89
Bio 7	10.91	20.35
Bio 15	25.97	20.02
DEM	7.71	18.30
NDVI STD	11.81	17.74
NDVI MEAN	0.23	12.78
Bio 9	4.20	11.30
LSTD	1.61	10.36
LSTN	0.42	9.72
Bio 3	1.06	9.52
IRRIGATED_LAND	8.42	8.94
Bio 8	0.87	8.02
HYDROZONES	0.61	3.10

TABLE 6 Variation in the climatic parameters

Climatic parameter	In higher probability areas	In lower probability areas
Average temperature	22.92°C (Sahel area) to 35.4°C (oases of the south)	29.13°C (Kasserine) to 36.26°C (Siliana)
Precipitation seasonality (Bio 15)	54.6% (oases of the south) to 64.23% (Sahel)	29.13% (Kasserine) to 36.26 % (Siliana)
The vegetation layer (NDVI STD)	0.022 (oases of the south) to 0.24 (north of the country)	0.022 (Kasserine) to 0.11 (Siliana)
The vegetation layer (NDVI MEAN)	0.13 (oases of the south) to 0.42 (the north of the country)	0.14 and 0.22
Elevation (DEM)	23 m (north of the country) to 53 m (Sahel)	520 m (south west of Kef) to 917 m (south of Siliana)

FIGURE 4. MaxEnt (left) and RF (right) maps reporting the predicted probability of Cx. pipiens presence

FIGURE 5. Classification of the data points used (TP = true positive – red, TN = true negative – green, FP = FP = false positive – blue, FN = false negative – black) with the shaded probability as background

DISCUSSION

In a previous study, a model was developed based on the Mahalanobis distance methodology using human cases for the prediction and the identification of high-risk areas of West Nile in Tunisia (Ben Hassine et al., 2017). No vector distribution data were available at that time. In this study, we investigated the relationships between environmental, climatic and socioeconomic variables and mosquito presence collected across the country, by comparing two modelling approaches: RF and MaxEnt.

Our models showed consistent agreement between RF and MaxEnt methods in selection of habitat with the highest probability of occurrence of Cx. pipiens. The average predicted environmental suitability was lower for MaxEnt than for RF. It is visible in the predominance of red colour in the probability maps.

In terms of model performance, RF shows the best discrimination skills. Also in other studies, this technique was consistently reported to outperform other traditional modelling techniques (Cutler et al., 2007; Peters et al., 2007). The AUC for the two models showed good predictive power with AUC ≥0.8. RF sensitivity and specificity are excellent with values superior to 0.9. However, we need to consider that the training data are also used to evaluate the model, meaning that the accuracy measures will be overestimated (Hastie et al., 2004). To account for this, we also evaluated the AUC within a k-fold approach, resulting in an average value greater than 0.85 for both models and confirming the prediction power of these methods.

Figure S1 shows the geographical distribution of human cases of 2011–2012 epidemic in Tunisia from Ben Hassine et al. (2017) overlapped to the predicted probability maps for Cx. pipiens occurrence.

The geographical distribution of mosquitoes is influenced by several factors related to the modification of habitats and the active dispersion of mosquitoes (the displacement made by the mosquito in search of suitable breeding sites could reach 5 km) (Verdonschot and Besse-Lototskaya, 2014).

There are subtle associations between mosquito population dynamics and highly heterogeneous environmental drivers (due to a lack of clear shift between dry/wet rainfall regimes). Big cities like Tunis, Monastir, Sousse also exhibit a degree of anthropogenic command on water availability (diffuse water impoundments, irrigation) and human density. Most of the variables highlighted by the models as important are in agreement with field experience, existing biological knowledge, and known habitat preference of Cx. pipiens in Tunisia. For the variable importance, population density was the most important predictor, recorded as the most influential variable in both models. Cx. pipiens is known to breed in a wide variety of habitats and to be associated with areas with human activity. Culex mosquitoes are known to breed in areas with poor sanitary and housing conditions; thus human factors may play a more important role than precipitation.

Environmental changes and urbanisation could partially explain the concentration around big cities (Chekir and Ben Salem, 2021; Wang, 2020). The north of Tunisia is characterised by natural wetlands, big cities like Sousse, Kairouan or Sfax include bodies of water flooded by human activity, such as peri-domestic containers, agricultural channels and wells. The Sahel coast is characterised by the presence of advanced degradation due to pollution from several factories. The localised eutrophication and uncontrolled dumping of municipal wastes in several regions create favourable conditions for mosquitoes activity. Many industries discharge their effluents on the open land and into surrounding water bodies creating a favourable microclimate for mosquitoes activity. Anthropogenic sources as industrial wastewater discharges, sewage wastewater and other industrial activities can also explain the distribution map. The greatest densities of Cx. pipiens are mainly described in stagnant waters rich in organic matter.

The temperature (Bio 7) ranked second with RF and fourth for MaxEnt. West Nile virus incidence and transmission responded unimodally to average temperature (Shocket et al., 2019), since replication growth rates of pathogens and vectors are in time overburdened by hastening decreases in vector survival at high temperatures (Mordecai et al., 2013). Actually, metabolic reaction rates tend to rise significantly up to an optimal temperature, then decrease owing to protein degeneration and other processes (Dell et al., 2011). According to Ruybal et al. (2016), temperature had consistent non-linear effects. Adult female development time, daily juvenile and adult female survival decreased with increasing temperature. Increased temperature in our study corroborates these earlier findings. Hence, thermal increase or decrease can affect the distribution and abundance of mosquitoes or the pathogens they transmit, by affecting the life cycle of the mosquito. They need a certain temperature range to survive and develop. In addition, the increase in temperature can influence the reproduction rates and the incubation period (Bunyavanich et al., 2003).

The precipitation seasonality (Bio 15) was classified second with MaxEnt and third with RF. Risen precipitation could either expand or decrease mosquito abundance by creating breeding sites or washing container-breeding mosquitoes, depending upon the intensity. Our observations may support the hypothesis that increased precipitation provides more habitats for Cx. pipiens in Tunisia. Moreover, mosquito populations could either decline during droughts due to reduced breeding habitat, or increase in amplitude due to increased habitat quality or diminished predators. Besides, drought that is influenced by both precipitation and temperature could increase WNV prevalence in mosquitoes through increased contact (due to host movement to mosquito natural element), or higher vector-to-host proportions (due to aridity-induced reductions in juvenile birds) (Shaman et al., 2011; Paull et al., 2017). Wet conditions (spring and winter) facilitate an early increase of Culex mosquitoes; drier summers might then favour the collection of birds and mosquitoes around persisting water areas that hold Cx. pipiens breeding. This might increase WNV prevalence in mosquitoes and increase host–vector contact (Shaman et al., 2011; Paull et al., 2017).

NDVI is an index of photosynthetic activity, absorption and reflection of photosynthetically active radiation over a given period of time for the characterisation of the vegetation. The vegetation layer (NDVI STD) was classified third for MaxEnt and fifth for RF. These findings are in accordance with findings reported in Spain, NDVI had a positive relationship with annual and monthly presence of Cx. pipiens, the mean and standard deviation were 0.07 ± 0.06 for monthly NDVI (Roiz et al., 2015).

The (NDVI MEAN) layer ranked 13th for MaxEnt and 6th for RF. This is consistent with the work of Jian et al. (2014) as the NDVI ranged from 0.17 to 0.55 suggesting that areas with better vegetation conditions are preferred habitats with larger carrying capacity. In Egypt, Nagy et al. (2022) observed that the maximal NDVI mean in seepage water was 0.26 while the sewage water had a NDVI mean value of 0.16. The levels observed in this investigation are far above those observed in the averaged model built by Fornasiero et al. (2020), with a NDVI coefficient of 0.040. It is worth mentioning that Myer et al. (2017) indicated that the NDVI range was from 0.33 to 0.78. Our results also show that higher NDVI values may condition higher Culex presences in Tunisia.

RF picked elevation (DEM) above sea levels as the fourth variables and MaxEnt ranked the variable sixth. Indeed, Moua et al. (2021) found that the Habitat Suitability Index (HSI) increases rapidly for an altitude between 10 and 180 m but decreases to reach 0 after an altitude higher than 500 m. Our models indicated a contribution lower than 5%. Moreover the paved areas, which are linked to human presence, were less suitable for Cx. pipiens than mixed areas. In the study of Brown et al. (2008), though Cx. pipiens is an urban mosquito, its numbers were very low in highly urbanised areas. While they could be located in paved areas, the lack of feeding and resting sites owing to their ornithophilic nature makes this environment somewhat less suitable than mixed areas. Also, Gangoso et al. (2020) found an heterogeneous form of habitat suitability, with most suitable areas positioned in the southern and northeastern coastal areas of Spain, and unsuitable areas located at higher altitude and in colder zones. Actually, unsuitable areas mostly coincided with upland areas at altitudes of over 2000 m in the Pyrenees in the northeast and Sierra Nevada in the southeast, and at 1700–2000 m in the Cordillera Cantábrica in the north and Sierra de Gredos in central Spain. Suitability increasing as altitudes decreased down to 200 m with altitude having an overall negative effect. In accordance with the present results, elevation above 600 m has already been found to negatively influence the presence of the species in Spain (Alarcón-Elbal et al., 2012; Bravo-Barriga et al., 2017).

The irrigated areas (IRRIGATED_LAND) and hydrozone areas (HYDROZONES) respectively ranked fifth in the MaxEnt model and eleventh in the RF model. In the zones with higher probability of occurrence, the percentage of irrigated areas varied from 0% (north of the country, Bizerte, Jendouba) to 7.76% (Sahel). The percentage of hydrozone areas fluctuated from 0% (oasis of the south) to 6.55% (north of the country). In lower probability, the percentage of irrigated areas oscillated from 0% (south west of Kasserine) to 3.76% (south of Siliana). The percentage of hydrozone areas was 0%. In the oases of the south, groundwater use for the irrigation provides suitable water surfaces for mosquito larval development. Besides, the establishment of a wet microclimate due to the stratification and high density of plants produces suitable conditions for mosquitoes. These areas are also used as resting sites by birds during the night, increasing the probability of contact between birds and mosquitoes.

Several factors could explain the probability distribution such as human activities (urbanisation), accumulation of water in artificial containers and the pesticides use causing the mosquitoes transfer to other areas (Tabbabi et al., 2018; Tabbabi et al., 2019). The relationship between Cx. pipiens probability was positively associated with human population, precipitation, standard deviation of NDVI, mean of NDVI and irrigated areas, and negatively associated with temperature and elevation. During the mosquito trapping, historical data can be slightly different due to many reasons, but the overall spatial pattern of climatic–environmental characteristics can assumed to be kept constant in two decades across Tunisia. Changes can also occur at local scale (e.g. urbanisation areas, desertification areas) compared to the historical conditions. Making allowances also for the current environmental conditions and to characterise the recent vegetation dynamics, we considered other vegetation indices and land surface temperature related to the period years 2016 and 2017.

CONCLUSION

This study presents a description of the spatial distribution of Cx. pipiens in Tunisia. Our results revealed that the environmentally suitable areas of Cx. pipiens in Tunisia are widely distributed throughout the country with higher probability at the north of the country, in the Sahel zone and in the oases of the south. Suitable areas with lower probability correspond to the regions of Kasserine (southwest), Kef (southwest) and Siliana (southwest). Our modelling showed that RF produced higher evaluation metrics compared to MaxEnt. Human populations, temperature, precipitation seasonality as well as the standard deviation of the NDVI parameters are the major predictors of the areas characterised by the suitable conditions for Cx. pipiens distribution.

By describing the relationship between these predictors and Cx. pipiens distribution, we obtain a better understanding of how changes in land use or switching patterns of human habitat might influence disease transmission. Considering the re-emergence of mosquito borne diseases as West Nile, understanding how choices may influence disease ecology through population growth, urbanisation, is a subject of crucial importance. Unfortunately, variables on wild bird populations were not available for Tunisia and could therefore not be used in the models. Our methods do not guarantee the full removal of biases associated to sampling effort. Currently, there is no method that could.

Overall, maps that present areas of both high probability and low probability of Cx. pipiens occurrence could be used to pilot first-priority control interventions. In areas where habitat suitability evidently is low with respect to disease transmission, other factors must play essential roles which should be taken into consideration when planning health interventions. Moreover, seen in the context of limited resources for interventions, targeted WNV control could concentrate on areas suitable for both Cx. pipiens and human cases in Tunisia.

AUTHOR CONTRIBUTIONS

Jihane Amdouni: Conceptualisation, formal analysis, investigation, methodology, writing – original draft, writing – review & editing. Annamaria Conte: Conceptualisation, formal analysis, methodology, software, validation, visualisation, writing – original draft, writing – review & editing. Carla Ippoliti: Formal analysis, methodology, software, validation, visualisation, writing – original draft, writing – review & editing. Luca Candeloro: Conceptualisation, formal analysis, methodology, software, validation, visualisation, writing – original draft, writing – review & editing. Susanna Tora: Formal analysis, methodology. Soufien Sghaier: Conceptualisation, investigation, project administration. Thameur Ben Hassine: Formal analysis, methodology, validation. Emna Saida Ayari Fakhfakh: Project administration. Giovanni Savini: Project administration. Salah Hammami: Project administration, supervision, validation, writing – review & editing.

ACKNOWLEDGEMENT

This study was partially funded by IAEA (RAF50/68 regional project), Italian Ministry of Health IZS AM 01/10 RC and EU grant FP74613996 Vmerge and is catalogued by the Vmerge Steering Committee as VMerge0011.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

ETHICAL STATEMENT

The authors confirm that the ethical policies of the journal, as noted on the journal's author guidelines page, have been adhered to. The study included sampling of animals for diagnostic purposes and therefore did not require approval from the ethic committee.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1002/vms3.897.