Background {#Sec1}
==========

Mosquitoes (Diptera:Culicidae) are known to be vectors of a large number of pathogens around the globe and are considered as prime candidates for transmitting (re-)emerging vector-borne diseases (VBDs) in Europe \[[@CR1]\]. The increased mobility of humans, that has also increased the mobility of livestock and pathogens, as well as environmental modifications and climate changes can contribute to the (re-)emergence of vector-borne diseases \[[@CR2]\]. Furthermore, mosquito bites can cause a considerable nuisance for humans and mammals. Severe nuisance can have negative economic consequences (e.g., in tourism, work productivity outdoors, meat and dairy production) \[[@CR3]\]. These nuisance situations can eventually lead to autochthonous VBD cases, when in non-endemic areas infectious reservoirs, either humans (travellers, temporary workers) or animals (livestock, migrating animals) come in contact with high density of mosquito vectors.

In 2010, in the Netherlands, the Centre for Monitoring of Vectors (CMV) started a nationwide inventory of indigenous mosquitoes to acquire basic information on the composition, geographical distribution, biodiversity and environmental preferences of mosquito species. In this survey, the nuisance mosquito species and potential malaria vector *Anopheles plumbeus* (Stephens, 1828) (Figure [1](#Fig1){ref-type="fig"}) was one of the most frequently collected mosquito species \[[@CR4]\]. This species has been proven to be able to transmit *Plasmodium falciparum* (Welch, 1897), the causative agent of malaria tropica \[[@CR5]\]. Circumstantial evidence for local transmission of *P. falciparum* malaria by *An. plumbeus* has been reported for Germany \[[@CR6]\]. *Anopheles plumbeus* has also been incriminated as a vector of *Plasmodium vivax* (Grassi & Feletti, 1890) \[[@CR7],[@CR8]\] and has been proven to be a laboratory vector of West Nile virus \[[@CR9]\]. Even though the health care system is likely to rapidly identify malaria patients and thereby prevent the building up of an infectious human reservoir of *Plasmodium* parasites \[[@CR10]\], it is important to gain information on the spatial distribution of *An. plumbeus*, in order to inform the health care system on the areas at risk.Figure 1*Anopheles plumbeus* female (source: A. Ibañez-Justicia).

*Anopheles plumbeus* is a mosquito species commonly found in forests, where larvae are usually found in water in rot-holes of trees with high salinity and deficiency of oxygen \[[@CR11]\]. They can also be found in containers with stagnant rain water and groundwater, such as tyres, rainwater casks and cemetery vases \[[@CR5],[@CR12]\]. In the last decade, this species has also been associated with abandoned stables where it breeds in the rain water collected in the manure cellars \[[@CR13]\]. This species is known to be a particularly aggressive biter, feeding at any time of the day on different mammalian hosts (including humans), and to a lesser degree on birds and reptiles \[[@CR14]\]. In June 2006, nuisance caused by *An. plumbeus* was reported for the first time in the Netherlands, near the city of Nijmegen \[[@CR1]\]. Since then, *An. plumbeus* nuisance have been reported every year in the Netherlands, mostly in proximity to abandoned pig stables (Ibanez-Justicia, unpublished).

An understanding of the spatial extent of potential vector species, their abundance and seasonal activity, is important for estimating levels of risks of VBDs and enabling better targeting for surveillance and control. In order to develop basic reproduction number (R0) models and construct risk maps that indicate the risk for an outbreak after an introduction, abundance data of vectors are an essential parameter \[[@CR15],[@CR16]\]. Although vector presence and abundance are not the only factors determining whether or not a pathogen can spread in an area, determining the distribution of the vector is an essential step in studying the risk of transmission of a pathogen. Given the nuisance and potential risk for the human health, such information on *An. plumbeus* is needed. Currently, no information on the potential spatial distribution of this species is available for the Netherlands.

In this study, we modelled the potential spatial distribution, expressed in occurrence (predicted probability of presence) and abundance of *An. plumbeus* in the Netherlands, based on data collected during the National Mosquito Survey and environmental data. The occurrence was modelled to predict the environmental suitability of the species using a random classification forest model. The abundance was modelled using a random regression forest model with the aim to identify areas where mosquito peaks could be expected. Random (classification and regression) forest models allow external validation through a bootstrapping procedure. The occurrence model was validated also with an external dataset. The resulting maps are in agreement with the reported nuisance for this species and the predictions show a good matching with an external dataset used to validate the model.

Methods {#Sec2}
=======

Species distribution modelling links the occurrence or the abundance of species with environmental data and estimates the similarity of the conditions at any site based on the conditions at the locations of known occurrence/abundance of a species. Here we describe the mosquito data collection, the environmental data used and the statistical methods applied in this study.

Mosquito data {#Sec3}
-------------

Mosquito data used for the modelling were obtained from the national mosquitoes survey that was carried out from April to October 2010--2013 by the Dutch Centre for Monitoring of Vectors. Mosquitoes were captured using CO~2~ baited Mosquito Magnet Liberty Plus MM3100 (Woodstream® Co., Lititz, USA). Traps were randomly distributed in the country following a cross-sectional study design, with the following constraint: 40% of the traps were placed in urban areas, 40% in rural-agricultural areas and 20% in natural areas \[[@CR4]\]. Urban and agricultural areas were sampled more intensively, because of the potential higher human and veterinary health risk in those areas due to higher exposure.

Data consisted of mosquito abundance data, sampled at 778 locations. For this study the abundance data were reclassified into data of presence (when at least one mosquito was found in the trap) and absence (when no mosquitoes were found in the trap). Each of the locations was sampled only once and each trap was active for one week. The content of the traps was collected weekly and sent to the CMV laboratory. In the laboratory, mosquitoes were morphologically identified using, among others, the Culicidae key specifically designed for rapid field-identification of Dutch adult Culicidae (modified key after Snow \[[@CR17]\], Schaffner et al. \[[@CR9]\], Verdonschot \[[@CR18]\], Becker et al. \[[@CR19]\]). Twenty-seven mosquito species were found in the National Mosquito Survey and *An. plumbeus* was the 7th mosquito species most commonly found in the Netherlands. This species was active in the whole period of the survey, from April until October \[[@CR4]\]. When a presence and an absence point were in the same square kilometre only the presence point was used because presences inform about the places that are environmental suitable for a species, but absences do not necessarily indicate the opposite \[[@CR20]\]. This reduced the number of locations used in the analysis from 778 to 766.

For the validation of the occurrence model, data on *An. plumbeus* presence from confirmed nuisance notifications and data from other mosquito surveys carried out by the CMV in the Netherlands during the years 2010--2014 were used (Table [1](#Tab1){ref-type="table"}). The mosquito data from these surveys were collected with various trapping methods: dippers, pooters, CDC miniature light traps Model 512 (John W. Hock Company, Gainesville, USA), BG Sentinel traps (Biogents AG, Regensburg, Germany) and Mosquito Magnet traps. These data were extracted from VecBase, a tailor-made application built for CMV in 2010 for vector surveillance data.Table 1**Surveys used for the validationSurvey nameYearSampling strategyCapturing deviceTotal nr locations**EMS-Used tires2010-2014Target longitudinal samplingLarval sampling, manual aspirator, BG-Sentinel trap16EMS-Lelystad2013Target samplingMM-Liberty Plus trap3NVS-Limburg2009Cross-sectionalMM-Liberty Plus trap14NVS-Mosquitoes longitudinal2011Target longitudinal samplingMM-Liberty Plus trap1Projects2011, 2012Target longitudinal samplingCDC light trap, manual aspirator3Nuisance2010, 2011, 2013, 2014Check at locations of reported nuisanceLarval sampling, manual aspirator6West-Nile-Virus Wetlands2010Target longitudinal samplingMM-Liberty Plus trap, CDC light trap2The predictions obtained with the occurrence model that used National Mosquito Survey data were validate with data from these surveys.EMS: Exotic Mosquito Survey.NVS: National Vector Survey.

Environmental variables {#Sec4}
-----------------------

The environmental data included in the analysis as predictor variables are 1 km^2^ resolution satellite images and meteorological data in raster file format, commonly used for mosquito distribution modelling \[[@CR21]\]. The images were obtained from the MODIS sensor on NASA's Terra and Aqua satellites \[[@CR22],[@CR23]\] for 2000--2012 and subjected to temporal Fourier transformation \[[@CR24],[@CR25]\] to summarise the images and to produce sets of data that capture characteristics of the annual seasonality: the mean, the annual, bi-annual and tri-annual amplitudes and phases, the maxima, minima and variances of the middle infra-red, day and night-time land surface temperature, the enhanced vegetation index and the normalized difference vegetation index signals \[[@CR26]\]. Other environmental data used in this study are precipitation (WorldClim \[[@CR27]\] and CMORPH \[[@CR28]\] 1950--2000), population density (compiled from the Gridded Population of the World Dataset 2000 \[[@CR29]\]), the digital elevation model (MODIS \[[@CR23]\] 2012) and land cover (Corine land cover map of 2006 \[[@CR30]\]). The Fourier components used are provided in Table [2](#Tab2){ref-type="table"} and the environmental data in Table [3](#Tab3){ref-type="table"}. For each trap location the pixel values of the environmental variables were extracted.Table 2**Fourier components from temporal Fourier analysis of an imagery time seriesComponentDescription**A0Fourier mean for entire time seriesMNMinimum valueMXMaximum valueA1Amplitude of annual cycleA2Amplitude of bi-annual cycleA3Amplitude of tri-annual cycleVRTotal varianceP1Phase of annual cycleP2Phase of bi-annual cycleP3Phase of tri-annual cycleD1Proportion of total variance due to annual cycleD2Proportion of total variance due to bi-annual cycleD3Proportion of total variance due to tri-annual and cycleDAProportion of total variance due to all three cyclesComponent is the name used in Vecmap.Table 3**Environmental predictor variablesSourceVariables**MODISMiddle Infra-red (MIR)MODISDay-time land surface temperature (DLST)MODISNight-time land surface temperature (NLST)MODISEnhanced vegetation index (EVI)MODISNormalised difference vegetation index(NDVI)CMORPHPrecipitationWorldClimPrecipitationMODISDigital elevation model (DEM)Gridded Population of the WorldHuman population densityEuropean Environment AgencyCorine land cover

Statistical analysis {#Sec5}
--------------------

### Occurrence model {#Sec6}

Three distribution modelling techniques suitable for occurrence data were applied: non-linear discriminant analysis \[[@CR25]\], random classification forest \[[@CR31]\] and generalised linear model \[[@CR32]\]. For each model, the accuracy was assessed using (i) sensitivity, i.e. the ability of a model to correctly identify known positive sites; (ii) specificity, i.e. the ability of a model to correctly identify known negative sites; (iii) the area under the curve, (AUC) that can be roughly interpreted as the probability that a model will correctly predict positive and negative sites \[[@CR33]\]. Of the three techniques, random forest provided the best accuracy and therefore the results of this model are presented.

A random classification forest model consists of an ensemble of trees. To create a reliable model, it is generally considered necessary to have the same number of presence and absence points as input. This is because having a different number will create a bias in the model prediction towards the more prevalent category (presence or absence) \[[@CR33]\]. For this reason, a 'balanced' subset of the data, i.e., a dataset with the same number of presences and absences, was selected. The output produced by the model is an environmental suitability indicator, expressed as a value between 0 (low suitability) and 1 (high suitability). The predictions are visualised in a map with colours ranging from red (high suitability) to blue (low suitability). A list of the most important variables used in the model is given based on the mean decrease in Gini index \[[@CR31],[@CR34]\]. Random forest allows external validation through a bootstrapping procedure: for each tree, a random subset of the full dataset is sampled with replacement. The model validation is carried out for each tree using the points not used from the full dataset. This validation method is referred to as external, because the model is validated using data that are not used to build the tree. The comparison of the observed and predicted results enables us to calculate accuracy statistics, such as sensitivity and specificity. These measures are calculated for each tree and then averaged to give the overall values.

The predictions produced by the random classification forest were also externally validated against 45 observations from other surveys (Table [1](#Tab1){ref-type="table"}) that reported only the presence of *An. plumbeus.* Comparing the observations obtained with the other surveys and the predictions made by the model using National Mosquito Survey data, the error rate was calculated as the proportion of incorrectly predicted pixels to the total number of points used in the validation.

### Abundance model {#Sec7}

The abundance of the species was modelled using a random regression forest model. The abundance data were transformed according to the formula *log*~*10*~*(abundance + 1)* \[[@CR35]\]*.* Because the aim was to identify areas where mosquito peaks could be expected, only the data collected in months in which peaks were observed were selected (June-September). The predicted environmental suitability obtained with the occurrence model described above, was included as one of the predictor variables for modelling the abundance of the species, as it is frequently done in this type of analysis \[[@CR36]-[@CR40]\]. The predicted abundance is interpreted as the expected maximum number of mosquitoes caught in a trap in a certain pixel. The predictions are visualised in a map with colours ranging from light green (low abundance) to dark green (high abundance). The importance of the predictors was assessed using the Increase in Node Purity (INP). The difference between observed and predicted values was expressed as the mean squared error. The analysis has been performed with the software Vecmap demo version \[[@CR41]\]. The maps have been produced with Quantum GIS \[[@CR42]\].

Results and discussion {#Sec8}
======================

The probability of occurrence (environmental suitability) and the abundance of *An. plumbeus* have been predicted using mosquito field data and environmental data. The estimated environmental suitability and abundance are shown in maps. The important environmental variables used in the models and the accuracy of the models are discussed. The fact that out of three different modelling techniques for occurrence data random forest model was selected based upon its higher classification accuracy is consistent with earlier findings; random forest has been reported to outperform other traditional modelling techniques \[[@CR43]-[@CR45]\].

*An. plumbeus* was found in 100 locations out of 766, and it was observed in particular in the eastern part of the Netherlands (Figure [2](#Fig2){ref-type="fig"}a). The percentage of traps containing *An. plumbeus* per week is shown in Figure [3](#Fig3){ref-type="fig"}. For the modelling, 97 presence points and 97 absence points were selected. Using the random forest model, forest-rich areas in the centre of the Netherlands (e.g. National Park Hoge Veluwe and National Park Utrechtse Heuvelrug) are predicted as suitable for *An. plumbeus* (Figure [2](#Fig2){ref-type="fig"}b). Also the eastern parts of the country and in particular the southeastern provinces (Limburg and North Brabant) are predicted to be suitable. In these two provinces nuisance is often reported, especially close to abandoned and un-cleaned pig stables, where mosquitoes breed in manure pools \[[@CR13]\]. Based on the environmental characteristics included in the analysis, the model was capable of identifying areas where *An. plumbeus* is truly present, meaning that these characteristics can be a good proxy for abandoned and un-cleaned stables.Figure 2Observed presence and absence points and map of the estimated environmental suitability for *An. plumbeus*. **A**- Presence and absence observed during the National Mosquito Survey program carried out from April to October 2010--2013. **B**- Environmental suitability map of *An. plumbeus* produced using classification random forest. Environmental suitability is depicted using a gradient fill: blue indicates low environmental suitability, red indicates high suitability. Locations where other surveys took place are also shown on the map (black squares).Figure 3Percentage of positive sites of *An. plumbeus* per week in 2010--2013.

Elevation, precipitation, day and night temperature, vegetation indices and middle infra-red (index sensitive to changes in the vegetation) were found to be important predictors for environmental suitability for *An. plumbeus* (Table [4](#Tab4){ref-type="table"}). Precipitation and vegetation are likely to be biologically relevant, since this species typically breeds in water-filled tree holes with high organic material content \[[@CR19]\]. Presence of tree-holes is related to the age of the tree or to the tree species. Eggs of this mosquito species are not laid on the water surface but on the sides of potential breeding sites, just above the waterline, so the number of generations produced each year are often dependent upon hydrological conditions \[[@CR19]\]. The occurrence model predicts environmental suitability for *An. plumbeus* in areas where old trees with tree-holes are known to occur (e.g., National Park Hoge Veluwe). Even though the species is considered to be a tree-hole breeding species, results obtained using the random forest occurrence model depict environmental suitability for this species in areas without forests in the Netherlands. The results indicate the potential successful use of unforested environments for *An. plumbeus* populations, and imply a similar trend to that seen in continental Europe and UK, where *An. plumbeus* is shifting habitats from almost exclusively breeding in tree-holes to exploiting a wider array of novel man-made larval breeding sites \[[@CR4],[@CR46]\]. Day and night temperature were found to be related to environmental suitability also in another study in Belgium \[[@CR1]\].Table 4**List of the top 10 most important variables in the occurrence modelRankVariables**1DEM2CMORPH precipitation, phase of bi-annual cycle3CMORPH precipitation, phase of annual cycle4Worldclim precipitation, phase of annual cycle5Worldclim precipitation, proportion of total variance due to annual cycle6MIR, phase of annual cycle7NTLS temperature, minimum value8DTLS temperature, amplitude of annual cycle9NDVI mean10CMORPH precipitation, maximum valueThe lowest ranking number indicates the most important variable (e.g., rank = 1 is the most important variable). The ranking is based on the mean decrease in Gini index.

Fair accuracy was obtained with the model (AUC = 0.73), which showed a better ability in identifying suitable environments (sensitivity 0.71) than unsuitable environments (specificity 0.66). The accuracy is improved as compared with a first attempt of predicting the environmental suitability for *An. plumbeus* in the Netherlands, where the environmental suitability was extrapolated from Belgium to the Netherlands (sensitivity = 0.50, specificity = 0.49) \[[@CR1]\]. The error rate calculated to compare the predicted values to data of other surveys was low (0.29); 71% of the pixels were correctly predicted, meaning that the model could make good predictions in non-sampled areas. However, this is only a partial validation because it considers only presence points and does not give information about the performance of the model in predicting absence points.

The observed and estimated abundance are shown in Figure [4](#Fig4){ref-type="fig"}. In the summer, the observed field abundance was low; rarely more than 10 mosquitoes per trap were found (Table [5](#Tab5){ref-type="table"}). The maximum value observed was 1701 mosquitoes followed by 62 mosquitoes per trap. Given the huge difference between the maximum value and the numbers of mosquitoes observed in the other traps, the maximum value was considered as outlier and excluded from the analysis. However, the reason of this high abundance was investigated and it turned out that the trap was located in a rural area where cattle farms with abandoned stables, not in use anymore, are present. In total, 505 mosquitoes were used in the analysis and they were captured in 80 traps/weeks from 2010 to 2013. The predicted abundance, produced with the random forest technique, was also low (with a maximum of 15 individuals per trap) and the highest abundance was predicted in the eastern part of the country and especially in the south-eastern provinces Limburg and Brabant. These findings are in agreement with the suitable areas identified by the occurrence model. This is not surprising because the probability of occurrence was the most important factor among the environmental factors included in the abundance model (Table [6](#Tab6){ref-type="table"}). Similar results, where the probability of occurrence appears to be the most important factor when using this technique, were already observed in a study conducted on other mosquito and biting midges species in the Netherlands \[[@CR35]\]. Other important predictors for abundance of *An. plumbeus* were precipitation and vegetation, in accordance with the biology of the species, as it is explained above. The difference between the observed abundance and the predicted abundance was small, with a mean squared error value of 0.12, meaning that the model predictions matched the observation.Figure 4Observed and estimated abundance of *An. plumbeus*. **A** -- Observed abundance represented as *log* ~*10*~ *(abundance + 1)*. **B** -- Map of the estimated abundance produced using a regression random forest. A darker colour indicates higher abundance.Table 5**Observed abundance used in the modelCountFrequency**0661-107011-20421-30131-40241-50251-60061-701Table 6**List of the top 10 most important variables in the abundance modelRankVariables**1Occurrence2Worldclim precipitation, phase of annual cycle3Worldclim precipitation, proportion of total variance due to bi-annual cycle4NDVI, amplitude of annual cycle5Worldclim precipitation, amplitude of bi-annual cycle6MIR, amplitude of annual cycle7DEM8NTLS temperature, phase of bi-annual cycle9Worldclim precipitation, total variance10CMORPH precipitation, phase of bi-annual cycleThe lowest ranking number indicates the most important variable (e.g., rank = 1 is the most important variable). The ranking is based on the Increase in Node Purity (INP).

Conclusions {#Sec9}
===========

The aim of this study was to investigate the potential spatial distribution of *An. plumbeus* in the Netherlands. Using random (classification and regression) forest models, we identified areas with high environmental suitability and high abundance of this species in south-eastern provinces of Limburg and Brabant. These areas coincide with the areas where in recent years most nuisances have been reported. The predictions of the occurrence model were accurate and matched the external dataset used for validation. The abundance model predictions also matched the observation.

The output of species distribution modelling method can be used as an input for risk assessment of establishment and spread of vector-borne diseases \[[@CR47],[@CR48]\]. Understanding and depicting the potential spatial distribution of mosquito species with modelling techniques is of increasing importance, especially for nuisance mosquito species that can cause economic implications or impact on human health.

Adolfo Ibañez-Justicia and Daniela Cianci contributed equally to this work.

**Competing interests**

The authors declare that they have no competing interests.

**Authors' contributions**

DC and AIJ designed the study. AIJ collected the data. DC did the data preparation and performed the statistical analysis. DC and AIJ wrote the paper. Both authors read and approved the final version of the manuscript.

The authors would like to thank Nienke Hartemink and Marieta Braks for a critical reading of the manuscript and useful comments. The authors would also like to thank Els Ducheyne and Wesley Tack for helpful discussion and suggestions. This study was funded by EU grant FP7-261504 EDENext and is catalogued by the EDENext Steering Committee as EDENext311 ([http://www.edenext.eu](http://www.edenext.eu/)). The contents of this publication are the sole responsibility of the authors and do not necessarily reflect the views of the European Commission.
