Mapping the abundance of endemic mosquito-borne diseases vectors in southern Quebec
BMC Public Health volume 23, Article number: 924 (2023)
Climate change is increasing the dispersion of mosquitoes and the spread of viruses of which some mosquitoes are the main vectors. In Quebec, the surveillance and management of endemic mosquito-borne diseases, such as West Nile virus or Eastern equine encephalitis, could be improved by mapping the areas of risk supporting vector populations. However, there is currently no active tool tailored to Quebec that can predict mosquito population abundances, and we propose, with this work, to help fill this gap.
Four species of mosquitos were studied in this project for the period from 2003 to 2016 for the southern part of the province of Quebec: Aedes vexans (VEX), Coquillettidia perturbans (CQP), Culex pipiens-restuans group (CPR) and Ochlerotatus stimulans group (SMG) species. We used a negative binomial regression approach, including a spatial component, to model the abundances of each species or species group as a function of meteorological and land-cover variables. We tested several sets of variables combination, regional and local scale variables for landcover and different lag period for the day of capture for weather variables, to finally select one best model for each species.
Models selected showed the importance of the spatial component, independently of the environmental variables, at the larger spatial scale. In these models, the most important land-cover predictors that favored CQP and VEX were ‘forest’, and ‘agriculture’ (for VEX only). Land-cover ‘urban’ had negative impact on SMG and CQP. The weather conditions on the trapping day and previous weather conditions summarized over 30 or 90 days were preferred over a shorter period of seven days, suggesting current and long-term previous weather conditions effects on mosquito abundance.
The strength of the spatial component highlights the difficulties in modelling the abundance of mosquito species and the model selection shows the importance of selecting the right environmental predictors, especially when choosing the temporal and spatial scale of these variables. Climate and landscape variables were important for each species or species group, suggesting it is possible to consider their use in predicting long-term spatial variationsin the abundance of mosquitoes potentially harmful to public health in southern Quebec.
In Canada, the arbovirus most commonly transmitted to humans by mosquitoes are West Nile virus (WNV), California serogroup (CSG) viruses, and, less frequently, Eastern equine encephalitis virus (EEEV) [1,2,3,4,5,6,7]. Quebec is one of the provinces that is most concerned by the transmission of these viruses. Since 2002, the number of human cases of WNV has varied considerably from year to year, with an average of 27 cases by year and some particularly extreme year like 2018 with 201 confirmed cases, including 15 deaths . No human cases of infection by EEEV have been observed so far, but some positive horses and emus have been observed sporadically, mainly in the Laurentides region, north of Montréal . The number of human cases caused by the CSG viruses is hard to compare between years because cases with any kind of symptoms have been notifiable only since 2019 (before that date, only encephalitis cases were notifiable). Following enhanced surveillance in 2017, 82 cases were detected in Quebec, whereas in 2018 and 2019, without any enhanced surveillance, only 17 and 18 cases respectively were reported . In Quebec, CSG viruses, EEEV and WNV are on a watch list for health professionals, and are all notifiable diseases [8,9,10]. Better characterization of the risk associated with mosquitoes likely to transmit diseases to humans is therefore a public health issue for Canada in general, and Quebec in particular.
One tool helping public health assess the situation and for decision-making is the threat exposure risk map, or, in this case, the spatial distribution of the mosquito population . It allows for optimization of control (larviciding) or prevention (message to the public) programs. Modelling mosquito abundances is one way to predictively map exposure risks [12,13,14,15,16,17]. The variables employed for such a mapping often involve the use of meteorological data, such as temperature and precipitation [12, 14, 15, 18,19,20,21,22,23,24,25] and land-cover data, because the populations being modelled are arthropods, organisms that are highly dependent on their environment and the prevailing weather conditions. The environmental data used as explanatory variables vary widely across studies and cover the following major themes: land cover and land use, vegetation indices, hydrographic indices, drainage indices, wetlands, slope and altitudes [18, 20, 24,25,26,27,28].
In Quebec, Culex pipiens-restuans (CPR) is the main vector species for WNV . It is also suspected to have some vector competence for EEEV [30, 31]. CPR preferences in terms of habitat have been studied in several regions comparable with our study region. According to studies conducted in Illinois and New York state, CPR is found in woodlands in particular and, to a lesser extent, in fields, grasslands and urban areas, laying their eggs in artificial and natural containers [20, 26, 32, 33]. More recent studies performed in Ontario and southern Quebec highlighted the importance of urban dwellings, proximity of suburban backyards and edge vegetation and mixed or paved areas for the bioecology of CPR [17, 34, 35]. Similarly, favorable climatic conditions for CPR have been discussed in several studies. The climate variables used vary from one study to another, but in general, CPR abundance is associated with lagged variables, representing a mean or cumulative effect over several days prior to the day of capture, of temperature and precipitation variables [20, 23, 36, 37]. As an example, in southern Quebec, CPR abundance was associated with the average over 22 days of the degree days (Celsius degree over 9 °C only, developmental threshold for CPR) before mosquito capture and mean precipitation over 71 days before capture .
Aedes vexans (VEX), which is very abundant in southern Quebec, has some vector competence for WNV [38,39,40], EEEV , SSHV  and JCV [39, 41]. It spawns primarily in the temporary or more permanent aquatic areas that form in agricultural and wetlands areas, in proximity to urbanized landscapes [17, 34, 42, 43]. Other sources indicate that VEX prefer to develop in pools, shallow pasture ponds, woodlands and grasslands . As with CPR, we observe that VEX abundance is associated with both temperature and precipitation lag variables [22, 36, 37].
Coquillettidia perturbans (CQP) larvae require aquatic plants, such as reeds or water lilies; they get oxygen from them by piercing their roots. It is therefore a species that requires permanent bogs or marshes for the larval stage of its life cycle [42, 45]. However, while CQP is very common in wetlands and ‘forest’s [46, 47], not much more is known about the other environmental conditions necessary for their development, other than the importance of precipitation . CQP has some known vector competence for EEEV  primarily and for WNV  and JCV [41, 48] anecdotally.
The environmental preferences of the Oc. stimulans group (SMG) are poorly documented. It seems that SMG prefers to spawn in vernal pools, especially in woodlands [17, 42, 49]. The group is identified as a member of the “spring aedes/ochlerotatus” group and is therefore very abundant in spring . Ochlerotatus stimulans (Oc. stimulans), included in the Oc. stimulans group (SMG), has known vector competence for JCV  and SSHV . The SMG also include some rarer species such as Oc. fitchii and Oc. riparius. Nothing is reported in literature about their vector competence for JCV, SSHV, WNV or EEEV.
Not all of these mosquito species have the same preferences in terms of larval habitat and climate sensitivity. It is therefore essential that the models developed are species-specific and combine spatial and temporal aspects [12, 14, 25, 51].
In the present work, we study the variations of the spatiotemporal daily abundance of CPR, SMG, CQP and VEX mosquitoes as a function of land-cover and weather variables in southern Quebec by developing of series of spatio-temporal regression models. The purpose of these models was to map the daily mean mosquito abundance, for each of the four species, to guide campaigns to prevent and control the spread of the diseases borne by the selected species.
The study area is located in southern Quebec (Fig. 1) and follows the human ecumene of Greater Montréal. The study area is crossed from southwest to northeast by the St. Lawrence River. Most of the areas south of the river are dense agricultural areas with some open woodlands and the odd wetland. North of the river, past a narrow agricultural strip, there is a ‘forest’ area populated by mixed and coniferous forest whose density increases toward the north. According to the Köppen-Geiger classification , most of our study area has a warm-summer humid continental climate (Dfb boreal climate), with the warmest month of the year (July) averaging below 22 °C, but more than four months averaging above 10 °C. The temperature in the coldest four months varies between − 38 °C and 0 °C. Precipitation is distributed evenly throughout the year. In these types of climates, the season of mosquito activity generally extends from May to September.
Mosquito abundance data
We selected four species or groups of species that pose a potential public health risk because of their respective vector competence for arboviruses: CPR, CQP, VEX and SMG. Culex pipiens and Culex restuans were lumped into a single group Cx. pipiens-restuans (CPR) due to their morphological similarity. Oc. stimulans, Oc. fitchii and Oc.riparius are also grouped together for this reason (SMG group). In this last group, because of the geographic distribution of the three species and their relative abundance, we can reasonably consider that the main or most abundant species in the SMG group is Oc. stimulans . The mosquito catch data came partly from the provincial program under the Institut national de santé publique du Québec (INSPQ) and Quebec’s Ministère de la Santé et des Services sociaux (MSSS) for the surveillance of mosquito-borne diseases in Quebec, and partly from the Public Health Agency of Canada (PHAC) using the database of a private company involved in mosquito control (GDG Environnement Ltd.). Trapping was conducted with Center for Disease Control (CDC) CO2 light traps and mosquitoes were collected at each site at night. The dataset covers the years 2003 to 2016, but the number of traps deployed from year to year varied and fewer traps were installed in some years and in more geographically restricted areas (Fig. 2). The number of traps deployed also varied from week to week (Fig. 3). Trapped mosquitoes were identified and counted by species.
Variables used for modelling mosquito abundances included weather and land-cover class variables (see below for a description). Seasonality in mosquito abundance patterns was taken into account by using day-of-year as an explanatory variable in all models as it integrates a lot of information relevant to mosquito development and abundance cycles (see the Statistical Modelling section below for details of how it was entered in the models).
Weather variables were obtained from the Daymet estimates of daily surface weather data on a 1-km grid for North America  using the daymetr package . Prior to extraction, we aggregated the 1 × 1 km cells to 2 × 2 km cells to reduce computation time and to fill in the few cells with missing data. Since only the daily maximum and minimum temperature are available from the Daymet data, we estimated the daily mean temperature by calculating the mean of the daily maximum and minimum temperatures. For each trap catch, we extracted the daily mean temperature and the daily precipitation (prcp). Because of the strong seasonal trend in the daily mean temperature and the high correlation between day-of-year and temperature, we removed the seasonal trend in the mean temperature to consider temperature anomalies instead (anom). Specifically, we calculated for each pixel the difference between the daily mean temperature and the mean daily mean temperature for the period from 2003 to 2016. We calculated the pluriannual daily mean temperature across the study period using a generalized additive model (GAM) with the daily mean temperature as our response variable and a smooth function with a cyclic basis on the day-of-year as our predictive variable. Temperature anomalies were obtained by subtracting the daily mean temperatures from the predicted daily mean temperatures. The R package mgcv  was used to build the GAM models. We did not consider precipitation anomalies as the seasonal trend in precipitations was much less pronounced. Finally, to take into account the effect of previous weather conditions on mosquito abundance, we considered lagged effects of daily temperature anomalies and precipitations by calculating the mean of the variables over the previous 2, 7, 30 and 90 days counting from the date of each catch. We thus had the current (two days) and the previous (7, 30 and 90 days) weather conditions as predictors of mosquito abundance. The weather conditions over two days correspond to the weather conditions on the day the trap is set and the day it is lifted.
The first category of environmental data used is the land-cover class (LCC) layers derived from previously classified satellite (raster) data and vector data (Table 1). Eight LCCs were selected initially for our exercise based on our literature review: aquatic area (shallow and deep water), bare soil, urban/builtup, wetland, agricultural land, pasture/grassland, vegetated non-treed, and woodland. The raster layer processed by Agriculture and Agri-Food Canada (AAFC)  was used as a general base, and its classification was improved by the addition of multi-source data: the classified SPOT 4–5 raster land cover of Canada south of the treeline provided by the Government of Canada (GoC)  and two other layers in vector format were used for this purpose, namely the Ducks Unlimited Canada (DUC)  wetland map and Natural Resources Canada’s National Hydrographic Network (NHN) . The different sources of LCC data were combined using the following decision rule: in non-urban areas, the information provided by the AAFC layer prevailed. In ‘urban’ areas, the information provided by the SPOT image prevailed over the AAFC image. The LCC was classified as a wetland if and only if the data provided by AAFC, SPOT and DUC all indicated the presence of a wetland. The LCC was classified as an aquatic area based on the NHN information, which prevailed over all other data sources.
Although we had access to several land cover classes, trap locations were not necessarily chosen to cover the range of possible values for each class within the study area and most traps were installed near inhabited areas. Hence, we were not able to use some likely important classes for explaining mosquito abundance, such as wetland cover, which would have led to extrapolation beyond the range of observed values. We instead opted for using the three most abundant and variable classes in the study area, namely ‘urban’, ‘forest’ and ‘agriculture’ cover (see Fig. 1) for which we had good coverage and thus avoided extrapolation. We considered percentage of land cover around each trap at both a local (50 m buffer) and a regional (1 km buffer) scale. Although the environment immediately surrounding each trap likely has an important impact on the number of captures, the 30 m resolution of the land cover data prevented us from considering more local effects. Since the land-cover data used were not available for the years from 2003 to 2010, we used the land cover from the year 2011 and applied it to all traps from 2003 to 2016 on the assumption that changes in land cover were relatively small.
We modeled for each species/group the number of mosquitos captured for each trap night as a count using a negative binomial likelihood with a log link. We included a spatial component to account for unexplained variability in mosquito abundance across space and to take into account spatial autocorrelation. Spatial models were built using INLA  and the SPDE approach . Briefly, INLA is an alternative to MCMC for estimating latent Gaussian models in a Bayesian framework. The SPDE approach facilitates the inclusion of spatial effects by approximating Gaussian random fields through a link with Gaussian Markov Random Field (GMRF). The method proved very useful for estimating large spatial and spatiotemporal models  and for disease mapping [65,66,67].
Penalized complexity (PC) priors were used to reduce complexity of the spatial component whenever the data did not warrant it [68, 69]. The PC priors used for the range parameter controlling the extent of the spatial dependence in mosquito counts represented a probability of 0.1 that the range is smaller than 5 km. The PC prior for the standard deviation of the spatial fields represented a probability of 0.1 that the standard deviation of the spatial fields is over 0.5. Since this standard deviation has to be exponentiated to go from the link to the response scale, it represents substantial variations in mosquito counts. Basically, unless warranted by the data, the PC priors used shrink the parameters toward a model where the spatial effect is small and operating at a large scale. The shape parameter of the Matérn covariance function was set to the default value of two. The maximum side length of the mesh triangles used to approximate the spatial component of the model in the SPDE approach was set to 0.5 km. The details of mesh construction can be found in the R code available upon request.
All our models had day-of-year and a mix of weather and land-cover variables as predictors (see below for details of the model selection process). Natural splines were used to model the effect of day-of-year and to account for non-linearities in mosquito abundances across the season (function ns in R package splines). To facilitate the process of modelling and generating predictions with INLA, we specified the location of knots used for the splines. We subjectively chose a number of knots (10) that seemed to show a good compromise between capturing the main pattern of variation across the season while avoiding unlikely complexities. We used relatively tight normal priors on spline coefficients to facilitate model convergence (mean = 0, sd = 1). Prior to model formulation, multicollinearity between predictors was assessed using pairwise correlations and variance inflation factors (VIF) with the car R package . We avoided using the three land-cover classes as VIF values were over three when all classes were used in the same models. Thus, for each species, we chose the two land-cover classes most likely to have an influence on abundance according to the literature and to the biology of the species, naming ‘agriculture’ and ‘forest’ for VEX, and ‘urban’ and ‘forest’ for all the other species or group of species. Predictors included in a given model either had a correlation lower than 0.7 with any other predictor and had VIF values below three. Due to high correlation between mean temperature at longer time lags (i.e. 30, 90 days) and day-of-year, we avoided including both day-of-year and mean temperatures in models and instead used anomaly values which lacked a seasonal trend.
All variables were scaled to have a mean of 0 and a standard deviation of 1 prior to analysis. Wide normal priors were used for all variables (mean = 0, sd = 30). An offset was used to account for the number of nights a trap was deployed, although the vast majority of traps were deployed for a single night each time. Models were estimated with the R package INLA (v. 22.01.25) using the experimental mode  and in R 4.1.2 (R Core Team 2021). The following packages were use for spatial tasks and figure constructions: sp, sf, terra, exactextractr magick, rnaturalearth. The code is available on request.
Model selection was used to determine at which spatial or temporal scale weather and land-cover variables operate. For each species, we built a series of models that tested whether local and/or regional land cover and current and/or previous weather were better predictors of mosquito abundance. We built 3 sets of landcover variables (‘urban50 + forest50’, ‘urban1000 + forest1000’ and ‘urban50 + forest50 + urban1000 + forest1000’) where 50 refers to local scale and 1000 to regional scale. We also built 7 sets of weather variables (‘anom2 + prcp2’, ‘anom7 + prcp7’, ‘anom30 + prcp30’, ‘anom90 + prcp90’, ‘anom2 + prcp2 + anom7 + prcp7’, ‘anom2 + prcp2 + anom30 + prcp30’,’ anom2 + prcp2 + anom90 + prcp90’) where 2, 7, 30 and 90 refers to the previous days counting from the date of each catch. The combination of these 2 sets of variables brings us to 3 × 7 = 21 models for each species.
We compared each model within a set using the deviance information criterion (DIC), choosing the best model whenever a model clearly stood out or the most complete model when several models had similar DIC scores. Although the Watanabe-Akaike information criterion (WAIC) is considered superior to the DIC as it uses the entire posterior distribution, it is not recommended for models where there is dependence between observations . We therefore used the DIC, keeping in mind that it is one of many methods and that model selection remains an open question . To evaluate the usefulness of the spatial component, we also included in the set of models a model with the same explanatory variables as in the best model, but with the spatial component removed. To determine whether the importance of the spatial component was partly driven by the many explanatory variables we could not consider within a single model, we also considered a model with all explanatory variables, again with the spatial component removed.
Model checking and validation
Model fit was assessed through simulations by comparing data simulated from the models with the observed data. Histograms of simulated and observed count values were visually compared and showed an acceptable fit. The DHARMa package  was also used to assess zero-inflation and model fit through randomized quantile residuals . Only the best model for each species was used to assess model fit. For all models, observed values were comparable and did not differ substantially from simulated values. Although many trap counts were zero, there was no indication of substantial zero-inflation with our data and preliminary analyses showed zero-inflation probabilities very close to zero. It is our interpretation that the negative binomial distribution used was able to cope with the number of zeros along with the smooth function used on day-of-year and the spatial component. We used the correlation between observed and predicted counts as a pseudo-R square measure and as a rough estimate of the explanatory power of the models. We compared the best spatial model with its non-spatial version to estimate the importance of the spatial component in explaining variations in abundance. We also compared both models with the non-spatial model containing all explanatory variables to determine whether the spatial component importance was mainly due to missing explanatory variables that we had considered for this study.
More than 3 million mosquitos were captured in the ~ 14,000 trap nights (Fig. 4). The maximum count for a single trap night for each species/group varied from 3,104 (CPR) to 33,024 (CQP).
For most species, models containing both local and regional land-cover variables and both current and previous weather variables were selected as the best model, except for CQP, where the best model only contained the regional land-cover variables (Table 2). The weather conditions during the trapping day and previous weather conditions summarized over 30 or 90 days were preferred over a shorter period of seven days, suggesting current and long-term previous weather conditions effects on mosquito abundance. In all cases, the best model with the spatial component removed performed worse than any other model including the spatial component. This is also reflected in the pseudo-R square, where the best spatial model always had a pseudo-R square at least three times higher than the corresponding non-spatial model (Fig. 5). Models with the spatial component have interesting pseudo-R square values, especially for CQP, which is having a pseudo-R square of more than 45%. In general, non-spatial models with the best explanatory variables had relatively low pseudo-R square values, suggesting that variations in mosquito abundance were marginally explained by our explanatory variables. Except for CQP, non-spatial models with all explanatory variables had pseudo-R square values much lower than the best spatial model, suggesting that the importance of the spatial component was not entirely due to the limited amount of variables we could consider in each model.
Some of the weather variables retained in the best model were averaged for either the 30 (CPR, VEX) or 90 (CQP, SMG) previous days, indicating that previous weather can have a cumulative effect on mosquito abundance (Table 2). However, most effects were relatively weak, except for VEX, where higher precipitations had a strong positive effect on counts (Fig. 6). Current weather also seemed to have some effects on abundance, especially for VEX abundance, which is positively associated with high temperature during the day of capture. High precipitation on the day of capture seems to have a weak effect on diminishing the abundance for CPR and VEX. Precipitation over 90 days before capture has a strong positive effect on SMG abundance.
In general, land-cover variables had a relatively weak effect on abundance. In addition, land-cover variables measured at a regional scale (1 km) had a stronger effect on abundance when compared with the land cover measured on a local scale (50 m). For CPR, abundance decreases slightly with the presence of forest at the local scale. For CQP, abundance strongly increased with the percentage of forests and strongly decreased with the percentage of ‘urban’ land cover within 1 km. SMG abundances decreased with ‘urban’ land cover and the regional scale but increase with the same landcover class, ‘urban’, at the local scale. Finally, VEX abundances increased with ‘forest’ and ‘agriculture’ land cover at the regional scale.
For all species, seasonality (day-of-year) seemed to have the largest effect on abundance. For CQP and SMG, the peak in abundance appears more concentrated compared with CPR and VEX, for which the peak in abundance appears more spread out across the season (Fig. 7). Contrary to other species, SMG abundances are characterized by an earlier and shorter season compared with other species.
All mapped predictions were characterized by very wide credible intervals, which translates to very high uncertainties in predictions, especially when looking at predictions far from where the traps were installed (Fig. 8). Indeed, the standard deviation is high and more or less constant across a large part of the study area. Coupled with the relatively small range of the spatial field, predictions outside of areas where traps were deployed are characterized by large uncertainties. This is also easily seen in Fig. 7, where the spatial component accounts for the larger part of the uncertainty around predictions and, in most cases, the uncertainty associated with the spatial component appears more important than the variation associated with the explanatory variables. Nevertheless, maps show that there are very large variations in abundance across the study area associated with land-cover variables and the spatial component. For all species, abundances appear much lower than expected in the eastern and central part of Montréal, while abundances are higher in various areas in the western part and around the Island of Montréal. Figure 9 and Additional files 1, 2, 3 and 4 shows both the intra-seasonal and spatial variations, for the year 2003, of mosquito abundance, and it illustrates how the abundance of mosquitoes can be variable during the season and in space, depending on the landscape.
This paper led to the creation of maps of mean daily abundance for four mosquito species or species groups, namely CPR (Culex pipiens-restuans group), CQP (Coquillettidia perturbans), SMG (Ochlerotatus stimulans-fitchii-riparius group) and VEX (Aedes vexans), that are key to public health in southern Quebec. This type of map is of immediate interest to support decision-making for the management, control and prevention of mosquito-borne diseases. However, our modelling work also highlights challenges in modelling and mapping mosquito abundances. Despite incorporating environmental predictors at several spatial and temporal scales, the results obtained show that the variance explained by these variables remains low and that the unexplained variance remains predominant, as indicated by the strength of the spatial component and the high uncertainty in predicted abundance across the study area.
Globally, in the best models, we show the need for both current (for weather) or local (for land-cover), which describe the conditions of the day of capture and at the capture site, and the variables with a regional effect (for the land-cover variables) or cumulative effect (for weather variables), which describe what happened in the general habitat surrounding the trap or during the weeks prior to capture. Considering environmental variables measured at different temporal and spatial scales is already well established, although not always applied, and is a key aspect of mosquito abundance modelling . We can easily explain the role of these different scales, in particular the temperature conditions or the precipitation over several days before the capture, which is likely linked to availabilities of breeding sites and to the speed of development of the pre-adult stages, and the temperature and precipitation during the day of capture, which influence the daily activity of adult mosquitoes.
VEX is the mosquito species for which the environmental predictors appeared most important. They include the positive role of high precipitation over the month preceding capture as well as the temperature on the day of capture. This is coherent with what is described in other studies realized in Quebec  and Ontario . Furthermore, for VEX, increased abundance in agricultural landscapes is also well documented [34, 43]. The positive effect of ‘forest’ is less frequently found in the literature, although Hopkins et al. described VEX’s presence in forest areas subject to human disturbance , which corresponds to the type of forests we have in our study region.
The most important predictors for CPR had a relatively weak effect but underline the positive effect of higher-than-normal temperature during a period of one month prior to the day of capture, as well as the negative effect of precipitation on the day of capture. This is consistent with previous work realized in Quebec and Ontario [12, 22, 23, 36]. Negative impacts of ‘forest’ cover at the local scale can be compared with results observed in other studies that showed how an urbanized landscape increases the abundances of this species [17, 32, 33]. Indeed ‘forest’ could be considered as opposite of ‘urban’ in term of landscape. But again, this shows the difficulties of correctly describing these types of relationships with various spatial scales and data sources. We will come back to this point later.
For CQP, the environmental predictors with the largest effect are ‘forest’ land cover, at the regional scale, which tends to increase the mosquito abundance, and ‘urban’ land cover, again at the regional scale, which diminishes its abundance. These observations are hard to directly reconcile with the literature, as CQP is described as a mosquito species dependent on wetland, mainly because its larva requires permanent bogs and marshes to complete its life cycle [42, 45]. Unfortunately, wetland was not a land cover we were able to include in our analyses, because of the trapping design, as explained in the methods section. However, assuming that urban is negatively correlated to wetland, our results appear very coherent with the mosquito biology.
Finally, the most important predictors for the SMG group are the ‘urban’ land cover, at the regional and local scale, which respectively decreases and increases abundance, and the temperature anomalies over 90 days before capture, which negatively impact abundance. The SMG group consists of three species that are designated as “spring” species, abundant early in the season, just after the snow melts, and using forested vernal pools for breeding [17, 42, 49]. These biological characteristics are partially coherent with our best model if we hypothesize that forested vernal pools are less frequent in urban areas, at the regional scale. The positive correlation with ‘urban’ land cover at the local scale is more difficult to explain, although it could be associated with the existence in the mosquito environment of artificial breeding habitat that the mosquito is able to use.
Despite the coherence of our results with mosquito biology and other studies, the larger part of the variance was unexplained by our variables. This suggests that the weather and land-cover variables we used were not sufficient and we may have missed the most important variables for explaining variations in mosquito abundance, or simply have not characterized these variables at the right spatial or temporal scale.
By restricting the number of models and variables used, we may have failed to identify the most important combination of timescale at which weather conditions operate and the scale at which land-cover variables are the most important for mosquito habitat. Although we considered scales in this study, it was not our main objective and we were limited by our model selection approach. More targeted or exploratory approaches such as random forests  may prove useful in studying a wider range of scales, although identifying at which range remains a challenge [79, 80].
Since traps were not specifically deployed for the present study, the range of land-cover classes that could be studied without relying on extrapolation prevented the use of possibly more relevant land-cover classes that are important for mosquito habitats. For example, previous analyses using wetland cover as an explanatory variable led to extreme abundance values predicted in areas with important wetlands, although no traps were installed in such areas.
The large variance and the small range of the spatial field suggest that important variables for explaining mosquito abundance may operate on a very short scale, likely the local habitat variables surrounding traps . However, we did not have access to such detailed habitat information, which likely prevented us from explaining a large part of the variations in mosquito abundances.
The number of traps varied greatly from year to year and trap locations changed across years. We initially sought to use a spatiotemporal model with autocorrelated weekly spatial components, but this led to convergence problems and very large uncertainties. Moreover, the variability in trap location across seasons and years coupled with the small range of the spatial component made predictions vary wildly across the study period. We expect a certain level of variation of abundances across years and seasons within the study area, however, it was unclear whether this variation was genuine or just an artefact of having trap locations changing across years and seasons. We thus opted for a simpler spatial model with which we could more easily focus on the main abundance pattern across the study area. Using a single seasonal effect common to each year facilitated the study of the mean seasonal pattern of mosquito abundance. Our recommendations for conducting mosquito abundance assessments are as follows:
Place the traps so they cover the range of possible values of the presumed most important habitat characteristics to avoid extrapolating mosquito abundance using predictor values not seen by the model.
Identify the most important local-scale habitat variables affecting mosquito abundance, which should facilitate the study of regional mosquito abundance patterns. Otherwise, strong local effects coupled with few trap locations may make the study of regional patterns of mosquito abundance more difficult. This will also likely help in reducing the large uncertainties associated with abundance predictions.
In situations where explained variances seem low, consider using a spatial component to better represent uncertainties in predicted abundances. Although adding a spatial component to a mosquito abundance model may make the prediction of abundance across space difficult, we feel it more adequately represents the uncertainty associated with predictions compared with a model without a spatial component. The importance of the spatial component in this study shows that despite relatively narrow uncertainties around predictions when considering only predictor variables, the uncertainty gets much larger when considering the spatial component, which indicates that abundance may vary wildly across space around what is predicted simply by predictors. Ignoring this uncertainty may give a false sense of confidence in the mapped pattern of mosquito abundances.
The method presented in this paper to map abundances for various mosquito species in southern Quebec is a first step toward developing a practical tool that can be used in decision-making to identify areas of risk of exposure to infected mosquitoes. Indeed, no method had yet been developed in this study area that could be applied on a large scale and to four distinct species. The simplicity of the method used and the use of relatively generic land-cover and climate variables are encouraging and could easily be applied elsewhere. This type of model can also help to study the impact of climate and environmental change on the distribution of these mosquito species in time and space from a public health perspective. Indeed, temperature, precipitation and the presence of ‘urban’ areas are very important variables in our models and are undergoing profound changes: growing urban consumption of land, in tandem with an increase in mean temperatures and a change in precipitation patterns . These three major effects can therefore be expected to significantly alter the landscape of mosquito abundances in the coming years and, consequently, the risk of exposure to mosquito-borne diseases.
The data that support the findings of this study are available from GDG Environnement via the Public Health Agency of Canada, the Institut National de la santé publique du Québec and the Ministère de la santé et des services Sociaux du Québec but restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the corresponding author upon reasonable request and with permission of GDG Environnement via the Public Health Agency of Canada, the Institut National de la santé publique du Québec and the Ministère de la santé et des services Sociaux du Québec.
Agriculture and Agri-Food Canada
Akaike Information Criterion
Center for Disease Control
Culex pipiens-restuans group
Cache Valley virus
Ducks Unlimited Canada
Eastern Equine Encephalitis virus
Government of Canada
Institut national de santé publique du Québec
Jamestown Canyon virus
Ministère de la Santé et des Services sociaux
National Hydrographic Network
Public Health Agency of Canada
St. Louis encephalitis virus
Ochlerotatus stimulans group
Snowshoe hare virus
Variance inflation factor
West Nile virus
Canada PH. A.o. West Nile virus and other mosquito-borne disease national surveillance report 2016—final- summary, P.H.A.o. Canada, Editor. 2017, Government of Canada: Ottawa, Canada. p. 13.
Canada PHA. o. West Nile virus and other mosquito-borne disease national surveillance report 2015, P.H.A.o. Canada, Editor. 2017, Government of Canada: Ottawa, Canada. p. 14.
Canada PH. A.o. West Nile virus and other mosquito-borne diseases in Canada, annual national surveillance report—2017, P.H.A.o. Canada, Editor. 2018, Government of Canada: Ottawa, Canada. p. 4.
Canada PH. A.o. Mosquito-borne disease surveillance report—September 15 to 28, 2019 (week 40 & 41), P.H.A.o. Canada, Editor. 2019, Government of Canada: Ottawa, Canada. p. 4.
Canada PHA. o. West Nile virus and other mosquito-borne diseases surveillance report: annual edition, 2018, P.H.A.o. Canada, Editor. 2020, Government of Canada: Ottawa, Canada. p. 11.
Canada PH. A.o. Mosquito-borne disease surveillance report—September 27 to October 24, 2020 (week 40 to 43), P.H.A.o. Canada, Editor. 2020, Government of Canada: Ottawa, Canada. p. 5.
Canada PH. A.o. Mosquito-borne diseases surveillance report - September 27 to October 24, 2021 (week 39 to 42), P.H.A.o. Canada, Editor. 2021, Government of Canada. p. 5.
Québec. I.N.d.s.p.d. Virus du Nil occidental – tableau des cas humains - archives 2002–2020. Institut National de santé publique du Québec: Montréal, Québec, Canada; 2021.
Québec IN. d.s.p.d. Surveillance des maladies d’intérêt transmises par les moustiques au Québec - Encéphalite équine de l’Est. 2020 [cited 2021 30 November 2021]; Available from: https://www.msss.gouv.qc.ca/professionnels/zoonoses/surveillance-des-maladies-d-interet-transmises-par-des-moustiques-au-quebec/encephalite-equine-de-l-est/.
Québec IN. d.s.p.d. Surveillance des maladies d’intérêt transmises par les moustiques au Québec - Les virus du sérogroupe Californie. 2021 [cited 2021 30 November 2021]; Available from: https://www.msss.gouv.qc.ca/professionnels/zoonoses/surveillance-des-maladies-d-interet-transmises-par-des-moustiques-au-quebec/les-virus-du-serogroupe-californie/.
Kotchi SO, et al. Using Earth observation images to inform risk assessment and mapping of climate change-related infectious diseases. Can Commun Dis Rep. 2019;45(5):133–42.
Lebl K, Brugger K, Rubel F, Predicting. Culex pipiens/restuans population dynamics by interval lagged weather data. Parasites & Vectors, 2013. 6(1): p. 129.
Rochlin I, et al. Predictive mapping of human risk for West Nile virus (WNV) based on environmental and socioeconomic factors. PLoS ONE. 2011;6(8):e23280.
Yoo EH. Site-specific prediction of West Nile virus mosquito abundance in Greater Toronto Area using generalized linear mixed models. Int J Geogr Inf Sci. 2014;28(2):296–313.
Roiz D, et al. Climatic effects on mosquito abundance in Mediterranean wetlands. Parasites & Vectors. 2014;7(1):333.
Sun H, et al. Spatio-temporal analysis of the main dengue vector populations in Singapore. Parasites & Vectors. 2021;14(1):41.
Giordano BV, Turner KW, Hunter FF. Geospatial analysis and seasonal distribution of West Nile virus vectors (Diptera: Culicidae) in Southern Ontario, Canada. Int J Environ Res Public Health, 2018. 15(4).
Cleckner HL, Allen TR, Bellows AS. Remote sensing and modeling of mosquito abundance and habitats in coastal Virginia, USA Remote Sensing, 2011. 3(12).
El Adlouni S, et al. Effects of climate on West Nile virus transmission risk used for public health decision-making in Quebec. Int J Health Geogr. 2007;6:40–0.
Jacob BJ, et al. Developing operational algorithms using linear and non-linear squares estimation in Python for the identification of Culex pipiens and Culex restuans in a mosquito abatement district (Cook County, Illinois, USA). Geospat Health. 2009;3(2):157–76.
Wang X, et al. Clustering of the abundance of West Nile virus vector mosquitoes in Peel Region, Ontario, Canada. Environ Ecol Stat. 2014;21:651–66.
Ripoche M, et al. Short-term forecasting of daily abundance of West Nile virus vectors Culex pipiens-restuans (Diptera: Culicidae) and Aedes vexans based on weather conditions in southern Québec (Canada). J Med Entomol. 2019;56(3):859–72.
Wang J, Ogden N, Zhu H. The impact of weather conditions on Culex pipiens and Culex restuans (Diptera: Culicidae) abundance: a case study in Peel Region. J Med Entomol. 2011;48(2):468–75.
Lim A-Y, et al. Mosquito abundance in relation to extremely high temperatures in urban and rural areas of Incheon Metropolitan City, South Korea from 2015 to 2020: an observational study. Parasites & Vectors. 2021;14(1):559.
Fornasiero D, et al. Inter-annual variability of the effects of intrinsic and extrinsic drivers affecting West Nile virus vector Culex pipiens population dynamics in northeastern Italy. Parasites & Vectors. 2020;13(1):271.
Trawinski PR, Mackay DS. Identification of environmental covariates of West Nile virus vector mosquito population abundance. Vector Borne Zoonotic Dis. 2010;10(5):515–26.
Ha TV, et al. Spatial distribution of Culex mosquito abundance and associated risk factors in Hanoi, Vietnam. PLoS Negl Trop Dis. 2021;15(6):e0009497.
Chen L, Zhu H, Wang X. Modeling spatiotemporal distribution of mosquitoes abundance with unobservable environmental factors. J Med Entomol. 2019;56(1):65–71.
Reisen WK. Ecology of West Nile virus in North America. Viruses. 2013;5(9):2079–105.
Andreadis TG, Anderson JF, Tirrell-Peck SJ. Multiple isolations of eastern equine encephalitis and highlands J viruses from mosquitoes (Diptera: Culicidae) during a 1996 epizootic in southeastern Connecticut. J Med Entomol. 1998;35(3):296–302.
Oliver J, et al. Twenty years of surveillance for eastern equine encephalitis virus in mosquitoes in New York State from 1993 to 2012. Parasites & Vectors. 2018;11(1):362.
Gardner AM, Lampman RL, Muturi EJ. Land use patterns and the risk of West Nile virus transmission in central Illinois. Vector-Borne and Zoonotic Diseases. 2014;14(5):338–45.
Holmes CJ, Cáceres CE. Predation differentially structures immature mosquito populations in stormwater ponds. Ecol Entomol. 2020;45(1):97–108.
Cloutier CA, Fyles JW, Buddle CM. Diversity and community structure of mosquitoes (Diptera: Culicidae) in suburban, field, and forest habitats in Montréal, Québec, Canada. The Canadian Entomologist, 2021. 153(4): p. 393–411.
Moua Y et al. Mapping the habitat suitability of West Nile virus vectors in southern Quebec and eastern Ontario, Canada, with species distribution modeling and satellite earth observation data. Remote Sens, 2021. 13(9).
Dussault C, et al. Evaluating the impact of Aedes japonicus invasion on the mosquito community in the Greater Golden Horseshoe region (Ontario, Canada). PLoS ONE. 2018;13(12):e0208911.
Chuang TW, et al. Cross-correlation map analyses show weather variation influences on mosquito abundance patterns in Saginaw County, Michigan, 1989–2005. J Med Entomol. 2012;49(4):851–8.
Turell MJ, et al. Potential north american vectors of West Nile virus. Ann N Y Acad Sci. 2001;951(1):317–24.
Anderson JF et al. Arboviruses in North Dakota, 2003–2006. The American Society of Tropical Medicine and Hygiene, 2015. 92(2): p. 377–393.
Goddard LB, et al. Vector competence of California mosquitoes for West Nile virus. Emerg Infect Disease J. 2002;8(12):1385.
Main AJ et al. Arbovirus surveillance in Connecticut. II. California serogroup [Aedes species, insect vectors]. 1979. v. 39.
Maire A, Aubin A. Les moustiques du Québec (Diptera: Culicidae) essai de synthèse écologique, in Groupe de recherche sur les insectes piqueurs, mémoires de la société entomologique du Québec. 1980, Université du Québec à Trois-Rivières, Québec, Canada:Québec, Canada. p. 107 pages.
Rocheleau JP, et al. Characterizing environmental risk factors for West Nile virus in Quebec, Canada, using clinical data in humans and serology in pet dogs. Epidemiol Infect. 2017;145(13):2797–807.
Strickman D. Stimuli affecting selection of oviposition sites by Aedes vexans (Diptera: Culicidae): light. J Med Entomol. 1982;19(2):181–4.
Bosak PJ, Crans WJ. The structure and function of the larval siphon and spiracular apparatus of Coquillettidia perturbans. J Am Mosq Control Assoc. 2002;18(4):280–3.
Darold PB, et al. Phenology of Coquillettidia perturbans and Culiseta melanura (Diptera: Culicidae) in East-Central Georgia, USA: implications for the ecology of eastern equine encephalitis virus. J Entomol Sci. 2020;55(2):156–62.
Bosak PJ, Reed LM, Crans WJ. Habitat preference of host-seeking Coquillettidia perturbans (Walker) in relation to birds and eastern equine encephalomyelitis virus in New Jersey. J Vector Ecol. 2001;26(1):103–9.
Andreadis TG, et al. Isolations of Jamestown Canyon virus (Bunyaviridae: Orthobunyavirus) from field-collected mosquitoes (Diptera: Culicidae) in Connecticut, USA: a ten-year analysis, 1997–2006. Vector Borne Zoonotic Dis. 2008;8(2):175–88.
Maire A, Tessier C, Picard L. Analyse éecologique des populations larvaires de moustiques (Diptera: Culicidae) des zones riveraines du fleuve Saint-Laurent. Québec Naturaliste canadien. 1980;105(4):225–41.
Walker ED, Grayson MA, Edman JD. Isolation of Jamestown Canyon and snowshoe hare viruses (California serogroup) from Aedes mosquitoes in western Massachusetts. J Am Mosq Control Assoc. 1993;9(2):131–4.
Watts MJ, et al. The rise of West Nile virus in southern and southeastern Europe: a spatial-temporal analysis investigating the combined effects of climate, land use and economic changes. One health (Amsterdam Netherlands). 2021;13:100315–5.
Kottek MG, Jürgen G, Beck C, Rudolf B, Rubel F. World Map of the Köppen-Geiger climate classification updated. Meteorol Z. 2006;15(3):259–63.
Thornton MM, Shrestha R, Wei Y, Thornton PE, Kao S, Wilson BE. Daymet: daily surface weather data on a 1-km grid for North America, version 4. 2020; Available from: https://doi.org/10.3334/ORNLDAAC/1840.
Hufkens K, et al. An integrated phenology modelling framework in r. Methods Ecol Evol. 2018;9(5):1276–85.
Wood SN. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J Royal Stat Society: Ser B (Statistical Methodology). 2011;73(1):3–36.
AAFC, Agriculture and Agri-Food Canada. ISO 19,131 AAFC Annual crop inventory—data product specifications—revision A. 2019. p. 26.
Canada Go. Government of Canada. Land cover in government of Canada. Ottawa; 2019.
Canada DU. Our work/impact area—Wetlands in Ducks Unlimited Canada, Conserving Canada’s wetlands. 2018.
Canada. Go. Land cover in government of Canada. Ottawa; 2019.
Geomatics PCIPCI. Geomatica. 2019; Available from: https://www.pcigeomatics.com/.
Esri ArcGIS. Desktop: Release 10.6. Environmental Systems Research Institute. 2019; Available from: http://desktop.arcgis.com/en/.
Rue H, Martino S, Chopin N. Approximate bayesian inference for latent gaussian models by using integrated nested Laplace approximations. J Royal Stat Society: Ser B (Statistical Methodology). 2009;71(2):319–92.
Lindgren F, Rue H, Lindström J. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. J Royal Stat Society: Ser B (Statistical Methodology). 2011;73(4):423–98.
Bakka H, et al. Spatial modeling with R-INLA: a review. WIRE Comput Stat. 2018;10(6):e1443.
Musenge E, et al. 2010 in rural north east South Africa. Int J Appl Earth Obs Geoinf. 2013;22(100):86–98. Bayesian analysis of zero inflated spatiotemporal HIV/TB child mortality data through the INLA and SPDE approaches: applied to data observed between 1992.
Myer MH, Campbell SR, Johnston JM. Spatiotemporal modeling of ecological and sociological predictors of West Nile virus in Suffolk County, NY, mosquitoes. Ecosphere. 2017;8(6):e01854.
Myer MH, Johnston JM. Spatiotemporal bayesian modeling of West Nile virus: identifying risk of infection in mosquitoes with local-scale predictors. Sci Total Environ. 2019;650(Pt 2):2818–29.
Fuglstad G-A, et al. Constructing priors that penalize the complexity of gaussian random fields. J Am Stat Assoc. 2019;114(525):445–52.
Simpson D, et al. Penalising model component complexity: a principled, practical approach to constructing priors. Stat Sci. 2017;32(1):1–28.
Fox J, Weisberg S, An. R companion to applied regression, Third edition, ed. Sage. 2019, Thousands Oaks (CA).
Van Niekerk J, Krainski ET, Rustand D, Rue H. A new avenue for Bayesian inference with INLA. arXiv preprint, 2022: p. arXiv:2204.06797.
Gelman A, Hwang J, Vehtari A. Understanding predictive information criteria for bayesian models. Stat Comput. 2014;24(6):997–1016.
Hooten MB, Hobbs NT. A guide to bayesian model selection for ecologists. Ecol Monogr. 2015;85(1):3–28.
Hartig F. DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 0.4.4. 2021; Available from: https://CRAN.R-project.org/package=DHARMa.
Dunn PK, Smyth GK. Randomized quantile residuals. J Comput Graphical Stat. 1996;5(3):236–44.
Laporta GZ, Sallum MAM. Coexistence mechanisms at multiple scales in mosquito assemblages. BMC Ecol, 2014. 14(1).
Hopkins MC et al. Influence of forest disturbance on La Crosse virus risk in southwestern Virginia. Insects, 2020. 11(1).
Bradter U, et al. Identifying appropriate spatial scales of predictors in species distribution models with the random forest algorithm. Methods Ecol Evol. 2013;4:167–74.
Jackson HB, Fahrig L. Are ecologists conducting research at the optimal scale? Glob Ecol Biogeogr. 2015;24(1):52–63.
Stuber EF, Fontaine JJ. How characteristic is the species characteristic selection scale? Glob Ecol Biogeogr. 2019;28(12):1839–54.
Ogden NH, Gachon P. Climate change and infectious diseases: what can we expect? Can Commun Dis Rep, 2019. 45(4).
We wish to thank the Public Health Agency of Canada and the Mitacs Fund for funding this study and awarding a master’s fellowship, the Institut national de santé publique du Québec (Alejandra Irace-Cima) and Quebec’s Ministère de la Santé et des Services sociaux for sharing entomological surveillance data from the provincial West Nile virus surveillance program, and GDG Environnement for sharing some of their mosquito sampling data for southern Quebec.
The work was supported by a MITACS grant established between the University of Sherbrooke and GeoMont.
Open Access provided by Health Canada.
The authors declare that they have no competing interests.
Ethics approval and consent to participate
Consent for publication
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Ludwig, A., Rousseu, F., Kotchi, S.O. et al. Mapping the abundance of endemic mosquito-borne diseases vectors in southern Quebec. BMC Public Health 23, 924 (2023). https://doi.org/10.1186/s12889-023-15773-x
- Public health
- Vector-borne diseases
- Space-time modelling
- Southern Quebec
- Culex pipiens-restuans
- Aedes Ochlerotatus stimulans
- Coquillettidia perturbans
- Aedes vexans