Study sites and design
The field surveys were carried out in six small to medium-sized towns (3000 to 16,000 inhabitants) around the border area between Ticino in Switzerland and the Lombardy region in Italy (Fig. 1). The municipalities are located in the historical-geographical region of Insubria. The climate of this region is characterized by dry and sunny winters, with periods of foehn wind from the North with occasional heavy snowfall, rainfall, especially in the transitional seasons (spring and autumn), and sunny summers interrupted by downpours that can also be violent. The landscape of the region features foothills and typical components of Lombard agriculture next to residential, industrial and commercial urbanized areas. The particular geographical position of Insubria has been an incentive to build strong economic relations between Ticino and the neighbouring Italian provinces, resulting in intense traffic across the border, with on average over 67,900 Italian workers commuting to Switzerland in a single workday [16].
The six municipalities surveyed are located within a radius of 7 km, have similar dimension and urban structure, with a small-town centre surrounded by residential areas, and similar climatic (Additional file 1: Figure S1) and altitudinal characteristics (255–414 m a.s.l.). The three municipalities (i.e. Balerna, Coldrerio and Mendrisio) in the Mendrisiotto district in Ticino follow systematically the cantonal IVM since 2009 with monthly or weekly treatments of catch basins with diflubenzuron- or Bti-based products, respectively [5]. The treatment period starts at beginning mid-May, depending on the precipitation pattern, and lasts until mid-end September. In 2019, to encourage citizens to treat permanent breeding sites on private properties, Bti-based granules (VectoBac® G) were made freely available in each municipality. We categorized these three municipalities as “intervention” areas. The three municipalities in the provinces of Como (i.e. Maslianico and Uggiate-Trevano) and Varese (i.e. Malnate), in Lombardy, to our knowledge, did not follow an IVM and only applied adulticides irregularly. They were, therefore, categorized as “non-intervention” areas.
The territory of each municipality was divided into a grid of 250 × 250 m cells [5, 14]. Six cells, called sampling sites, were selected at random in urban context in each municipality. An ovitrap and a Gravid Aedes Trap (GAT, Biogents, Germany) were placed in each sampling site at a distance of 20–100 m from each other to avoid as much as possible interference in mosquito attraction on the ground at hidden, shaded, wind-protected locations close to vegetation. All traps were geo-referenced and uniquely labelled. The ovitraps were the same as used for the surveillance of Ae. albopictus in Ticino [5, 14]. Both ovitraps and GATs mimic breeding sites, attracting container-breeding mosquitoes in search of an oviposition site [13]. Ovitraps allow the invasive mosquitoes to deposit their eggs on a wooden slat and to fly away, while GATs capture the mosquitoes by means of an adhesive plastic sheet [17].
Sample collection and processing
Field surveys were carried out from mid-end May (calendar week 21) to the beginning of October (calendar week 41) 2019. The slats of ovitraps and the sheets of the GATs were replaced at the same time every 14 (range 10 to 19) days with new ones [14]. As a result, ten collection rounds were executed for each ovitrap, except for ovitraps in Uggiate-Trevano, where the survey started 2 weeks later as for the survey with all GATs.
In addition to the already established Ae. albopictus and the indigenous species Aedes geniculatus, two other invasive container-breeding mosquito species, i.e. Aedes japonicus and Aedes koreicus, have started spreading across the study area since 2013 [18]. Although Ae. japonicus and Ae. koreicus have a lower vectorial competence and therefore lower public health significance compared to Ae. albopictus [19], their presence can introduce a confounding factor in the surveillance of Ae. albopictus. Indeed, while eggs of Ae. geniculatus can be easily distinguished by morphology from the other Aedes species, it is not possible to distinguish morphologically the eggs of Ae. albopictus, Ae. japonicus and Ae. koreicus without resorting to special microscopy equipment and expertise [20].
After transportation to the laboratory, the ovitrap wooden slats and GAT adhesive plastic sheets were examined using a stereo microscope (EZ4 D, Leica Microsystems, Germany) for the presence of Aedes eggs and adults, respectively. Adult mosquito females in GATs were identified to the species level by morphology and enumerated. To evaluate the presence of Ae. koreicus and Ae. japonicus in ovitraps, most (87%) of the positive ovitraps were analysed by matrix-assisted laser desorption/ionization time of flight mass spectrometry (MALDI-TOF MS). The rest of the positive ovitraps (13%) could not be analysed because of the low number (1–3) and low quality (dryness or other types of damage) of eggs present on the wooden slat. For this analysis, each wooden slat was divided into ten sectors. For each sector where eggs were present, three to five intact eggs were randomly picked and identified through MALDI-TOF MS with an AXIMA Confidence mass spectrometer (Shimadzu Biotech, Kyoto, Japan) following the method described in Schaffner et al. [20].
Data analysis
The numbers of Ae. albopictus eggs and adult females were recorded in an Excel sheet together with additional information such as the sampling site, date, category of area, etc. (Additional file 2: Table S1, Additional file 3: Table S2 and Additional file 4: Table S3). Statistical analysis was performed through the freely available software R, version 4.0.3 [21]. All analyses are fully reproducible as data, code and package version control tools are available (Additional file 2: Table S1, Additional file 3: Table S2, Additional file 4: Table S3, Additional file 5: Dataset S1, Additional file 6: Dataset S2, Additional file 7: Dataset S3, Additional file 8: Text S1, Additional file 9: Text S2 and Additional file 10: Text S3). All statistical analyses are documented in Additional file 11: Text S4, Additional file 12: Text S5, Additional file 13: Text S6 and Additional file 14: Text S7.
A Spearman’s rank order correlation was used to evaluate the relationship between the number of eggs per ovitrap and the number of Ae. albopictus adult females in the GAT deployed in the same sampling site. Both variables were square-root transformed.
Three modelling analyses were performed. The first two analyses aimed at modelling the number of Ae. albopictus eggs found in ovitraps and the number of Ae. albopictus adult females caught with GATs in 2019, respectively. The third analysis focused on the number of eggs found in ovitraps in 2012, 2013 and 2019. The data for the years 2012 and 2013 come from the study published by Suter et al. [14]. These data are freely available at https://doi.org/10.1371/journal.pntd.0004315.s001. The area monitored in 2012 and 2013 (Additional file 15: Figure S2) was larger compared to the one monitored in 2019 but it included five of the six municipalities surveyed in 2019, namely Balerna, Coldrerio, Mendrisio, Maslianico and Uggiate-Trevano. Additionally, the non-intervention municipality Malnate, surveyed in 2019, was located just outside the southwest boundary of the area studied in 2012 and 2013. Therefore, we believe that the results of both surveys are comparable as well as representative of the situation in the urban communities across the Swiss-Italian border. Note that the study performed by Suter et al. [14] collected data in both sylvatic and urban environments. Since the present study focused on urban environments, for the third analysis we only used the “urban” data from Suter et al. [14] and the data collected in 2019.
The graphical analysis was performed with the ggplot2 package, version 3.3.3 [22]. All three analyses were performed with the glmmTMB function from the glmmTMB package version 1.0.2.1 [23]. Inference was performed with likelihood ratio tests (for P-values) and profiling likelihood methods (to estimate confidence intervals) [24]. The level of significance was set at α = 0.05. Different distributional families and non-nested models were compared with information criteria [25]. Model assumptions were assessed via usual residuals analyses. Quadratic effects were modelled via orthogonal polynomials.
First model
The response variable “number of eggs” (No..eggs.AEDES) was modelled with a generalised mixed-effects model. In particular, to account for the nature of the data, a negative binomial distribution was assumed. This allowed accounting for the fact that we are dealing with count data and that overdispersion is present. Alternative families where compared (Additional file 11: Text S4). The predictor of main interest “AREA” defined whether the trap was to be found in a sampling site under IVM (i.e. in intervention area) or not (i.e. in non-intervention area) and was included as a fixed effect. The other predictors were “Municipality”, “TRAP.ID.fac” (i.e. trap identity), “Day of the year”, “No..Days.ovitrap.in.field” (i.e. number of days that the trap was deployed in the field) and “Altitude” (i.e. altitude of the trap in meters a.s.l.). Municipality and TRAP.ID.fac were taken as random effects. No..Days.ovitrap.in.field was included to account for the “exposure” effect, as not all traps stayed exactly 14 days in the field (range 10 to 19 days). Traps that are left longer in the field are expected to contain more eggs. The seasonal effect of time (i.e. date when ovitrap collected, “Day of the year”) was modelled with a quadratic effect.
The equation used to fit the model is:
$$\begin{gathered} {\text{glmmTMB}}({\text{No}}..{\text{eggs}}.{\text{AEDES }}\sim {\text{ AREA }} + \hfill \\ {\text{poly}}\left( {{\text{Day}}.{\text{ovitrap}}.{\text{collected}},{\text{ degree }} = \, 2} \right) \, + \hfill \\ {\text{scale}}\left( {{\text{ALTITUDE}}} \right) \, + \hfill \\ {\text{No}}..{\text{Days}}.{\text{ovitrap}}.{\text{in}}.{\text{field }} + \hfill \\ \left( {1 \, |{\text{ TRAP}}.{\text{ID}}.{\text{fac}}} \right) \, + \, \left( {1 \, |{\text{ MUNICIPALITY}}} \right), \hfill \\ {\text{family }} = \, {\text{nbinom}}1, \hfill \\ {\text{data }} = {\text{ d}}.{\text{eggs}}.2019) \hfill \\ \end{gathered}$$
Second model
The same analysis as in the first model was applied to the response variable “number of Ae. albopictus adult females in GAT” (No..Ad..Albo.in.GAT hereafter). As the graphical analysis indicated that there might be an interaction between “Day of the year” and “AREA”, we fitted a model that included this two-fold interaction (Additional file 12: Text S5).
Third model
The number of eggs in urban areas in 2012 and 2013 from Suter et al. [14] and the number of eggs collected in 2019 were also modelled together. About half (45%) of the observations in this dataset were zeros. Therefore, the negative binomial used for the two other analyses was extended such that the excess of zeros could be accounted for. To do that, we fitted a hurdle model, where a part of the model focused on the presence-absence part of the data and another part of the model focused on the abundance (Additional file 13: Text S6). The presence-absence part of the model was modelled with binomial family, while the abundance part of the model was modelled with the truncated negative binomial family. The predictors used here were the same as in the previous two analyses, except for the “No..Days.ovitrap.in.field”, not present in Suter et al. [14] data and therefore not included. In addition, the variable “Year” (2012, 2013 and 2019) was added to the model. This model is also a generalised mixed-effects model.