Mapping of Aedes albopictus Abundance at a Local Scale in Italy

Given the growing risk of arbovirus outbreaks in Europe, there is a clear need to better describe the distribution of invasive mosquito species such as Aedes albopictus. Current challenges consist in simulating Ae. albopictus abundance, rather than its presence, and mapping its simulated abundance at a local scale to better assess the transmission risk of mosquito-borne pathogens and optimize mosquito control strategy. During 2014–2015, we sampled adult mosquitoes using 72 BG-Sentinel traps per year in the provinces of Belluno and Trento, Italy. We found that the sum of Ae. albopictus females collected during eight trap nights from June to September was positively related to the mean temperature of the warmest quarter and the percentage of artificial areas in a 250 m buffer around the sampling locations. Maps of Ae. albopictus abundance simulated from the most parsimonious model in the study area showed the largest populations in highly artificial areas with the highest summer temperatures, but with a high uncertainty due to the variability of the trapping collections. Vector abundance maps at a local scale should be promoted to support stakeholders and policy-makers in optimizing vector surveillance and control.


Introduction
The increasing spread of invasive mosquito species (IMS) and the (re-)emergence of mosquito-borne pathogens, such as the dengue, chikungunya, and Zika viruses, are a major threat to public health worldwide [1,2].Among IMS, Aedes (Stegomyia) albopictus (Skuse) is rapidly expanding around the globe because of international travel and trade, and climate change [3,4].In Europe, Ae. albopictus has colonized almost all of the Mediterranean countries, and has been involved in the local transmission of the dengue and chikungunya viruses in Italy, France, and Croatia [5].Given the growing risk of arbovirus outbreaks in Europe, it is crucial to better describe the distribution of IMS such as Ae.albopictus based on species distribution models and mapping.Correlative species distribution models (SDMs) relate species distribution data at known locations with information on the environmental and/or spatial characteristics of these locations [6]; these models differ from the mechanistic ones that mathematically describe the biological processes underpinning population performance [7].In addition, Remote Sensing (RS) and Geographic Information Systems (GIS) are helpful tools for predicting and mapping species distribution [8,9].
To evaluate the transmission risk of mosquito-borne pathogens, one main challenge of SDMs is to simulate the abundance rather than the presence of mosquito species; mosquito abundance data related to host-seeking females are more reliable than occurrence data (i.e., the presence or presence/absence of data), since mosquito abundance is directly related to the intensity of pathogen transmission [7,10,11].However, most of the previous studies have proposed SDMs for Ae.albopictus based on occurrence data, i.e., presence only or presence-absence data, and have produced suitability maps at a global scale that do not fit well with the scale at which surveillance and control measures may be implemented [7,12].On the contrary, maps at a local scale should be promoted to stakeholders and policy makers [7,13].The identification of areas highly infested by Ae. albopictus, hereby referred to as hot spots, is crucial in terms of pathogen transmission risk and mosquito control strategy.Indeed, mosquito control interventions in hot-spot areas optimize integrated mosquito management by increasing their cost-effectiveness [14].
The distribution of IMS depends on multiple factors, including suitable climate (e.g., for adult survivorship and the overwintering survival of eggs) and land use-land cover, hereafter LULC (e.g., for the availability of larval habitats and resting sites).Previous suitability maps of Ae. albopictus have often been built based on SDMs including climate variables as predictors, mainly temperatures (e.g., the mean temperature of the coldest month, the mean temperature of the warmest quarter), but also precipitations (e.g., annual precipitation, precipitation in spring or summer) [12,15].By contrast, LULC has been only marginally considered as a predictor of Ae. albopictus distribution, even though Ae. albopictus population is strongly associated with urbanized areas [4,10,[16][17][18][19][20][21][22][23].
In the present study, we aimed to model the abundance of Ae. albopictus females collected from adult traps, including both meteorological and LULC variables, and produce a map of the simulated abundance of this species in a specific area of north-eastern Italy.We hypothesized that LULC variables, particularly the proportion of artificial areas, could be an important driver of Ae.Albopictus abundance.This study was based on original sampling data collected at a local scale, which is a good prerequisite to characterize the distribution of a species using SDMs where well-designed survey data and functionally relevant explanatory variables are analyzed with appropriate models [6].Thus, we could use a correlative SDM which us enabled to assess the effect of land use on Ae. albopictus abundance, and we converted the sample points into surface data in order to produce a continuous-surface map.

Study Area and Mosquito Data
The study area was located in the provinces of Belluno and Trento in north-eastern Italy.This is a mountainous area influenced by a continental climate, with more than 70% of the territory lying over 1000 m above sea level, and 55% covered by coniferous and deciduous forests [11].It includes well developed agricultural, industrial, and commercial areas with main transport axes [24].The human population as of 2014 is around 208,000 in the Belluno Province and 537,000 in the Trento Province (http://www.istat.it),with most of the population concentrated in the valley floors.Aedes albopictus is well established in north-eastern Italy since its first detection in 1991.Populations of Ae. albopictus were found for the first time in the provinces of Trento and Belluno in 2001 and 2005, respectively (Simone Martini, personal communication, 7 June 2016) [25].Previous monitoring activities in the study area have shown that Ae. albopictus had colonized the main valleys and areas up to 800 m above sea level.[26].
A total of 72 sampling locations were selected in areas well-colonized by Ae. albopictus, ranging from 76 to 663 m above sea level in 2014 and 2015 (68 locations were selected in 2014, and 4 out of the 68 were changed in 2015) (Figure 1).We chose these locations in order to sample different municipalities in the provinces of Belluno and Trento along an altitudinal gradient and with different land use patterns.We ensured that no insecticide treatments were applied around the sampling locations to avoid biases due to mosquito control measures.Adult mosquitoes were collected using BG-sentinel traps (Biogents AG, Regensburg, Germany).Adult trapping is a reliable method to estimate a mosquito population's size [27,28], and the BG-sentinel trap is the most effective trapping device to collect Ae. albopictus adults, as observed in many studies conducted worldwide [29][30][31][32].The BG-sentinel traps were positioned by skilled entomologists to optimize mosquito collections.They were mainly located in artificial areas with surrounding vegetation as recommended by the supplier, such as private properties, public buildings, cemeteries, and farms.The minimum distance between two sampling locations was around 200 m to avoid any trap competition.Each BG-sentinel trap was baited with a BG Lure (Biogents AG, Regensburg, Germany) and CO 2 from a thermos bottle filled with 750 g of dry ice pellets.The traps were operated for 24 h fortnightly during the whole mosquito breeding season, i.e., from 18 April to 4 November 2014 and from 21 April to 3 November 2015.Practically, we sampled half of the locations one week and half of the locations the week after, so that each location was sampled every two weeks.The trapping sessions were carried out considering weather conditions (i.e., sunny days with low wind speed), since heavy rainfall or strong winds can disturb mosquito collection.Each time, the traps were set between 9.00 a.m. and 12.00 p.m. and collected the day after between 9.00 a.m. and 12.00 p.m. following the same setting order.
Remote Sens. 2017, 9, 749 3 of 14 municipalities in the provinces of Belluno and Trento along an altitudinal gradient and with different land use patterns.We ensured that no insecticide treatments were applied around the sampling locations to avoid biases due to mosquito control measures.Adult mosquitoes were collected using BG-sentinel traps (Biogents AG, Regensburg, Germany).Adult trapping is a reliable method to estimate a mosquito population's size [27,28], and the BG-sentinel trap is the most effective trapping device to collect Ae. albopictus adults, as observed in many studies conducted worldwide [29][30][31][32].The BG-sentinel traps were positioned by skilled entomologists to optimize mosquito collections.They were mainly located in artificial areas with surrounding vegetation as recommended by the supplier, such as private properties, public buildings, cemeteries, and farms.The minimum distance between two sampling locations was around 200 m to avoid any trap competition.Each BG-sentinel trap was baited with a BG Lure (Biogents AG, Regensburg, Germany) and CO2 from a thermos bottle filled with 750 g of dry ice pellets.The traps were operated for 24 h fortnightly during the whole mosquito breeding season, i.e., from 18 April to 4 November 2014 and from 21 April to 3 November 2015.Practically, we sampled half of the locations one week and half of the locations the week after, so that each location was sampled every two weeks.The trapping sessions were carried out considering weather conditions (i.e., sunny days with low wind speed), since heavy rainfall or strong winds can disturb mosquito collection.Each time, the traps were set between 9.00 a.m. and 12.00 p.m. and collected the day after between 9.00 a.m. and 12.00 p.m. following the same setting order.The adult mosquitoes were killed by freezing at −20 °C and then identified using taxonomic keys [26,33].Over the 2-year survey, the beginning and the end of the sampling could vary among locations, resulting in differences in the total number of trap nights.Therefore, we estimated Ae. albopictus abundance per sampling location per year with the sum of females only collected during the eight trap nights between June and September, corresponding to the maximum activity period for mosquitoes in the study area.Five samples with less than eight trap nights during this period (i.e., missing data because of field problems) were not included in the statistical analysis.The sum of females collected in the summer can be considered an estimator of the Ae.albopictus infestation level in a sampling location, and the number of Ae. albopictus females trapped by BG-sentinel traps may be related to the number of females biting humans [34].

Environmental Data
To model the abundance of Ae. albopictus females, we selected several meteorological and LULC variables according to the literature [10,12,15].
We obtained precipitation data from meteorological stations located close to the sampling locations (Table 1).The total annual precipitation, precipitation of the second quarter (i.e., April-May-June), and precipitation of the third quarter (i.e., July-August-September) for 2014 and 2015 were calculated from the dataset freely available at the ARPA Veneto (Agenzia Regionale per la Prevenzione e Protezione Ambientale del Veneto) and Meteotrentino websites The adult mosquitoes were killed by freezing at −20 • C and then identified using taxonomic keys [26,33].Over the 2-year survey, the beginning and the end of the sampling could vary among locations, resulting in differences in the total number of trap nights.Therefore, we estimated Ae. albopictus abundance per sampling location per year with the sum of females only collected during the eight trap nights between June and September, corresponding to the maximum activity period for mosquitoes in the study area.Five samples with less than eight trap nights during this period (i.e., missing data because of field problems) were not included in the statistical analysis.The sum of females collected in the summer can be considered an estimator of the Ae.albopictus infestation level in a sampling location, and the number of Ae. albopictus females trapped by BG-sentinel traps may be related to the number of females biting humans [34].

Environmental Data
To model the abundance of Ae. albopictus females, we selected several meteorological and LULC variables according to the literature [10,12,15].We obtained precipitation data from meteorological stations located close to the sampling locations (Table 1).The total annual precipitation, precipitation of the second quarter (i.e., April-May-June), and precipitation of the third quarter (i.e., July-August-September) for 2014 and 2015 were calculated from the dataset freely available at the ARPA Veneto (Agenzia Regionale per la Prevenzione e Protezione Ambientale del Veneto) and Meteotrentino websites (http://www.arpa.veneto.it/;http://www.meteotrentino.it/).Rainfall is essential to provide water-filled containers for egg deposition and larval development; an annual precipitation of less than 500 mm is considered a limiting factor for Ae.albopictus [35].We used two sets of Land Surface Temperature (LST) data: a 2-year dataset (2014-2015) based on the years of sampling to carry out the abundance model using mosquito collections from 2014-2015, and a wider dataset based on a 13-year long-term average (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) to map the simulated abundance of Ae. albopictus females in the study area (because the long-term averages better characterized the temperatures of the region than the 2-year averages, which only considered 2014 and 2015).The temperatures for the period 2003-2012 gathered from the EuroLST bioclim dataset (spatial resolution to 250 m) are freely available at the EuroLST website (http://www.geodati.fmach.it/eurolst.html).The temperatures for 2013, 2014, and 2015 were obtained from the MODIS version 5 LST products MOD11A1 and MYD11A1 with four records per day and a spatial resolution of 1000 m (Table 1); we applied outlier filtering and gap filling procedures to these data, and enhanced the spatial resolution to 250 m [36].LST data were extracted from the 250 m pixels corresponding to the trap locations.We computed the mean temperature of the coldest month (i.e., January) and the mean temperature of the warmest quarter (i.e., July-August-September) in 2014, 2015, and for the period 2003-2015.According to the literature, an average winter temperature higher than 0 • C enables the overwintering of diapausing eggs of Ae. albopictus, and a mean summer temperature ranging from 25 to 30 • C is optimal for the establishment of its populations [15,37].
The LULC data were obtained from the dataset freely available at the open data Trentino and Veneto websites (http://dati.trentino.it/;http://dati.veneto.it/)(Table 1).The LULC datasets from Trentino and Veneto were elaborated in 2003 and 2013, respectively, at a scale of 1:10000.Buffers of 250 m were overlaid on the sampling locations using QGIS [38]; the size of the buffer was based on the maximum daily dispersal of Ae. albopictus, ranging from 200 to 300 m [39,40].We intersected the 250 m buffer with the LULC data to extract the surface areas of the different LULC classes included in each buffer, and we aggregated the data to the first-level class.For each buffer, we calculated the percentages of artificial areas (urban fabrics, industrial, commercial, and transport units . . .), agricultural areas (arable lands, permanent crops, and pastures . . .), and forests and semi-natural areas (forests, scrubs, and open spaces with little vegetation . . . ) (See Figure S1 in the supplementary materials).Previous studies have shown that urban areas, open agricultural areas, and forested areas affected Ae. albopictus presence [10,16,23,41].

Statistical Analysis
We investigated the association between Ae. albopictus female abundance (i.e., the sum of females collected during eight trap nights between June and September) and a range of meteorological and LULC variables based on the collection carried out during 2014-2015 using R version 3.2.2[42].Collinearity among environmental variables was tested using Pearson's correlation test (See Figure S2 in supplementary materials).We excluded variables with an absolute value of Pearson's coefficient >0.7, and the selected variables were standardized by subtracting their mean values and dividing by their standard deviations.Afterwards, we performed negative binomial mixed models using the R package 'glmmADMB' [43].The explanatory variables included in the full model were the total precipitation of the second quarter, the mean temperature of the coldest month, the mean temperature of the warmest quarter, and the percentages of artificial areas, forests, and semi-natural areas.The trap location was included as a random effect to take into account the variability among sampling locations [44].We carried out multi-model inference to compare all possible models using the R package 'MuMIn' [45], as described in Rosà et al. [46].The models were ranked using AIC, and differences in AIC (∆AIC) between consecutively ranked models were used to calculate the weights and relative evidence ratios for each variable.The best models were selected using a threshold of ∆AIC ≤2, and the most parsimonious model was used to map the simulated abundance of Ae. albopictus females.
Model assumptions, for the most parsimonious selected model, were evaluated by checking the model's normalized residuals for any pattern or dependency following the protocol as in Zuur et al. [47].The model's performance was assessed using in-sample errors by computing the root mean squared error (RMSE), which represents the sample standard deviation of the differences between predicted and observed values, and could be interpreted as an estimation of the standard deviation of the unexplained variance.In addition, during the model validation process, a simulation study was carried out to assess if the simulated distribution of values generated by the model were in compliance with the observed values.
To map the simulated abundance of Ae. albopictus females across the study area, we built a 250 m cell raster grid as a reference spatial layer in GRASS GIS 7 [48].The grid was built to both best overlap the study area and match the EuroLST grid.The mean temperature of the warmest quarter was sampled for each cell of the grid.Then, we extracted the percentage of artificial areas for each grid cell from the LULC dataset.The percentage of artificial areas in each cell was defined as the number of pixels belonging to the level class 1 (according to the Corine Land Cover classification) over the total number of pixels in each grid cell.We then simulated the abundance of Ae. albopictus females in each grid cell using the selected model as described above.Simulations were performed using the R function 'predict.glmmadmb',which only accounts for the fixed effects of the model.Simulations were done using the link scale (logarithmic scale) in order to compute the uncertainty (i.e., the standard error of the simulated mosquito abundance) of our simulations.Afterwards, the simulated values were back-transformed (using the exponential function) to obtain the simulations on the scale of the response variable, while the standard error was transformed in the 95% Confidence Interval (CI) of the simulated abundance value through the formula: CI 95% = UL − LL with UL = exp(fit + 1.96 * se.fit) and LL = exp(fit − 1.96 * se.fit).The simulated values and their 95% CI were finally reported as two distinct maps: a map of the study area showing the simulated Ae. albopictus abundance per pixel, and a map reporting the uncertainty associated with the simulated abundance in each pixel.The most parsimonious model selected for simulating Ae. albopictus female abundance included the mean temperature of the warmest quarter and the percentage of artificial areas in a 250 m buffer around the sampling locations as significant explanatory factors (Table 2).The model indicated that the abundance of Ae. albopictus females would increase by 28.7% and 2.8% for an increase in temperature of 1 °C and an increase in artificial areas of 1%, respectively.Despite a strong association between Ae. albopictus abundance and the two fixed effects, we observed a relatively high standard deviation (SD) of the random effect for each 'trap' (SD = 1.974), leading to a marked uncertainty in the simulated value of Ae. albopictus abundance.

A
The statistical model's assumptions were satisfied and no spatiotemporal correlation was detected in the residuals (see Figure S3 in supplementary materials).The total RMSE between the observed and fitted values was equal to 452.The results of the simulation study, performed during The most parsimonious model selected for simulating Ae. albopictus female abundance included the mean temperature of the warmest quarter and the percentage of artificial areas in a 250 m buffer around the sampling locations as significant explanatory factors (Table 2).The model indicated that the abundance of Ae. albopictus females would increase by 28.7% and 2.8% for an increase in temperature of 1 • C and an increase in artificial areas of 1%, respectively.Despite a strong association between Ae. albopictus abundance and the two fixed effects, we observed a relatively high standard deviation (SD) of the random effect for each 'trap' (SD = 1.974), leading to a marked uncertainty in the simulated value of Ae. albopictus abundance.
The statistical model's assumptions were satisfied and no spatiotemporal correlation was detected in the residuals (see Figure S3 in supplementary materials).The total RMSE between the observed and fitted values was equal to 452.The results of the simulation study, performed during the model validation process, indicate that the 93% of the observed values of mosquito abundance are within the 95% CI of the negative binomial simulated distribution obtained from 10,000 model simulations using covariates measured at each observation (see Figure S4 in supplementary materials).
In Figure 3   The coloured lines represent the simulated abundance (applying the most parsimonious model) of Ae. albopictus females collected in eight trap nights between June and September with respect to the proportion of artificial areas, while varying the temperature of the warmest quarter across its range in the study area (the curves take into account the two explanatory fixed effects, but not the trap random effect included in the negative binomial mixed models).In the background, we reported (as points) the total abundance of adult females observed in each sampling location as a function of the artificial area in a 250 m buffer.The y axis was limited at 300 adult females for graphical reasons.
Figure 4 shows the map of the simulated abundance of Ae. albopictus females in the provinces of Belluno and Trento using the mean temperature of the warmest quarter and the percentage of artificial areas in 250 m pixels.The colours ranged from purple-blue for areas with low abundances of Ae. albopictus to orange-red berry for areas with high abundances.The simulated abundance of Ae. albopictus was mainly related to artificial areas.We simulated the largest populations of Ae. albopictus in municipalities located along the two main valleys, Val d'Adige (where Trento and Rovereto are the two biggest towns) and Val Belluna (where Belluno and Feltre are the two biggest towns), and the area close to the Lake of Garda.A high uncertainty of Ae. albopictus female abundance was observed especially in areas with a high percentage of artificial areas (Figure 5).

Figure 3.
The coloured lines represent the simulated abundance (applying the most parsimonious model) of Ae. albopictus females collected in eight trap nights between June and September with respect to the proportion of artificial areas, while varying the temperature of the warmest quarter across its range in the study area (the curves take into account the two explanatory fixed effects, but not the trap random effect included in the negative binomial mixed models).In the background, we reported (as points) the total abundance of adult females observed in each sampling location as a function of the artificial area in a 250 m buffer.The y axis was limited at 300 adult females for graphical reasons.
Figure 4 shows the map of the simulated abundance of Ae. albopictus females in the provinces of Belluno and Trento using the mean temperature of the warmest quarter and the percentage of artificial areas in 250 m pixels.The colours ranged from purple-blue for areas with low abundances of Ae. albopictus to orange-red berry for areas with high abundances.The simulated abundance of Ae. albopictus was mainly related to artificial areas.We simulated the largest populations of Ae. albopictus in municipalities located along the two main valleys, Val d'Adige (where Trento and Rovereto are the two biggest towns) and Val Belluna (where Belluno and Feltre are the two biggest towns), and the area close to the Lake of Garda.A high uncertainty of Ae. albopictus female abundance was observed especially in areas with a high percentage of artificial areas (Figure 5).

Discussion
In the present study, we mapped the simulated abundance of Ae. albopictus females in an area of north-eastern Italy using the mean temperature of the warmest quarter and the percentage of artificial areas in a 250 m pixel.To the best of our knowledge, this is one of the first efforts to produce an Ae.albopictus abundance map based on original sampling data collected from adult traps, whereas previous studies mostly described suitability maps based on occurrence data.By including LULC variables in the model selection, we improved the accuracy of Ae. albopictus distribution in the study area with respect to previous suitability maps of Ae. albopictus in Europe mainly based on climatic predictors only, especially temperatures [12,15,35,49,50] that did not allow the studies to differentiate artificial, agricultural, or forest areas for Ae.albopictus habitat suitability.In our study, we simulated high abundances of Ae. albopictus in highly artificial areas, i.e., urbanized areas.This can be explained by the presence of numerous man-made containers, which are the main sources of Ae. albopictus populations in urban or peri-urban environments [51].Urbanization has been related to an increase in the larval habitats, mosquito density, larval development rate, and adult survivorship of Ae. albopictus [52].Furthermore, several studies conducted worldwide on land use effect demonstrated a strong affinity of Ae. albopictus for urbanized and anthropized environments [10,16,19,22,23].
On the other hand, our model analysis did not show any significant relationship between the percentage of forest and semi-natural areas and Ae.albopictus abundance.In our survey, traps were mainly located in areas surrounded by vegetation.These small patches of vegetation represented suitable resting sites for Ae.albopictus; however, at the scale of our analysis (i.e., 1:10,000), these patches could not be differentiated from artificial areas.In addition, half of the sampling locations were characterized by a low percentage of forest and semi-natural areas (median = 10.5%), and this could explain why no effect of forest areas was detected.By contrast, in previous studies, vegetation areas have been positively associated with Ae. albopictus abundance, especially inside urban environments.In Hawaii, Vanwambeke et al. [18] predicted the highest abundance of Ae. albopictus in places with a mix of built-up areas and vegetation, while in Italy, Manica et al. [53] described small green areas as hot spots for Ae.albopictus in a metropolitan environment.
Despite the well-characterized Ae. albopictus distribution in the study area, we found that a high uncertainty of the simulated abundance of females was due to the variance of the 'trap' effect.This finding reflects the high variability in mosquito collections observed among traps set in a nearby area (e.g., within a 1 km radius).The high variability in the collected Ae. albopictus may be explained by several characteristics of the sampling locations, including the availability of breeding sites for oviposition and larval development, the availability of resting sites for mosquito adults, the co-occurrence of competing mosquito species in larval habitats, and the presence of vertebrate hosts as blood meal sources [53].Among these characteristics, it can be assumed that the spatial heterogeneity of breeding sites underpins the high variability in mosquito local abundance.However, the relationship between breeding sites' availability and LULC is difficult to estimate at a fine-scale, limiting the accuracy of the abundance map [54].Therefore, reliable LULC data linked to breeding site availability might improve SDMs and mosquito abundance maps considering the heterogeneity of artificial areas [28,54].
Our model showed that the temperature of the warmest quarter was positively associated with Ae. albopictus abundance, in agreement with previous studies [15,50,55].Temperatures strongly affect mosquito population growth, and thus the density of Ae. albopictus adults [37,56,57].In the range of the mean temperature of the warmest quarter observed in the study area, i.e., 17-25 • C, Ae. albopictus populations were expected to increase with increasing temperatures.Indeed, a mean temperature of the warmest quarter equal to 25 • C represents an optimum temperature for Ae.albopictus [58].On the other hand, the mean temperature of the coldest month was not related to Ae. albopictus abundance, whereas it has been demonstrated that it might play a main role on Ae. albopictus occurrence.In fact, all of the sampling locations were suitable for Ae.albopictus, with mean temperatures of the coldest month above 0 • C, allowing for the survival of diapausing eggs during winter conditions [59].Likewise, the precipitation of the second quarter was not associated with Ae. albopictus abundance, probably because of the low variability of this parameter among the sampling locations.In our study area, the average total annual precipitation is above 500 mm, and so it is not a limiting factor for Ae.albopictus [35].On the contrary, in Southern Europe, especially Spain, low summer precipitation seems to be the most limiting factor for Ae.albopictus [15].However, human water supply in urban areas may provide many breeding sites, so that these areas may become suitable for Ae.albopictus despite low precipitation.Thus, it appeared that several meteorological variables (e.g., the temperature of the coldest month), which could have a quantitative effect on Ae. albopictus presence, did not affect its abundance as suggested by Roche et al. [22].
To model Ae.albopictus abundance, we used adult collections from BG-sentinel traps.These traps are designed to collect host-seeking females, and so they provide a direct observation of the abundance of Ae. albopictus females that can be related to the human biting rate [32,34].BG-sentinel traps have been proved to be very effective for Ae.albopictus; however, adult collections may be affected by trap performance and positioning as well as weather conditions.Furthermore, these traps are costly and require a source of power.Other methods, such as ovitrap collections or a pupal survey, have been recommended to estimate IMS abundance [27,28].In the Emilia-Romagna region of Italy, Carrieri et al. [51,60] demonstrated a positive correlation between the number of Ae. albopictus eggs collected in ovitraps, female density estimated by a pupal survey, and human landing collections (HLC); and in Rome, Italy, Manica et al. [61] estimated an increase of one biting female per person collected by HLC for every five additional eggs observed in ovitraps.Ovitraps are a low-cost and sensitive method to estimate Ae. albopictus abundance.However, egg collections may be influenced by the position of the trap, the proximity of breeding sites, and micro-environmental changes such as the cutting of the vegetation around the trap [62]; in some cases (e.g., Rimini, Italy), ovitrap data may not correctly estimate Ae. albopictus abundance simulated by HLC [62].Therefore, monitoring activities should be carried out or at least supervised by skilled entomologists to minimize the sampling error.Good sampling data is an important prerequisite to model mosquito species distribution.

Conclusions
A projection of climatic suitability for Ae.albopictus indicated an increase in suitable areas in central western Europe within the first half of the 21st century [12,35].In northern Italy, increasing temperatures due to climate change might favour the expansion of Ae. albopictus northward and upward, especially in a mountainous environment [50].Ten years ago, the altitudinal limit for Ae.albopictus was considered to be around 600 m above sea level in Italy, while nowadays this species is observed up to 900 m in the province of Belluno (Fabrizio Montarsi, personal communication, 8 November 2016).Increasing summer temperatures predicted by climate scenarios [63] and the growing urbanization observed during the last 50 years [64] might increase Ae. albopictus populations, and consequently the transmission risk of pathogens such as the chikungunya and dengue viruses [11].Furthermore, northern Italy is also highly suitable for two other IMS, Ae. japonicus and Ae.koreicus, whose invasion success may be facilitated by a spatiotemporal niche segregation with Ae. albopictus [24,31,58].Further research is needed to lower the high uncertainty of prediction in artificial areas and to produce more accurate IMS abundance maps that are required by stakeholders for vector surveillance and control.

Supplementary Materials:
The following are available online at www.mdpi.com/2072-4292/9/7/749/s1, Figure S1: Histograms of the proportions of artificial, agricultural and forest areas in 250 m buffers around the sampling locations (N = 72), Figure S2: Scatterplot of the temperature of the coldest month, the temperature of the warmest quarter, the total annual precipitation, the precipitation of the second quarter, the precipitation of the third quarter and the proportions of artificial, agricultural and forest areas.The upper panels contain the Pearson correlation coefficient with the font proportional to its estimated value.The lower panels show scatterplots with a smoother added to visualize the pattern, Figure S3: Assessment of model assumptions for the most parsimonious model selected.The term 'Site' corresponded to different localities [Belluno (BL), Lamen (LA), Riva del Garda-Arco (RA), Sospirolo (SO), Tisoi (TI), Tenno-Pranzo (TP)].In each locality, 4 to 6 sampling location were selected, Figure S4: Results of the simulation study.Observed values (filled circles) vs. mean values (open circles) and 95% CI (vertical bars) of the negative binomial simulated distribution obtained from 10,000 model simulations using covariates measured at each observation point, Table S1: Original data by sampling location including Aedes albopictus female sum, mean, standard error and range with altitude, weather and LULC variables.

Figure 1 .
Figure 1.Trapping locations with BG-sentinel traps baited with BG lures and CO2.The different colours indicate the year(s) of sampling (2014 and/or 2015) for each sampling location (Microsoft Bing Aerial map).

Figure 1 .
Figure 1.Trapping locations with BG-sentinel traps baited with BG lures and CO 2 .The different colours indicate the year(s) of sampling (2014 and/or 2015) for each sampling location (Microsoft Bing Aerial map).
total of 17,695 Ae. albopictus, including 66.8% females, were collected in 2014 (N = 5188) and 2015 (N = 12,507), during a total of 928 and 846 trap nights in 2014 and 2015, respectively.Aedes albopictus was the most abundant species, representing 75% of the total number of collected mosquitoes.The mean number of specimens collected per trap night are plotted for each year in Figure 2. The sum of Ae. albopictus females collected per trap during a sampling season ranged from 0 to 345 in 2014 and 0 to 870 in 2015.For modelling, we included 131 out of 136 samples representing a total of 9206 females collected during 1048 trap nights (seeTable S1 in supplementary materials).Remote Sens. 2017, 9, 749 6 of 14 0 to 345 in 2014 and 0 to 870 in 2015.For modelling, we included 131 out of 136 samples representing a total of 9206 females collected during 1048 trap nights (see Table S1 in supplementary materials).

Figure 2 .
Figure 2. Mean number of Aedes albopictus females collected per trap night in 2014 and 2015 with BGsentinel traps baited with BG lures and CO2.

Figure 2 .
Figure 2. Mean number of Aedes albopictus females collected per trap night in 2014 and 2015 with BG-sentinel traps baited with BG lures and CO 2 .
is shown the simulated abundance of Ae. albopictus across the range of possible artificial areas for different values of the mean temperature of the warmest quarter across the range estimated in the study area, i.e., from 17 to 25 • C. At the highest temperatures, the simulated female abundances increase exponentially with increasing artificial areas.Remote Sens. 2017, 9, 749 7 of 14 the model validation process, indicate that the 93% of the observed values of mosquito abundance are within the 95% CI of the negative binomial simulated distribution obtained from 10,000 model simulations using covariates measured at each observation (see Figure S4 in supplementary materials).In Figure 3 is shown the simulated abundance of Ae. albopictus across the range of possible artificial areas for different values of the mean temperature of the warmest quarter across the range estimated in the study area, i.e., from 17 to 25° C. At the highest temperatures, the simulated female abundances increase exponentially with increasing artificial areas.

Figure 3 .
Figure3.The coloured lines represent the simulated abundance (applying the most parsimonious model) of Ae. albopictus females collected in eight trap nights between June and September with respect to the proportion of artificial areas, while varying the temperature of the warmest quarter across its range in the study area (the curves take into account the two explanatory fixed effects, but not the trap random effect included in the negative binomial mixed models).In the background, we reported (as points) the total abundance of adult females observed in each sampling location as a function of the artificial area in a 250 m buffer.The y axis was limited at 300 adult females for graphical reasons.

Figure 4 .
Figure 4. Map (datum: ETRS89 European Terrestrial Reference System; projection: LAEA Lambert Azimuthal Equal Area) of the simulated abundance of Aedes albopictus females (sum of eight trap nights from June to September) in the provinces of Belluno and Trento (Italy) based on the most parsimonious model including the mean temperature of the warmest quarter and the percentage of artificial areas in 250 m pixels as predictors.The colour ranges from purple, corresponding to 0 specimens per trap per year, to red-berry, corresponding to 112 specimens per trap per year.

Figure 5 .
Figure 5. Map (datum: ETRS89 European Terrestrial Reference System; projection: LAEA Lambert Azimuthal Equal Area) of the simulation 95% Confidence Interval (CI) of the abundance of Aedes albopictus females (sum of eight trap nights from June to September) in the provinces of Belluno and Trento (Italy).

Figure 4 .
Figure 4. Map (datum: ETRS89 European Terrestrial Reference System; projection: LAEA Lambert Azimuthal Equal Area) of the simulated abundance of Aedes albopictus females (sum of eight trap nights from June to September) in the provinces of Belluno and Trento (Italy) based on the most parsimonious model including the mean temperature of the warmest quarter and the percentage of artificial areas in 250 m pixels as predictors.The colour ranges from purple, corresponding to 0 specimens per trap per year, to red-berry, corresponding to 112 specimens per trap per year.

Figure 4 .
Figure 4. Map (datum: ETRS89 European Terrestrial Reference System; projection: LAEA Lambert Azimuthal Equal Area) of the simulated abundance of Aedes albopictus females (sum of eight trap nights from June to September) in the provinces of Belluno and Trento (Italy) based on the most parsimonious model including the mean temperature of the warmest quarter and the percentage of artificial areas in 250 m pixels as predictors.The colour ranges from purple, corresponding to 0 specimens per trap per year, to red-berry, corresponding to 112 specimens per trap per year.

Figure 5 .
Figure 5. Map (datum: ETRS89 European Terrestrial Reference System; projection: LAEA Lambert Azimuthal Equal Area) of the simulation 95% Confidence Interval (CI) of the abundance of Aedes albopictus females (sum of eight trap nights from June to September) in the provinces of Belluno and Trento (Italy).

Figure 5 .
Figure 5. Map (datum: ETRS89 European Terrestrial Reference System; projection: LAEA Lambert Azimuthal Equal Area) of the simulation 95% Confidence Interval (CI) of the abundance of Aedes albopictus females (sum of eight trap nights from June to September) in the provinces of Belluno and Trento (Italy).

Table 1 .
Mean (standard deviation) and range values of the meteorological and land use variables used for abundance model selection (period 2014-2015).

Table 2 .
Estimated coefficients and statistics of the most parsimonious negative binomial mixed model selected for simulating the seasonal abundance of Aedes albopictus females (i.e., sum of females collected in eight trap nights between June and September).

Table 2 .
Estimated coefficients and statistics of the most parsimonious negative binomial mixed model selected for simulating the seasonal abundance of Aedes albopictus females (i.e., sum of females collected in eight trap nights between June and September).