Modelling the risk of being bitten by malaria vectors in a vector control area in southern Benin, west Africa.

BACKGROUND: The diversity of malaria vector populations, expressing various resistance and/or behavioural patterns could explain the reduced effectiveness of vector control interventions reported in some African countries. A better understanding of the ecology and distribution of malaria vectors is essential to design more effective and sustainable strategies for malaria control and elimination. Here, we analyzed the spatio-temporal risk of the contact between humans and the sympatric An. funestus and both M and S molecular forms of An. gambiae s.s. in an area of Benin with high coverage of vector control measures with an unprecedented level of resolution.
METHODS: Presence-absence data for the three vectors from 1-year human-landing collections in 19 villages were assessed using binomial mixed-effects models according to vector control measures and environmental covariates derived from field and remote sensing data. After 8-fold cross-validations of the models, predictive maps of the risk of the contact between humans and the sympatric An. funestus and both molecular M and S forms of An. gambiae s.s. were computed.
RESULTS: Model validations showed that the An. funestus, An. gambiae M form, and S form models provided an excellent (Area Under Curve>0.9), a good (AUC>0.8), and an acceptable (AUC>0.7) level of prediction, respectively. The distribution area of the probability of contact between human and An. funestus largely overlaps that of An. gambiae M form but this latter showed important seasonal variation. An. gambiae S form also showed seasonal variation but with different ecological preferences. Landscape data were useful to discriminate between the species' distributions.
CONCLUSIONS: These results showed that available remote sensing data could help in predicting the human-vector contact for several species of malaria vectors at a village level scale. The predictive maps showed seasonal and spatial variations in the risk of human-vector contact for all three vectors. Such maps could help Malaria Control Programmes to implement more effective vector control strategy by taking into account to the dynamics of malaria vector species.

Background

In response to the worldwide decrease in deaths due to malaria [1,2], malaria elimination is back on the global agenda [3]. However, evidences of reduced effectiveness of vector control interventions have been reported in some African countries [4-6], hence challenging confidence in malaria control efforts in the region. According to the authors of these studies, the diversity of malaria vector populations, expressing various resistance and/or behavioural patterns could explain this trend. A better understanding of the mosquito ecology and distribution is thus essential to take into account this diversity and to design more effective and sustainable strategies for malaria control and elimination [7-9].

The ~60 species of Anopheles which are considered as vectors of malaria [10] vary in terms of their capacity to transmit the disease because of variations in their susceptibility to the parasite (vectorial competence), host-seeking preferences (anthropophagy), and/or abundance (this latter being related to the local availability of their specific habitat) [11]. The difference in vectorial capacity can also be observed at the sub-species level. Indeed, M and S molecular forms of An. gambiae s.s. (sister taxa that are under speciation process [12]) have different ecological preferences [13-15] and susceptibility to the Plasmodium infection [16]. Regarding their ecological preferences, the S form was more frequently observed in the driest areas, in open fields that accommodate for temporary and rainfall-dependent breeding sites with low predation [13,14,17] whereas the M form showed a longer larval development [18,19] and a low sensibility to predation [20,21] highlighting a better adaptation to more permanent breeding sites. Moreover, major malaria vector species facing the wide use of residual insecticides in Africa have adapted by expressing different “strategies” including insecticide resistance [22] or behavioural avoidance [23-25]. These strategies seem to vary according to the vector population and/or the species [22,26-33].

In the Ouidah-Kpomasse-Tori Bossito (OKT) region in southern Benin, a recent cluster randomized controlled trial (RCT) did not report any benefits of using combined long-lasting insecticidal nets (LLIN) and carbamate indoor residual spraying (IRS) or Carbamate Treated Plastic Sheeting (CTPS) compared to selective coverage of LLIN [5]. In this area, the three major African malaria vectors, An. gambiae M and S forms and An. funestus are sympatric, with spatial and temporal variations [31]. High levels of physiological resistance have been recorded in both M and S forms of An. gambiae populations [34-36] whereas An. funestus seems to avoid the contact with the insecticide by modulating its biting behaviour [23]. The presence, in sympatry, of three major malaria vectors with different ecologies and resistance patterns gives a particularly interesting opportunity to study the bioecology of these vectors exposed to vector control interventions.

In the present work, we studied the spatio-temporal risk of the contact between humans and the sympatric An. funestus and both M and S molecular forms of An. gambiae s.s. according to a set of environmental covariates in 19 villages in southern Benin, with an unprecedented level of resolution. Using Binomial mixed-effect models, we assessed the environmental determinants (including vector control measures) of the probability of human-vector contact and we mapped the predicted risk of being bitten by each of these vectors during the dry and the rainy seasons.

Methods

Study area and vector control interventions

This study was carried out in the Ouidah-Kpomassè-Tori Bossito (OKT) health administrative region in southern Benin (on the Atlantic coast). Of the 58 villages screened at the baseline, 30 were excluded because they did not fulfil inclusion criteria i.e. distance between two villages greater than two kilometres, population size between 250 and 500 inhabitants with non-isolated habitations, and absence of any local health care centre [5]. Twenty-eight villages were randomly assigned to four groups (seven villages by groups) for implementation of four vector control interventions: (1) targeted LLIN coverage (TLLIN) to pregnant women and children younger than six years; (2) universal LLIN coverage of sleeping units (ULLIN), (3) targeted LLIN coverage plus full coverage of carbamate IRS (TLLIN+IRS), and (4) universal LLIN coverage plus full coverage of carbamate CTPS lined up to the walls of the household (ULLIN+CTPS). Of the 28 villages used during the clinical trial, nine (three TLLIN, two ULLIN, two ULLIN+CTPS, and two TLLIN+IRS) were excluded because they were not covered by the satellite image used for the landscape analysis. Details about the allocation of vector control intervention methods have been described elsewhere [5].

Mosquito collections

Mosquitoes were collected every six weeks during the year 2009 (eight surveys) using the human landing catch (HLC) technique. HLC were carried out from 22:00 to 06:00 both indoors and outdoors of four houses in each village on two successive nights for each survey (i.e., 16 collector-nights per village per survey). Sites were a minimum distance of 50 metres from each other and were homogeneously distributed across the village (sites situated near eucalyptus trees, smoke, etc. were discarded). Collectors were rotated hourly between collection sites and/or position (indoor/outdoor). Independent staff supervised rotations and regularly checked the quality of mosquito collections.

Malaria vectors collected on humans were identified using morphological keys [37,38]. All mosquitoes belonging to the Gambiae complex and the Funestus Group were identified to species by PCR [39,40]. Both M and S molecular forms of An. gambiae s.s. were identified by the method of Favia et al. [41].

Peri-domestic breeding sites inventory

Anthropogenic larval habitats inside the village borders were surveyed using a standard dipping method [42]. Sampling was repeated during each survey of adult collection and the number of breeding sites positive (or not) for the presence of Anopheles sp. larvae was recorded. The Breteau index, representing the number of positive containers for larvae per 100 houses was estimated in each locality [43].

Size of the buffer zone

To characterize the environment in the neighbourhood of the 19 villages, buffer zones were defined around each village. According to Service [44], two kilometres is the maximum flight range for Anopheles sp. and breeding sites located beyond that distance can be considered as insignificant. The distance of two kilometres was therefore selected as the radius of the buffer zones.

Spatio-temporal data

We used the Land-Surface Temperature (LST) at a spatial resolution of one kilometre measured by the Moderate Resolution Imaging Spectroradiometer (MODIS) sensors on the Terra satellite (https://lpdaac.usgs.gov/). The weekly nocturnal and diurnal temperatures during the week of the surveys and during each of both weeks preceding the surveys were extracted using ArcGis ArcInfo 9.3 software (ESRI, Redlands, CA) in the buffer zone of each village. When the data were not available at a specific date, we used the TiSEG software [45] to estimate missing data: according to the quality statement of the pixels provided as metadata with MODIS data, pixels that do not fit a "good" or an "acceptable" quality (because of clouds for example) were temporally interpolated (linear interpolation) between the values of the same pixels (if they fitted the required quality) of the closest preceding and following images.

We used the Normalized Difference Vegetation Index (NDVI) at a spatial resolution of 250 metres measured by the MODIS sensor. For each survey and each village, the average 16-day NDVI during the two-week period including the survey and the two-week period preceding the survey was extracted in the buffer zone. TiSEG was used to interpolate missing data using the same method that was used with LST data.

Daily rainfall data used was derived from the Tropical Rainfall Measuring Mission (TRMM) and other Satellite Precipitation products (3B42 Version 6). TRMM daily rainfall data were cumulated for the 15 days preceding each mosquito collection and spatially interpolated (using a kriging technique). The mean estimate of the cumulated rainfall in the buffer zone of each village was extracted. Cumulated rainfall was expressed in millimetres. The same procedure was used to extract the number of rainy days over the 15 days preceding each mosquito collection in each village.

Spatial only data

A Digital Elevation Model (DEM) derived from the Satellite Pour l'Observation de la Terre (SPOT) SPOTDEM product with a 30 meters resolution was used to compute mean elevation, mean slope, count of sink and mean theoretical flow accumulation (i.e. for each pixel, the surface of its drainage area giving the theoretical drainage network capacity). To describe the ability of soils to accommodate for breeding sites, hydromorphic soils were identified on a digitized soil map of south Benin [46]. Hydromorphic soils are poorly-drained and characterized by an excess of water due to temporary waterlogging or the presence or rise of groundwater [47,48]. The percentage of area of hydromorphic soils in the buffer areas was computed.

A land-cover map was obtained by carrying out a supervised object-oriented nearest-neighbour classification [49] of a SPOT-5 satellite image (pixel of five metres) acquired on 01/23/2010 (eCognition™ software, Definiens-imaging, Munich, Germany). More than 140 plots of known land-cover were identified and geolocated in the field. These plots were used as training plots to define the feature of each land-cover class and then classify the overall image using a nearest-neighbour algorithm. The accuracy of the final classification was assessed using a confusion matrix to compare the allocated land cover class to 70 supplementary ground-truth validation plots.

A landscape analysis was conducted on the land-cover map to characterize its spatial structure. Using the Fragstats freeware [50], the following metrics were calculated for each class of land-cover in the buffer zone of each village: the percentage of the surface, the number of patches, the edge length. Other metrics were used to describe the diversity present within the whole landscape (i.e. for the entire buffer regardless of the class): the Patch Richness Density, the Simpson's Diversity Index and the Modified Simpson's Eveness Index (all metrics are defined in the Fragstats documentation [50]).

Market gardening could create suitable breeding sites due to the irrigation and water storage. However, these activities are difficult to discriminate on SPOT imagery because of the small size of fields. We therefore carried out field surveys (in March 2011) to inventory market gardening areas. Roads were digitalized from 1/50000 maps from 1968 available at the Institut Géographique National of Benin and updated using the SPOT-5 image.

In order to describe attractiveness and penetrability for malaria vectors, village perimeters were extracted by on screen digitalization using visual interpretation of the SPOT image. Then, area and population density were computed and the distance from each collection site to the perimeter of the village was measured. The number of neighbourhoods (homogenous groups of houses separated from other groups by a vegetated strip) was recorded to describe the clustered or spread-out layout of the villages. Using the Shape Metrics Tool (http://clear.uconn.edu/tools/Shape_Metrics/), the normalized Spin and Depth indices describing the shape of the villages were computed. Cattle farms were inventoried because they can interfere with host seeking behaviour of malaria vectors that may bite on cattle rather than humans [51].

Statistical analysis

All analyses were performed using the 'R' software [52] and the additional 'lme4' [53] and 'pROC' [54] packages.

The probability of human-vector contact (PHVC) for each species and molecular forms considered during a survey (two consecutive nights) was assessed using a Binomial Mixed-Effect (BME) model with nested random effects at the village and collection site level. The dependent variable was the presence/absence of vectors collected during 1, 216 catches (19 villages * four collection sites * two places (indoors and outdoors) * eight surveys). The BME was adjusted for the vector control intervention and the collector's position (indoor or outdoor). The independent variables were tested both as continuous and categorical variables (after recoding). They were kept as continuous variables whenever possible, and were stratified into terciles (or into presence/absence if more relevant) when necessary to accommodate for non-monotonic and non-linear responses. Only variables found to be significantly correlated (p<0.05) with the PHVC of malaria vector species were kept for analysis of collinearity. Collinear covariates were deleted based on empirical knowledge until the Variance Inflation Factor (VIF) of remaining variables was below three [55].

All selected variables from the univariate analysis were introduced in the multivariate BME models, and a backward procedure was applied to select only those that remained significant with a p-value<0.05 in the final model.

The structure of the models was evaluated by 8-fold cross-validation with the data of each survey successively used for validation. The predictive accuracy of the models was assessed using the ROC (Receiver Operating Characteristic) curve method [56]. The area under the curve (AUC) of the ROC curve, its 95% confidence interval (95% CI), sensitivity and specificity were calculated to assess the quality of the models [54].

Risk maps of human-vector contact

Based on the final multivariate BME models, two seasonal maps of the risk of contact with the three species were computed for the 15-16/01/2009 (dry season) and the 30-31/06/2009 (rainy season). Covariates describing attractiveness and penetrability for which data were not available at all points of the area were set at a constant value equal to the mean calculated for overall villages.

Ethics statement

The IRD (Institut de Recherche pour le Développement) Ethics Committee and the National Research Ethics Committee of Benin approved the study (CNPERS, reference number IRB00006860). All necessary permits were obtained for the described field studies. No mosquito collection was carried out without the approval of the head of the village, the owner and occupants of the collection house. Mosquito collectors gave their written informed consent and were treated free of charge for malaria presumed illness throughout the study.

Results

Malaria vector collection

In the 19 villages (Table 1) during 2,432 human-nights of collection, 2,379 malaria vectors were collected: 1,091 were An. funestus and 1,288 were An. gambiae s.s. among which, 1,063 belonged to the molecular M form and 225 belonged to the molecular S form. Over 1,216 catches (two human-nights per site per survey) were used for modelling, only 252, 323, and 114 were positive for An. funestus, An. gambiae M form, and An. gambiae S form, respectively.

Table 1

Spatial coordinates of the studied villages and entomological data

	Village center coordinates		Proportion (%) of catches positive for			Total no. of vectors collected
Village name	Latitude	Longitude	An. funestus	An. gambiae M form	An. gambiae S form	An. funestus	An. gambiae M form	An. gambiae S form
Agouako	N 06°30'30.2"	E 002°09'47.7"	1.56	9.38	6.25	1	7	5
Amoulehoue	N 06°26'36.3"	E 002°07'19.9"	35.94	26.56	4.69	50	53	8
Assogbenou daho	N 06°27'19.5"	E 002°05'24.5"	4.69	18.75	6.25	3	22	6
Ayidohoue	N 06°31'35.8"	E 002°05'59.6"	1.56	6.25	0.00	1	8	0
Adjahassa	N 06°33'11.3"	E 002°09'34.5"	0.00	18.75	6.25	0	27	4
Dokanmey	N 06°33'13.9"	E 002°13'37.5"	3.13	10.94	4.69	2	9	5
Kindjitokpa	N 06°25'37.5"	E 002°07'31.3"	84.38	39.06	6.25	297	127	6
Tokoli-vidjinnagnimon	N 06°26'57.1"	E 002°09'36.6"	65.63	34.38	12.50	230	81	10
Hekandji	N 06°29'26.0"	E 002°05'51.0"	0.00	4.69	6.25	0	3	5
Hla	N 06°32'34.0"	E 002°11'54.8"	3.13	21.88	10.94	2	41	7
Hounkponouhoue	N 06°33'31.4"	E 002°05'27.7"	1.56	7.81	15.63	2	8	12
Manguevier	N 06°28'21.8"	E 002°09'29.3"	18.75	53.13	4.69	18	95	4
Tanto	N 06°33'13.1"	E 002°10'48.0"	0.00	29.69	23.44	0	88	64
Lokohoue	N 06°24'24.2"	E 002°10'32.1"	71.88	53.13	7.81	273	175	10
Todo	N 06°30'09.9"	E 002°08'09.2"	0.00	26.56	6.25	0	56	6
Tokoli-dozouzrame	N 06°27'35.2"	E 002°10'17.1"	50.00	43.75	25.00	112	100	22
Wanho	N 06°27'32.9"	E 002°11'19.6"	46.88	62.50	0.00	96	77	0
Agadon	N 06°33'35.8"	E 002°06'25.7"	3.13	23.44	25.00	2	67	43
Zoume	N 06°32'42.7"	E 002°06'50.5"	1.56	14.06	6.25	2	19	8

Latitudes and longitudes are given in the standard coordinate system WGS84. Proportion of catches positive for each species was calculated over 64 catches carried out in each village (4 collection sites * 2 places (indoors and outdoors) * 8 repetitions (surveys)).

Land-cover map

Fourteen land-cover classes were discriminated: surface freshwater (1.8% of the total surface), aquatic grassland (1.4%), herb swamp (2.1%), coco tree (1%), eucalyptus tree (4.5%), palm tree (11.5%), teak tree (6.3%), pineapple (1.8%), rainfed agriculture (28%), forest (6%), degraded riparian forest (3%), thicket (13.9%), savanna (5.9%), and degraded surface (bare soil and constructed areas, 12.4%). The confusion matrix revealed that the overall accuracy of the classification was excellent (98.9%). The resulting land-cover map of the study area is presented in Figure 1.

Figure 1

Land-cover map based on a supervised object-oriented classification of a SPOT-5 satellite image acquired on 01/23/2010. This map was obtained by carrying out a supervised object-oriented nearest-neighbour classification using the software eCognition™ (Definiens-imaging, Munich, Germany). More than 210 plots of known land-cover were identified and geolocated in the field to be used as training and validation (ground-truth) data.

Determinants of the probability of human-vector contact

Covariates that were kept in final multivariate models, their odds-ratio and 95% confidence intervals are presented in Tables 2 (An. funestus), Table 3 (An. gambiae M form), and Table 4 (An. gambiae S form).

Table 2

Multivariate binomial mixed-effect model of the risk of human-contact

		OR	95% CI		P-value
Surface water	absence	1
	presence	228.775	81.888	639.144	< 2e-16	***
Area of hydromorphic soils (per additional 100 ha)		3.823	2.136	6.843	< 5e-4	***
NDVI 2 weeks before the catch		0.001	0.000	0.067	0.001	***
Diurnal temperature 2 weeks before the catch (per additional °C)		1.321	1.175	1.485	< 5e-4	***
Nocturnal temperature (1 week before the catch; in °C)	<20.68	1
	20.68-21.64	0.090	0.040	0.200	< 5e-4	***
	≥21.64	0.129	0.053	0.316	< 5e-4	***
Cumulated precipitation 16 days preceding the catch (per additional mm)		1.012	1.009	1.016	< 5e-4	***
Number of neighbourhoods	<2	1
	2	2.187	0.767	6.233	0.143
	≥3	116.273	35.779	377.856	< 5e-4	***
Vector control intervention	TLLIN	1
	ULLIN	0.679	0.311	1.481	0.330
	ULLIN+CTPS	0.264	0.108	0.646	0.004	**
	TLLIN+IRS	0.398	0.151	1.045	0.061	.

OR: Odds-ratio; CI: Confidence Interval; ha: hectares; °C: degrees Celsius; mm: millimetres; m: metres; TLLIN: Targeted distribution of Long-Lasting Insecticidal Nets; ULLIN: Universal distribution of LLIN; CTPS: Carbamate Treated Plastic Sheeting; IRS: carbamate Indoor Residual Spraying.

Table 3

Multivariate binomial mixed-effect model of the risk of human-M form contact

		OR	95% CI		P-value
Proportion of unvegetated land (per additional %)		1.080	1.047	1.113	< 5e-4	***
Number of patches of unvegetated lands	<75	1
	75-92	6.321	3.472	11.509	< 5e-4	***
	≥92	1.174	0.421	3.275	0.759
Herbswamp	absence	1
	presence	0.276	0.135	0.565	< 5e-4	***
Area of hydromorphic soils (per additional 100 ha)		20.323	9.908	41.688	< 2e-16	***
Domestic breeding sites	absence	1
	presence	2.021	1.215	3.360	0.007	**
NDVI (during collection)	<0.591	1
	0.591-0.652	2.472	1.427	4.281	0.001	**
	≥0.652	3.042	1.643	5.635	< 5e-4	***
Diurnal temperature 1 week before the catch (per additional °C)		0.735	0.614	0.882	0.001	***
Nocturnal temperature the week of the catch (per additional °C)		0.755	0.651	0.875	< 5e-4	***
Nocturnal temperature 2 weeks before the catch(per additional °C)		0.835	0.714	0.977	0.024	*
Number of days with rainfall during the 16 days preceding the catch (per additional day)		0.842	0.751	0.943	0.003	**
Cumulated precipitation during the 16 days preceding the catch (per additional mm)		1.020	1.015	1.025	< 5e-4	***
Number of cattle (per additional individual)		0.990	0.981	0.998	0.021	*
Number of neighbourhoods	<2	1
	2	0.779	0.380	1.595	0.495
	≥3	10.644	3.114	36.386	< 5e-4	***
Vector control intervention	TLLIN	1
	ULLIN	3.214	1.241	8.324	0.016	*
	ULLIN+CTPS	1.902	0.820	4.416	0.134
	TLLIN+IRS	0.459	0.170	1.239	0.124

OR: Odds-ratio; CI: Confidence Interval; ha: hectares; °C: degrees Celsius; mm: millimetres; m: metres; TLLIN: Targeted distribution of Long-Lasting Insecticidal Nets; ULLIN: Universal distribution of LLIN; CTPS: Carbamate Treated Plastic Sheeting; IRS: carbamate Indoor Residual Spraying.

Table 4

Multivariate binomial mixed-effect model of the risk of human-S form contact

		OR	95% CI		P-value
Road length (m)	<18,000	1
	18,000 - 22,000	6.272	2.321	16.949	< 5e-4	***
	≥22,000	2.978	1.202	7.380	0.018	*
Number of patches of pineapple	<10	1
	10-25	5.668	1.135	28.309	0.035	*
	≥25	1.581	0.664	3.763	0.301
Elevation (per additional m)		1.058	1.028	1.088	< 5e-4	***
Slope (per additional degree)		0.551	0.291	1.042	0.067	.
Hydromorphic soil	absence	1
	presence	13.407	4.705	38.206	< 5e-4	***
Nocturnal temperature 2 weeks before the catch(per additional °C)		1.330	1.086	1.628	0.006	**
NDVI 2 weeks before catch		0.016	0.000	0.639	0.028	*
Cumulated precipitation 16 days preceding the catch (per additional mm)		1.005	1.002	1.009	0.003	**
Number of breeding sites per 100 houses (per additional site)		1.055	0.999	1.114	0.056	.
Domestic breeding sites	absence	1
	presence	2.668	1.485	4.792	0.001	**
Vector control intervention	TLLIN	1
	ULLIN	2.799	1.016	7.709	0.046	*
	ULLIN+CTPS	1.016	0.365	2.824	0.976
	TLLIN+IRS	0.615	0.152	2.483	0.495
Collection site	Indoor	1
	Outdoor	1.519	0.992	2.326	0.054	.

OR: Odds-ratio; CI: Confidence Interval; ha: hectares; C: degrees Celsius; mm: millimetres; m: metres; TLLIN: Targeted distribution of Long-Lasting Insecticidal Nets; ULLIN: Universal distribution of LLIN; CTPS: Carbamate Treated Plastic Sheeting; IRS: carbamate Indoor Residual Spraying.

Cumulated rainfall during the 15 days preceding the catch positively correlated to the PHVC of each vector species (p< 0.001), whereas the number of rainy days was negatively associated with the PHVC of An. gambiae M form. NDVI was correlated negatively with the PHVC of An. gambiae S form and An. funestus (NDVI two weeks before collection) but positively with An. gambiae M form (NDVI of the collection period).

The PHVC of the three species was explained by temperature data recorded two weeks before the catch. The PHVC of An. funestus and An. gambiae M form decreased respectively with the diurnal and the nocturnal temperatures recorded during the week or the week preceding the mosquito collection. They were also both positively correlated with the number of neighbourhoods.

PHVC of M and S forms of An. gambiae was positively correlated with the presence and the number of domestic breeding sites. Interestingly, the PHVC of the An. gambiae M form was associated with both the surface and an intermediate number of patches of unvegetated land but negatively correlated with the presence of herbswamp. The PHVC of An. gambiae M form decreased with the number of cattle recorded around the village.

The PHVC of An. gambiae S form was positively associated with the elevation, the length of roads, an intermediate number of patches of pineapple plantation (between 10 and 25) but negatively correlated with the slope. The presence of hydromorphic soils correlated with the PHVC of the three species and the surface of freshwater increased the PHVC of An. funestus.

The PHVC of An. gambiae S form was lower in the villages where ULLIN+CTPS were implemented compared to villages belonging to the reference group (TLLIN). However, the PHVC of An. gambiae M form or S form was higher in villages where ULLIN was implemented. The PHVC of An. gambiae S form was significantly higher outdoors than indoor.

Validation

The ROC curves of the 8-fold cross validation are presented in Figure 2. The AUC of the An. funestus BME was 0.92 (95% CI 0.90-0.94; Figure 2A) indicating an excellent prediction. For the probability threshold 0.13, specificity was 80% (95% CI 78–83) and sensitivity was 93% (95% CI 90–96).

Figure 2

ROC curve of the predictive model of the contact between human and (A) (B) M form and (C) S form. ROC (Receiver Operating Characteristic) curves were performed based on 8-fold cross-validations of the final multivariate models with the data of each survey successively used for validation. The curves plot the sensitivity versus one minus the specificity. Greater is the area under the curve, greater is the predictive accuracy of the model.

The AUC of the An. gambiae M form BME was 0.84 (95% CI 0.81-0.86; Figure 2B) indicating a "good" prediction. For the probability threshold 0.21, specificity was 65% (95% CI 63–69) and sensitivity was 86% (95% CI 83–90). The AUC of the An. gambiae S form BME was 0.77 (95% CI 0.72-0.82; Figure 3C) indicating an "acceptable" prediction. For the probability threshold 0.12, specificity was 81% (95% CI 79–83) and sensitivity was 65% (95% CI 56–74).

Figure 3

Risk maps of the contact between human and (A, B) , (C, D) M form and (E, F) S form for 2 nights during the dry and the rainy season. Risk maps were computed based on the final multivariate binomial mixed-effect models. Two seasonal maps of the probability of contact with the three species/forms were produced for the 15-16/01/2009 (dry season) and the 30-31/06/2009 (rainy season). Covariates describing attractiveness and penetrability for which data were not available at all points of the area were set at a constant value equal to the mean calculated for overall villages.

Risk maps

Maps of the predicted PHVC of An. funestus, An. gambiae M form, and An. gambiae S form for two nights during the dry and the rainy season are presented in Figure 3. During the dry season, An. funestus was highly concentrated around Toho Lake in the southern part of the study area. During the rainy season, we observed an increase of its PHVC along the downstream parts of the rivers that drain into Toho Lake and the lagoon. The PHVC of An. gambiae M form followed a similar distribution pattern in the dry season: areas of highest risk were found in the southern part of the study area, largely overlapping areas of high PHVC of An. funestus. Distribution of the PHVC of An. gambiae M form showed the most important seasonal variations spreading along the rivers and in the most elevated part of the study area during the rainy season. The PHVC of An. gambiae S form was lower than that of An. funestus and An. gambiae M form. Its distribution is very different as it is largely confined in the northern and eastern part of the study area. It represented a risk in the south only during the rainy season.

Discussion

The M and S molecular forms of An. gambiae s.s. and An. funestus occurred in sympatry in the OKT region [31] and this setting is particularly favourable to study the environmental determinants involved in the spatial and temporal distribution of each malaria vector species. Here, we fitted specific models at an unprecedented level of resolution to predict the spatio-temporal risk of the human-vector contact in 19 villages in southern Benin. The models were validated using an 8-fold cross-validation technique. The predictions given by BME models were considered as "good" and "excellent" for An. gambiae M form and An. funestus respectively, whereas the predictions for the An. gambiae S form model was "acceptable". The lower predictive value of the An. gambiae S form model was likely due to the lower number of mosquitoes collected during the collection period. Indeed An. gambiae S form was previously described as being at the edge of its distribution area in southern Benin [57]. To our knowledge, this study is the one considering the most important number of environmental covariates likely to impact on the ecology of these malaria vectors. Indeed, we paid great attention to take into account most of the relevant factors driving the creation of breeding sites, the vector dispersion, and the vector survival. However, other factors (e.g. human settlements other than the studied villages, predation, personal protection, humidity, wind, etc.) affecting the contact between human and malaria vectors were not included in the analysis because of unavailability, poor quality, or low resolution.

The present study seems to indicate that An. gambiae M form was not dependent on the presence of permanent wetlands (freshwaters and aquatic grassland). This result contrasts with previous studies conducted in Burkina Faso and Cameroon [13,14]. Instead, we found a positive correlation between the PHVC of An. gambiae M form and both the surface and the number of patches of unvegetated land. This might be explained by the fact that the SPOT image used to assess land use was acquired during the dry season when seasonal wetlands were almost absent. The PHVC of An. gambiae M form was negatively correlated to the number of rainy days hence suggesting that this species may be particularly sensitive to flush-out. Nevertheless, all species were positively associated with the surface and presence of hydromorphic soils that might provide suitable semi-permanent breeding sites during the rainy season. As expected, the PHVC of An. funestus was predicted by the presence of freshwater bodies that could provide ideal breeding sites for this species [58].

Length of roads in the buffer area around the villages positively predicted the PHVC of An. gambiae S form. Indeed, road ditches and puddles created by vehicles are known to provide good breeding sites for this species, particularly during the rainy season [59]. Interestingly, only the intermediate class of number of patches (10 to 25) of pineapple plantations increased the probability of contact with An. gambiae S form. Leaf axils of pineapple plants could provide breeding sites for An. gambiae as previously observed in the neighbouring Nigeria [60] and regular traffic of vans carrying fruits could create ruts in roads. However, the relationship was not linear and not related to the surface indicating that external factors probably alter the development of An. gambiae S form when the number of plantations increases. We presume that different species of pineapples, cropping techniques or the presence of more or less predators may impact on this trend. Despite a small range of elevations and slopes in the study area, these factors were significantly associated with the PHVC of An. gambiae S form as described elsewhere in Africa [13-15]. An. funestus and An. gambiae M form were correlated to the number of neighbourhoods that is an indicator of the villages' layout. This may indicate that scattered habitations in a village might increase the attractiveness of the village for these species. This suggests that the field of attraction of a scattered village may be larger than that of a clustered village of the same size.

Temperatures recorded two weeks before the catch (corresponding to the larval period) explained the PHVC of the three species. The PHVC of An. funestus and An. gambiae S form was higher when temperatures increased. A different scenario was observed with An. gambiae M form: the lower the temperature (over different time periods comprised between the week of the catch and two weeks before the catch), the higher the PHVC of An. gambiae M form. This result could indicate that the M form was more susceptible to high water temperatures that can induce high mortality rates at the larval stage [61-63] and more active when temperatures were the lower. However, there is no published data to support these assumptions and the real underlying reason for these associations remains unclear. Our models show that the relationship between temperature and distribution of malaria vectors is complex. Indeed, both diurnal and nocturnal temperatures, over different time periods (comprised between the week of the catch and two weeks before the catch), were significantly predictive for two of the three species studied. Temperature is known to influence many aspects of the mosquito lifecycle such as larval growth, survival of adults and biting activity. Moreover, our team found a negative correlation between elevated temperature and LLIN use in Benin that could also impact on mosquito longevity and density [64]. Further investigations are needed to better understand the impact of temperature on the ecology of these vector species.

Our results indicate that the PHVC of the M and S molecular forms of An. gambiae s.s. was greater in presence of domestic breeding sites. This is worrying because it might indicate that villagers “bred” malaria vectors in their villages and inside their houses and might therefore support part of the malaria transmission. Since the majority of containers involved as breeding sites were used for water storage, these breeding sites were not rain-dependent and remained full especially during the dry season. No information is available yet regarding the Anopheles species present in the domestic breeding sites. Larviciding may then represent a complementary strategy for targeting malaria vectors within semi-urbanized population centres in rural areas [65].

Regarding vector control interventions, the probability of contact between humans and An. funestus was reduced only in villages sharing the ULLIN+CTPS combined intervention. This trend might reflect a decrease in the density of An. funestus following the implementation of vector control intervention hence underlying a certain efficacy of the treatments. In contrast, the PHVC of both molecular forms of An. gambiae s.s. was higher in villages where ULLIN were implemented alone. This finding contrasts with that of Corbel and colleagues [5] as they did not find a significant difference between ULLIN and TLLIN in terms of mosquito density and malaria transmission. We assume that mosquito collectors were overexposed to the vector bite. Indeed, high levels of LLIN use as observed in the ULLIN arm [5] might have contributed to divert a part of the vector population to alternatives hosts (like cattle, as suggested by the An. gambiae M form model) or to unprotected hosts such as our collectors [51,66].

Conclusion

To conclude, the predictive maps showed seasonal and spatial variations of the risk of human-vector contact for the all three vectors. These variations were characterized by an increase in the distribution areas of the three species during the rainy season. Such maps could help Malaria Control Programs to implement more effective vector control strategies by taking into account the dynamics of malaria vector species. The use of larvicides in domestic breeding sites might be a cost effective strategy for malaria control [67,68] considering the possible ecology of An. gambiae s.s. in rural areas in Benin.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

NM, VC, FC and HG designed the study. NM, ASB, and AD collected field data. NM treated the data. NM and HG performed the analysis. NM, VC and HG wrote the article. All authors read and approved the final version of the article.