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 . 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 .
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. .
Peri-domestic breeding sites inventory
Anthropogenic larval habitats inside the village borders were surveyed using a standard dipping method . 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 .
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 , 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.
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  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 . 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  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 , 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 ).
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 .
All analyses were performed using the 'R' software  and the additional 'lme4'  and 'pROC'  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 .
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 . 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 .
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.
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.