Skip to main content
Download PDF

Spatio-temporal impacts of aerial adulticide applications on populations of West Nile virus vector mosquitoes

Parasites & Vectors volume 14, Article number: 120 (2021) Cite this article

Abstract

Background

Aerial applications of insecticides that target adult mosquitoes are widely used to reduce transmission of West Nile virus to humans during periods of epidemic risk. However, estimates of the reduction in abundance following these treatments typically focus on single events, rely on pre-defined, untreated control sites and can vary widely due to stochastic variation in population dynamics and trapping success unrelated to the treatment.

Methods

To overcome these limitations, we developed generalized additive models fitted to mosquito surveillance data collected from CO2-baited traps in Sacramento and Yolo counties, California from 2006 to 2017. The models accounted for the expected spatial and temporal trends in the abundance of adult female Culex (Cx.) tarsalis and Cx. pipiens in the absence of aerial spraying. Estimates for the magnitude of deviation from baseline abundance following aerial spray events were obtained from the models.

Results

At 1-week post-treatment with full spatial coverage of the trapping area by pyrethroid or pyrethrin products, Cx. pipiens abundance was reduced by a mean of 52.4% (95% confidence intrval [CI] − 65.6, − 36.5%) while the use of at least one organophosphate pesticide resulted in a mean reduction of 76.2% (95% CI − 82.8, − 67.9%). For Cx. tarsalis, at 1-week post-treatment with full coverage there was a reduction in abundance of 30.7% (95% CI − 54.5, 2.5%). Pesticide class was not a significant factor contributing to the reduction. In comparison, repetition of spraying over three to four consecutive weeks resulted in similar estimates for Cx. pipiens and estimates of somewhat smaller magnitude for Cx. tarsalis.

Conclusions

Aerial adulticides are effective for achieving a rapid short-term reduction of the abundance of the primary West Nile virus vectors, Cx. tarsalis and Cx. pipiens. A larger magnitude of reduction was estimated in Cx. pipiens, possibly due to the species’ focal distribution. Effects of aerial sprays on Cx. tarsalis populations are likely modulated by the species’ large dispersal ability, population sizes and vast productive larval habitat present in the study area. Our modeling approach provides a new way to estimate effects of public health pesticides on vector populations using routinely collected observational data and accounting for spatio-temporal trends and contextual factors like weather and habitat. This approach does not require pre-selected control sites and expands upon past studies that have focused on the effects of individual aerial treatment events.

Graphical Abstract

Background

West Nile virus (WNV; genus Flavivirus, family Flaviviridae) causes a potentially fatal, neuroinvasive mosquito-borne disease [1]. It is maintained in an enzootic cycle between birds and mosquitoes [2, 3], predominantly in the genus Culex [4], and can spill over to infect horses and humans, both of which are dead-end hosts vulnerable to disease [5]. Culex (Cx.) tarsalis and Cx. pipiens complex mosquitoes are the primary enzootic and epizootic vectors in California [6, 7]. While 80% of human infections are asymptomatic, clinical manifestations can include acute febrile illness, encephalitis, flaccid paralysis and death [8]. Often the severe form results in long-term physical and mental disabilities [9]. WNV invaded California in 2003 and has become endemic [7]. An average of 238 neuroinvasive cases occur statewide annually, with approximately one-third occurring in the Central Valley, where the landscape is dominated by large-scale agriculture punctuated by cities and small towns [10].

As no human vaccine exists against WNV, prevention of human diseases relies primarily on personal protective measures (i.e. wearing long sleeves, using insect repellent and avoiding the dawn/dust periods when mosquitoes bite) and vector control by local vector control districts or health departments [11, 12]. In periods of epidemic risk when large numbers of WNV-infected Culex mosquitoes are detected near human population centers, large-scale aerial applications of insecticides are utilized to rapidly reduce the abundance of adult mosquitoes and disrupt virus transmission cycles, thereby reducing zoonotic transmission risk [13].

Three main classes of pesticide products have been licensed for use in aerial spray applications in California: pyrethrins, pyrethroids and organophosphates [12, 14, 15]. Pyrethrins are naturally derived insecticides from chrysanthemum flowers (Chrysanthemum cineriaefolium) that inactivate sodium channels in the insect nervous system, resulting in paralysis and death [15]. Pyrethroids are synthetically derived pyrethrins with a similar mode of action and longer half-life. Organophosphates inhibit acetylcholinesterase, affecting neurotransmission and causing uncontrolled nerve activation and death in insects [14].

A standard method for evaluating the efficacy of an aerial spray event compares pre- to post-treatment mosquito trap counts inside the treatment zone versus changes for the same period in an adjacent unsprayed control area [16]. This method, first proposed by Mulla et al. [17], has been adapted to and widely used in evaluating the efficacy of aerial spraying for reducing the abundance of female mosquitoes and has been extended to assess changes in other indicators of risk, namely infection prevalence in mosquitoes, human cases and reported dead birds, with WNV infection [16]. However, reported estimates vary widely, with some studies even indicating occasional increases in trap counts following spray events [18,19,20,21].

Despite its wide use, the assumptions behind Mulla’s formula are often violated, resulting in confounded estimates. First, treatment and control sites are often not independent due to the spatial connectivity of populations with mosquito dispersal and immigration [22]. With the connectivity and potential drift of pesticides via wind, there is the potential that insecticide sprays have wider population impacts than just the targeted spray zone [20]. Second, the difference in pre- to post-trap count ratios in and between areas are not solely due to control measures, but rather are impacted by weather, seasonality in mosquito populations, differential presence of larval breeding sources or simply stochastic variation in trapping success [18, 20, 21]. Overall, Mulla’s formula neglects both the spatio-temporal structure of mosquito populations and the external factors impacting the random volatility of mosquito trapping success.

To overcome the limitations of assessing the efficacy of aerial sprays on the individual spray event basis, we paired long-term surveillance and vector control records (12 years) from Sacramento-Yolo Mosquito and Vector Control District (SYMVCD) in California to capture baseline spatio-temporal mosquito population dynamics and estimate the magnitude and duration of the impacts of aerial sprays on the abundance of Cx. pipiens and Cx. tarsalis, the predominant WNV vectors in California. We chose a generalized additive modeling (GAM) framework to capture the nonlinear population dynamics and associations inherent to mosquito collections.

Methods

Study area

The study area encompasses Sacramento and Yolo counties, California (Fig. 1) which have a combined area of approximately 5,126 km2 and a population of approximately 1.73 million people in 2016 [23]. Sacramento County is 34.12% urban and 65.88% rural, with the majority of urban areas consisting of the concentrated Sacramento urban center and surrounding suburbs. In comparison, Yolo County is 4.61% urban and 95.39% rural, with smaller, more dispersed urban areas [24]. These counties, located in the northern part of California’s Central Valley, are characterized by a Mediterranean climate, with hot, dry summers (July mean temperature: 25.8 °C, May-September mean total rainfall: 3.18 cm) and mild, rainy winters (January mean temperature: 9.6 °C, October–April mean total rainfall: 47.30 cm) [25] and extensive irrigated agriculture, especially rice and row crops, such as tomatoes. SYMVCD, established in 1946 to protect the public from nuisance mosquito biting and mosquito-borne diseases, manages mosquito populations in Sacramento and Yolo counties [26].

Fig. 1
figure 1

Land cover, aerial sprays and mosquito collections (2006–2017) in Sacramento and Yolo counties, California. a Distribution of cultivated crops (primarily rice), urban, and other natural land cover types across the study area. Land cover categories were derived from the 2011 National Land Cover database [34]. Inset highlights the location of these counties in the state of California. b, c Location of zones targeted for aerial treatment applications (b) and of CO2-baited mosquito traps (c) during 2006–2017 in Sacramento and Yolo counties. Each polygon (b) and point (c) represents a single spray or trapping event, respectively. A random shift of ≤ 1 km was applied to trap locations for visualization of repeated trapping at the same location across time

Full size image

Aerial treatments and mosquito collections

SYMVCD provided the spatial polygons (Fig. 1b) and associated data detailing the date, area targeted for spraying, number of consecutive nights of spraying in the same location and pesticide product used for all aerial sprays during the study period (1021 unique nights of spraying during 930 spray events).

Mosquito collection records for CDC CO2-baited EVS traps [27] from SYMVCD for the period 2006–2017 (Fig. 1c) were obtained with permission through the CalSurv Gateway [28], an online database hosting data from California vector control agencies. Any records that indicated trap malfunctions or which ran longer than 1 night were excluded. Each record (n = 24,344) contained latitude, longitude, date, number of traps employed and total female Cx. tarsalis and Cx. pipiens captured. The distribution of traps and spray events by year (Additional file 1: Figure S1) and season (Additional file 2: Figure S2) are presented in the Supplementary Information.

Any records that indicated trap malfunctions or traps which were operated for > 1 night were excluded. For each species separately, we removed the collection reports corresponding to those > 2 standard deviations above the mean in each week to remove the influence of outliers (i.e. large singular deviations from broader abundance trends) on smooth functions subsequently estimated by the models.

All geographic data were projected from geographic to planar coordinates (Albers conic equal-area; EPSG 3310, NAD83) using the rgdal package in R statistical software (version 3.3.2; [29, 30]) for all data processing and analysis.

Covariate development

To isolate the effects of aerial insecticide treatments within our final statistical model, we first developed a set of spatio-temporal and environmental covariates to explain the long- and short-term trends in Cx. tarsalis and Cx. pipiens abundance. Inclusion of these covariates established a counterfactual basis in the models for the expectation in abundance in the absence of controls, leaving the additional terms characterizing aerial insecticide sprays to explain any deviations attributable to the treatments.

For each remaining trap collection (n = 23,707 for Cx. pipiens; n = 23,678 for Cx. tarsalis), we derived a set of temperature variables to capture the effect of weather on trap collections. The mean temperature during the host-seeking period (dusk to dawn) and 30-year monthly average temperature were determined for each collection using 4-km resolution data provided by the PRISM Climate Group [31]. We calculated the deviation in temperature from the monthly average during the host-seeking period to capture activity rates on the night of trapping (i.e. warmer/colder than ‘normal’ resulting in higher/lower mosquito activity and resulting trap counts). As mosquito developmental rates are highly impacted by temperature [32, 33], we also calculated the average temperature during the 2-week period immediately preceding the trap collection to capture the short-term effects of weather on mosquito abundance. Rainfall was not considered because amounts were negligible in the study area during the season when aerial insecticide applications occurred.

To characterize the larval habitat present around a trap location and to incorporate the sharp changes in land use across the study area, we used the 30 × 30-m gridded land cover data from the 2011 National Land Cover database [34]. We used the classification of all pixels within a 5-km radius area surrounding each trap to determine the proportion of three non-overlapping land use categories: urban, cultivated crops and natural (Fig. 1a). Land use categories were chosen to represent larval habitat and bionomics of Cx. pipiens and Cx. tarsalis [35]. The radius was chosen based on known dispersal distances for the species [36, 37]. The ‘urban’ category encompasses all levels of developed land (i.e. structures, roads and constructed materials). The ‘crops’ category represents annual irrigated crops, of which rice is the predominant type in the study area. The remaining classifications were combined to create the ‘natural’ category.

We quantified the spatio-temporal intersections between spray zone polygons and trap locations to quantify the degree to which antecedent spray events impacted mosquito collections. Spatial coverage of each trap was quantified as the average proportion of the area from which a trap collects mosquitoes (‘collection area’) that was covered by aerial treatment zones during the 4 weeks preceding the collection. Using a conservative estimate on Culex flight distance [36,37,38,39] and to account for insecticide drift during application, we used a 5-km radius collection area for both species. The temporal sequence of overlapping sprays during the 4 weeks preceding a collection (modeled as a factor for each unique sequence) was used to capture any lagged effects and impacts of repeated spray events. For each week preceding a collection, the total proportion of the collection area overlapping with the treatment zone was assessed. When multiple treatment zones overlapped with a single collection area in a week, we assumed an additive effect, summing the proportions of overlap from each unique spray up to a maximum of 1.0 that represented complete coverage. We used the average spatial coverage of targeted spray zones during all weeks with at least one overlapping spray event to quantify the spatial impact of sprays for each specified temporal sequence of spraying. Traps > 5 km from all treatment zones had a spatial coverage of zero and corresponded to the reference level of the temporal sequence factor. These traps were included to capture baseline spatio-temporal mosquito dynamics in the absence of aerial treatments. Therefore, the impact of aerial spray events on collections was quantified with a twofold approach, namely with the factor corresponding to the sequence of weeks when spraying overlapped during the preceding 4 weeks and the average proportion of the collection area that overlapped under that sequence.

According to guidelines from the California Department of Public Health, periods of high risk for arbovirus transmission are characterized in part by abnormally high mosquito abundance [12]. In order to capture this dramatic deviation from ‘normal’ abundance that our smoothed modeling framework could not capture, but which precipitated aerial spray events, we also applied a prospective assessment to identify spray events closely following each collection. To account for the time required to respond to a high-risk period, we assessed the presence of overlapping treatment zones with a collection in the following 4 weeks on the weekly scale, similar to the above retrospective assessment of sprays.

To capture potential differences between broad classes of pesticides used (organophosphate vs pyrethrin and pyrethroids combined), we included a binary indicator variable for whether at least one spray event associated with a particular trap collection used an organophosphate. Sample sizes were too small to further investigate differences between pyrethrins and pyrethroids or between individual insecticide products.

We considered time in a variety of ways to capture trends in mosquito abundance in two parts: typical annual seasonality and coarser spatio-temporal trend over the 12-year study period. Using trap collection dates, ‘week’ (number of weeks from the start of the study period; range 1–626) and ‘day’ (range 1–365), variables were created to capture a continuous yearly effect and seasonality, respectively. We interacted the ‘day’ variable with each category of land use (i.e. ‘urban,’ ‘crops’ and ‘natural’) to capture the seasonal trend in these different habitats. Each individual seasonality curve was weighted by the proportion of land use in that category within 5 km of the trap collection to produce a single unified seasonal trend that reflected the specific habitat composition for that collection.

Statistical analysis

We developed GAMs to relate nightly trap counts of female mosquitoes, either Cx. tarsalis or Cx. pipiens, to aerial adulticide applications, adjusted for variation in trap counts due to spatio-temporal mosquito dynamics. We chose GAMs because of the flexible parameterization of smooth functions of covariates to explain spatial and temporal trends [40, 41]. Covariates considered to explain baseline mosquito dynamics were day or week of the year, year, location, land use, 2-week average temperature and nightly deviations from average temperature during trapping, the presence of a spray event in the following 1 to 4 weeks (high-risk period) and pesticide class used in the aerial spray. Either a smooth function or a factor was used in fitting the covariates, with the choice between these forms, along with the spline and basis dimension if a smooth function was chosen, based on model fit and biological relevance. Cyclic cubic regression splines were used to prevent discontinuity between the ends of the smooth representing the seasonal patterns (i.e. between 31 December and 1 January). Thin plate regression splines were chosen for most covariates because they are isotropic and have been shown to be the optimal smoother of any given basis dimension [42]. A cubic regression spline was used in the spatio-temporal surface due to its superior performance over thin plate regression splines for the large amount of observations [43].

We fit negative binomial GAMs using the gam function in R (version 3.3.2; package mgcv) [29, 44] with restricted estimation maximum likelihood (REML) as the smoothing parameter estimation method. We chose a negative binomial function to account for the over-dispersed nature of the trap count data and used backward selection guided by the Akaike information criterion (AIC) [45] to reach our final model. In each model, we included an offset term for the number of traps operated per trapping event. All covariates were included in the initial model, and choices of interactions between covariates entered into the initial model were guided by biological relevance. We used concurvity, a measure of collinearity for smooth functions (range 0–1; 43), and visually examined deviance residuals for consistency in space and time to assess the final model fit.

Using the other covariates to establish the expected abundance of each species in the absence of aerial spraying, we estimated the mean change in predicted abundance across the range of spray regimes observed in the data, using the Bayesian posterior covariance matrix for the parameters that accounted for smoothing parameter uncertainty [43]. We simulated 10,000 random draws from the posterior distribution of the fitted model, a multivariate normal distribution with mean equal to the estimated model coefficients and covariance matrix of the parameters, to predict the abundance of each species across the spatial and temporal sequences of sprays observed in the data. For each draw, we then calculated the mean change in abundance from the baseline no-spray scenario at each point on the spatio-temporal surface, along with the corresponding 95% confidence interval (CI). Estimates of efficacy from the model were compared with those derived from Mulla’s formula [16, 17].

An R script outlining the workflow of parameter development, model fitting and estimating change in abundance across the spatio-temporal surface is presented in Additional file 3: Text S1.

Results

Data overview and model selection

The relative abundance (number of females per trap-night) of Cx. pipiens and Cx. tarsalis varied spatially during the peak WNV season from late June to early October when aerial sprays occurred (Fig. 2). Typically, higher abundance of Cx. pipiens was observed in urban areas, whereas higher abundance of Cx. tarsalis was typically in non-urbanized areas near irrigated agriculture.

Fig. 2
figure 2

Collections of Culex pipiens (a) and Cx. tarsalis (b) during peak West Nile virus (WNV) season. Plus signs (+) indicate the location and number of female mosquitoes per trap-night for each collection during the period when aerial spraying occurred (late June to early October). Colors represent abundance quintiles by species for non-zero collections. A random shift of ≤ 1 km was applied to trap locations to aid visualization of repeated collections at the same locations during the study time frame

Full size image

The final model for each species included an offset for the number of traps run per collection event, and smooth functions of space by time (on the weekly timescale across the 12 years of the study), day of the year by each land use category (‘urban,’ ‘crops’ and ‘natural’), 2-week average temperature, nightly deviations in average temperature during trapping and spatio-temporal impacts of aerial spraying. Our choice of cutoff for removing outliers during model fitting did not significantly change our results. Removing the top 0, 4.5 or 10% of data in each week for each species resulted in only minor shifts in confidence interval widths and magnitude of some estimates, but no change to inference. The resulting smooth functions for each species are illustrated in Fig. 3 (for spatio-temporal surfaces for all years for each species, see Additional file 4: Figure S3; Additional file 5: Figure S4). The construction of smooth functions used in the final models is outlined in Additional file 6: Table S1. All smooth functions were highly significant (P < 0.0001). A random intercept for site location was included to account for repeated collections at the same location, fitted using coefficients penalized by a ridge penalty [46]. We also retained, based on reductions in AIC, parameters for the presence of sprays in the 1 & 4 and 1, 2, & 3 weeks following the trap collection for Cx. tarsalis and Cx. pipiens, respectively (Additional file 7: Table S2). Based on reduction in AIC, only the model for Cx. pipiens retained the term indicating the presence of at least one spray event with an organophosphate pesticide during the previous 4 weeks, as compared to all sprays using combinations of pyrethrin or pyrethroid products (− 48.9% change in abundance for ≥ 1 organophosphate; P < 0.001).

Fig. 3
figure 3

Smooth covariate functions explaining nightly abundance of Cx. pipiens (a) and Cx. tarsalis (b). Smooth functions from final generalized additive models are shown for the spatio-temporal surface (top), seasonality in a fully urban area, seasonality in a fully crop area, seasonality in a fully natural area, deviation (°C) from the 30-year monthly average temperature on the night of trapping and the average temperature (°C) during the 2 weeks prior to trapping (bottom). The shaded region represents the 95% confidence interval for one-dimensional functions. A representative slice of the three-dimensional spatio-temporal surface is presented for 2011 at the midpoint of the typical WNV season (week of 1 August). Spatio-temporal surfaces are plotted on individual axes for each species to resolve the spatial scale

Full size image

The largest magnitude of variability in the baseline abundance for both species was due primarily to the seasonality covariates, followed by the spatio-temporal surface. The temperature covariates contributed the smallest magnitude to establishing abundance, but all smoothed functions were highly significant (P < 0.0001). In all smooth functions, positive estimates correspond to increases in the population, negative estimates correspond to decreases in the population and 0 indicates no modulation in abundance at that covariate value.

The distribution of the final model residuals was right-skewed for both species, indicating the model underestimated extreme trap counts. However, no spatial or temporal pattern remained in the deviance residuals (Additional file 8: Figure S5; Additional file 9: Figure S6). Both models had low estimated concurvity values [43] for the aerial spraying smooth function with the rest of the model parameters (Cx. pipiens: 0.145; Cx. tarsalis: 0.148), indicating that the smooth estimates for the impact of aerial spraying were not confounded by other parameters. In addition, the relatively high deviance-explained value for both models (Cx. pipiens: 44.0%; Cx. tarsalis: 62.3%) indicates good model fit despite the complex dynamics inherent in mosquito populations.

Effects of aerial insecticide treatments

A smooth surface of the spatio-temporal impacts on Cx. pipiens abundance is presented for treatments with only pyrethrin or pyrethroid products (Fig. 4a) or with at least one organophosphate product (Fig. 4b). For Cx. tarsalis, the difference in impact by broad pesticide class was not retained in the final model, and therefore a single smooth is presented (Fig. 4c). Overall, the models estimated a lower magnitude of change in Cx. tarsalis abundance as compared to that in Cx. pipiens abundance. For example, following aerial spraying with full spatial coverage (i.e. 100% coverage of the area within 5 km of the trap), we estimated a mean 1-week change in Cx. pipiens abundance of − 52.4% (95% CI − 65.6, − 36.5%) if all spray events had used pyrethroid or pyrethrin products. If at least one organophosphate product had been used, we estimated a mean change in Cx. pipiens abundance of − 76.2% (95% CI − 82.8, − 67.9%). In contrast, Cx. tarsalis populations with full spatial coverage by aerial sprays showed an estimated mean 1-week post-spraying change in abundance of − 30.7% (95% CI − 54.5, 2.5%) regardless of the product class.

Fig. 4
figure 4

Mean changes in abundance following aerial spraying, as compared to no-spray baseline. Estimates are shown for changes in abundance of Cx. pipiens (a, b) and Cx. tarsalis (c) with respect to antecedent sequence and average spatial coverage of aerial treatments. a, b For Cx. pipiens, estimates are shown for sprays that used only pyrethrin or pyrethroid products (a) and for sprays that utilized an organophosphate product at least once (b). c For Cx. tarsalis, estimates are with any product class. Horizontal axes represent the average proportion of the 5-km buffer surrounding a trap covered by a spray event, and vertical axes represent the temporal sequence of aerial sprays during the 4 weeks preceding the trapping event. On the left vertical axis, the 4-digit sequence indicates presence (1) or absence (0) of sprays in the 1, 2, 3 and 4 weeks (going from right to left along sequence) prior to trapping. The sequences of spray events are ordered from the fewest number and temporally most distant spray events (bottom) to the largest number and temporally closest spray events (top). Estimates are truncated to the range observed with the available data (gray squares indicate points present in dataset). Areas enclosed in a black border represent the portion of the spatio-temporal surface with significant estimates (P < 0.05)

Full size image

For both species, larger reductions in abundances were estimated in areas with higher spatial coverage of aerial sprays (large proportion of spatial overlap) than those on the fringes (low proportion of spatial overlap). Sprays occurring closer in time to collections were generally estimated to result in larger reductions in abundance compared to those that occurred further back in the past. At longer time lags (i.e. 2–4 weeks post-spraying), higher than expected abundance for both species was estimated, with the increase only occurring in Cx. pipiens populations following sprays with pyrethrins and pyrethroids.

The majority of temporal spray sequences lacked data across the full range of spatial overlap (0–100%); we did not estimate the change in abundance for these areas. Data were sparser or lacking for higher spatial coverage during multiple weeks of spraying. For regions of the spatio-temporal surface with data support, the reduction at 1-week post-spraying with full spatial coverage was the largest reduction predicted for Cx. tarsalis. In contrast, results for Cx. pipiens indicated a similar reduction for populations on the fringes of spray events for the preceding 4 weeks and populations with full spatial coverage by aerial sprays 1 week previously (change all pyrethrin/pyrethroids: − 54.3%; 95% CI − 81.0, − 6.2%; change at least one organophosphate: − 77.2%; 95% CI − 90.4, − 53.9%).

Comparison to conventional estimates

As a comparison for our model results, we applied the conventional approach of Mulla’s formula [16, 17] to the combined trapping and control records from 2006 to 2017 to estimate the efficacy of sprays. Considering trapping 1 week before and 1 week after a spray event and using a 5-km buffer around the targeted spray zone as the adjacent comparison area, we were able to calculate the effect for 36 spray events (3.87%) for Cx. pipiens and Cx. tarsalis; the majority of spray events lacked traps in all of the required spatial and temporal locations for the calculation. Most of the estimates for the change in abundance indicated a reduction in trap counts, but estimates varied widely, ranging from complete population elimination (100% decrease) up to increased of 11,000% following a spray event (Additional file 10: Figure S7). Most estimates for Cx. pipiens indicated varying degrees of reduction while those for Cx. tarsalis spanned reductions to increases.

Discussion

This study found that aerial insecticide treatments achieve strong short-term reductions in both Cx. tarsalis and Cx. pipiens populations. Previous studies have assessed the short-term impact of aerial spraying on mosquito abundance and highlighted the volatility of estimates across space and time. In order to overcome the limitations of using single events to estimate the efficacy of aerial spraying on reducing the abundance of WNV vector mosquitoes, we used a large dataset of surveillance and control records together with GAM models. This modeling framework allowed us to establish baseline mosquito adult abundance and identify deviations from expected nightly abundance attributed to aerial spraying (i.e. counterfactual basis) as well as the spatial and temporal impacts of aerial applications. Our results indicate that aerial sprays do achieve population reduction for both Cx. pipiens and Cx. tarsalis with heterogeneity in the magnitude and pattern of reduction between species and pesticide product.

The differences in the magnitude and pattern of the estimated response between species can be attributed partially to the different bionomics of the individual mosquito species. In the study area, Cx. pipiens are predominantly peridomestic with larval habitats limited primarily to backyard sources and stormwater systems in urbanized areas [35, 47]. Thus, areas targeted by aerial insecticide treatments would span a large fraction of any particular population, leaving few adults to repopulate the treated area from proximal unsprayed locations. In contrast, Cx. tarsalis breed in agricultural areas and may disperse into surrounding agricultural and urban areas [35, 48, 49]. This species also has a larger typical dispersal distance and achieves higher population densities than Cx. pipiens [35,36,37]. Aerial insecticide treatments typically target only a small fraction of the total available habitat for Cx. tarsalis, often habitats near urbanized areas, and any effect of aerial sprays could be moderated by immigration of adult Cx. tarsalis from surrounding unsprayed locations, potentially from distant locations [38, 50]. Therefore, the best suppression of Cx. tarsalis populations would be achieved in isolated areas, as has been reported previously [22]. Additionally, repetition of sprays on the weekly scale is less effective at controlling Cx. tarsalis than Cx. pipiens because of the rapid immigration and emergence of new adults from large areas of productive larval habitat.

While our estimated reduction in the abundance of Culex mosquitoes shows some similarity to previous published estimates of aerial spray events from Sacramento and Yolo counties (Table 1), our methodology also accounts for contextual factors, resulting in more robust estimates of the average effect of spraying. Previous estimates exhibited spatial heterogeneity [19]. Utilizing covariates to capture the spatial structure, temperature deviations and varying distribution of larval habitats removed the confounding impact of these factors on our model results. Similarly, previous estimates have varied, in part because Mulla’s formula cannot fully capture spatio-temporal nuances of mosquito population dynamics. For example, Lothrop et al. [20] observed 73% increases in Cx. tarsalis abundances post-spraying despite observing mortality in caged sentinel mosquitoes and large reductions during previous spray events. The authors attributed the estimated increase to the dynamics of Cx. tarsalis populations at the time of the study, particularly the large emergence of Cx. tarsalis following the annual flooding of the nearby wetlands that was not captured by the Mulla’s formula framework and the fact that the sprays were not impacting mosquitoes in the productive larval habitats. Previous estimates also used differing time interval lengths to estimate mosquito abundance before and after spray events. The heterogeneity in these time intervals combined with the heterogeneity of the resulting estimates (Table 1) highlight the need to use consistent time intervals to improve generalizability of estimates between studies. As standardization of the time interval largely depends on the operational capacity of vector control districts for trapping and responding to epidemic conditions, no single recommendation may be feasible across all studies. However, we recommend that mosquito control agencies should keep the time frame consistent across their evaluations to increase comparability of intra-agency control efforts. A similar range of estimates for change in Culex abundance following adulticide treatments has been reported outside California; most are in broad agreement with our findings, although none used Mulla’s formula. In Chicago (Illinois), a reduction of 54% in Cx. pipiens trap counts within the spray zone versus the baseline pre-spray abundance was reported, in contrast with a 153% increase outside the spray zone, following two single-night aerial spray events with a pyrethroid 7 days apart [51]. An average 65.3% reduction in Cx. pipiens/restuans populations was observed within 24 h of truck-mounted applications of a pyrethroid [52]; in contrast, no significant changes in Cx. pipiens abundance were observed following single-night, truck-mounted applications of a pyrethroid in three communities near Boston (Massachusetts) [53]. Up to a 75% reduction in the 2-day counts of female Cx. quinquefasciatus was reported during a month-long period with truck-mounted pyrethroid sprayed 5 days a week in Dubai, United Arab Emirates [54]. Without untreated comparison locations, it is hard to directly compare these results to those of our study.

Table 1 Estimated change in Culex mosquito populations following aerial spray events in California using Mulla’s formula
Full size table

Our model structure most closely resembles the study design used by Elnaiem et al. [21] where a timescale of 1 week before and after a spray event with a pyrethroid pesticide was chosen when assessing mosquito abundance (Table 1). Our estimates for reduction for Cx. pipiens (− 52.4%) and Cx. tarsalis (− 30.7%) are lower than the observed reductions (Cx. pipiens: − 75%; Cx. tarsalis: − 48.7%). While qualitatively similar, the differences in magnitude may be due to differences in analytical methods, shifts from pyrethrin to pyrethroids over time or the longer 12-year time period of our study that could have yielded a more conservative estimate of average spray effects.

Our approach to causal inference using observational data builds on earlier contributions from the fields of environmental science and epidemiology. Mulla’s formula can be considered to be an extension of the Before–After–Control–Impact (BACI) analysis framework. Originally defined by Green [55] and extended and applied by others [56,57,58,59], the BACI analysis framework originated in environmental science literature to distinguish natural variability from the impact of an anthropogenic disturbance, and has been applied to mosquito larvicide evaluations [60, 61]. The BACI methodology compares an impact and at least one separate control location, sampled at various time points before and after the impact, to detect changes in the natural history of the environment due to the impacts [56,57,58]. An analysis of variance (ANOVA) test is used to detect a significant difference in the trajectories before versus after the disturbance in the impact area as compared to the control area. The location and timing of sampling in each trajectory is chosen to ensure each is independent across space and time and increase the evidence that a detected change was attributable to the disturbance itself [56, 59]. Our GAM framework extends BACI using a three-dimensional spatio-temporal function and other covariates to capture the entire spatio-temporal context as a strategy to estimate the expected mosquito abundance in the absence of spraying. The functions also capture trends in the impact over spatial and temporal combinations, while the BACI framework may miss significant changes due to the sampling time frame chosen [62]. Additionally, our methods do not depend on the ANOVA assumptions of independence and homoscedasticity of samples [63], as the spatial and temporal correlation and the non-normal distribution of trap collections are accounted for through the covariates in the negative binomial GAMs.

Our statistical approach for estimating the effects of mosquito control on abundance relies on counterfactual theory that has been applied in the field of epidemiology as a conceptual basis for understanding measures of effect [64,65,66]. Counterfactual theory as a basis for causal inference is premised on the assumption that for any unit being observed, there are multiple potential exposures but only one actually occurs, and outcomes under other alternative exposures exist only as potential outcomes that would have occurred if an alternative exposure had been applied. Because the alternative exposures did not occur, these are contrary to fact, or counterfactual. In this study, our units of study were trapping locations, and we sought to understand the effect of aerial spraying by statistically relating the observed mosquito abundance following spray events to the mosquito abundance in the same place and time that would have been observed in the absence of the spraying. BACI and Mulla’s formula approaches utilize untreated control sites to establish expectations for the unsprayed condition. Our approach instead aims to estimate the counterfactual expectation for mosquito abundance directly at the same place and time using spatio-temporal trends and contextual variables (i.e. weather and land use). This approach offers two clear advantages for estimating the effects of public health pesticide use: (i) it allows for use of rich observational datasets that already exist and which capture pesticide usage in real operational contexts, as opposed to experimental settings that are often closer to ideal conditions; and (ii) it does not rely on pre-selected, untreated control sites, which is helpful because vector management programs are rarely willing to withhold treatments in experimental control sites if their public health action thresholds are met.

The smooth functions associated with land use categories in the final GAMs accurately captured known seasonal and population dynamics. Cultivated crops were the primary source of Cx. tarsalis, with smaller contributions from other, non-urban land types during the warmest months of the year, accurately representing the presence of highly productive larval habitats in clean, recently created water sources characteristic of cultivated crops [35, 48, 67]. Urbanized areas did not produce large numbers of Cx. tarsalis during the WNV season as they contain few suitable larval habitats for this species. The estimated peak in abundance occurred in late July, but remained high through September, capturing the variation in timing of the peak across the years of the study. Populations of Cx. tarsalis in the Sacramento Valley are greatest from July to September [35, 68, 69]. The steep slope of the curve up to the peak mimicked the rapid increase in Cx. tarsalis observed at the start of the planting season [35, 69]. For Cx. pipiens, urbanized areas largely contributed to abundance throughout the year, with additional contributions in natural areas (i.e. non-cultivated croplands) later in the season, reflecting the presence of larval habitats in artificial structures like storm drains or dairy wastewater lagoons [70, 71]. Crops were generally associated with lower Cx. pipiens abundance across the season, reflecting the general lack of suitable high-quality larval habits in these areas.

Temperature anomalies during trapping and the 2-week average antecedent temperature prior to trapping contributed to the overall abundance of both species, albeit relatively weakly in the presence of the other spatio-temporal terms. Their inclusion in the model was required to fully account for mosquito dynamics and night-to-night fluctuations in trapping success. Concordant with previous experiments, extremes in the average temperatures reduced abundance for both species, illustrating the negative impacts on mosquito developmental rates and adult survivorship [32, 33]. The estimated region of positive contribution to abundance for both species (Cx. pipiens: 19.2–25.4 °C; Cx. tarsalis: 18.9–27.0 °C) was narrower than the thermal tolerance of the species, but it contains the observed regions of rapid developmental and high reproduction rates and the typical temperature ranges during the summer. As expected, small anomalies in average temperature on the night of trapping made relatively small contributions to change in the abundance while extreme deviations resulted in much more marked change, highlighting the non-linear relationship underlying temperature and trap success.

It is interesting to note the additional marked reduction in Cx. pipiens abundance when at least one organophosphate product was used, especially as compared to the lack of a similar difference in Cx. tarsalis populations. As mentioned in the preceding text, the class of product used may be more important for Cx. pipiens due to their focal distributions and more limited dispersal [35, 72]. As an aerial spray will likely impact a large proportion of the localized Cx. pipiens population at once, there will be limited immigration from unsprayed segments of the population in the nearby proximity. Therefore, the full effect of an aerial spray is discernable. In contrast, the dispersed nature of Cx. tarsalis populations facilitates rapid immigration from surrounding unsprayed locations [35, 36, 38], thus diluting any difference in effect between product classes; any difference is not discernable against the background population dynamics accounted for in our modeling framework. Another factor contributing to the difference by species could be insecticide resistance, as resistance to pyrethroids and organophosphates has been reported for both species in California [14, 73]. If Cx. pipiens populations in the study area were more resistant to pyrethroids than Cx. tarsalis, as has been previously reported in the Central Valley [74,75,76], this could explain the increased efficacy of organophosphates for Cx. pipiens. However, since we found a stronger effect of pyrethroids on Cx. pipiens as compared to Cx. tarsalis, resistance does not fully explain the observed difference. Alternatively, a single organophosphate spray may be insufficient to produce a marked difference in Cx. tarsalis populations; a repetition may be required. The underlying shape of the smooth function of spatio-temporal impacts of aerial spraying for either species likely differs between product classes and the specific timing and number of different products used, but sparse data prevented us from including an interaction to assess these dynamics.

A potential population rebound effect occurred for both species at more distant time lags from spraying in cases where abundance was estimated to be higher than expected at 2–4 (Cx. pipiens under pyrethrin and pyrethroid sprays) and 3–4 (Cx. tarsalis) weeks post-spray. Appropriately spaced treatments in time may be required to maintain a long-term reduction in population abundance. However, such a rebound does not negate the potential value of aerial treatments for achieving short-term reductions in the abundance of WNV-infected adult mosquitoes during periods of epidemic risk.

The increase in abundance at low spatial coverages of sprays for both species could reflect excito-repellency of pesticides at the fringes of targeted areas. Excito-repellency, a form of behavioral avoidance, combines two forms of sub-lethal exposure that results in mosquito movement away from a chemical source, namely contact excitation (increased activity upon contact) and non-contact spatial repellency [77,78,79]. These non-toxic behavioral impacts of pesticides were first identified in Anopheles mosquitoes in response to DDT and later with insecticide-treated bednets and indoor residual spraying [80,81,82]. Populations of Cx. quinquefasciatus, another species in the Cx. pipiens complex, exhibit strong contact excitation and poor spatial repellency to pyrethroid, organophosphate and carbamate pesticides [83, 84]. No study has investigated excito-repellency in Cx. tarsalis populations. The behavioral avoidance of Culex to pesticides could be pushing mosquitoes out of the spray zones, resulting in an increase in abundance around the fringes of sprays, consistent with our estimates for spatial coverages of < 40% for Cx. pipiens following sprays with pyrethrin or pyrethroids only in the previous week.

A limitation of choosing a GAM framework is that we were only able to capture the average effects of covariates on nightly mosquito trap counts and unable to fully account for large stochastic fluctuations inherent to mosquito populations. However, GAMs easily allowed us to incorporate nonlinear relationships between abundance and covariates without having to constrain relationships with a priori knowledge. In particular, we were able to capture the higher-order relationships and correlation across space and time with the three-dimensional spatio-temporal function. The form of the smooth relationships in the final model did appear to approximate what is observed in nature. Other strengths of our modeling approach are that it takes into account regional differences in mosquito populations, population dynamics and seasonality in different land use types, as well as the impacts of short-term (night) and longer-term (2 week) weather, resulting in robust estimation of the baseline expected abundance in the absence of spray effects, allowing us to isolate the deviations in abundance due to aerial spraying. This counterfactual basis of the model enables estimation in the absence of an independent control. However, even with the large amount of data available, data were still inadequate to estimate the impact of aerial spraying reliably for certain time lags or spatial coverages that were rare or absent in the data. This is primarily due to typical SYMVCD spraying and trapping practices due to logistical and financial constraints. SYMVCD concentrated their mosquito collection efforts near urban areas to maximize the sensitivity for assessing the risk in proximity to human populations while minimizing time and costs associated with large-scale mosquito surveillance. Additionally, in an effort to control mosquitoes in known problem areas (highly productive larval habitats in proximity to the margins of urban areas) and reduce aerial applications over urban areas, the majority of areas receiving repeated sprays across the WNV season are in more rural areas where the mosquito trapping is sparser. This limited the data and the statistical power to quantify the full range of spatial overlap with aerial sprays.

We were also unable to account explicitly for drift outside the target zones during sprays, as has been previously described [20]. Our use of a continuous variable to measure spray coverage within the 5-km radius collection areas surrounding each trap partially accounts for this effect. As such, however, we are unable to fully parse out the effect of aerial spraying on populations outside aerial spray zones and limited our analysis to only assess spatial coverage of traps within targeted spray zones.

Additionally, we were unable to estimate the relative effects of different lengths of multi-night spray events (1 vs 2 vs 3 consecutive nights) due to data limitations and our analytical choice to aggregate all sprays on the weekly scale to achieve the balance between robust estimates and operationally relevant information for vector control districts. These limitations of observational studies such as this one could be addressed in future experimental field trials.

Conclusions

Aerial adulticides were shown to achieve short-term reductions in the abundance of the primary WNV vectors Cx. tarsalis and Cx. pipiens. A greater reduction was estimated for Cx. pipiens, likely due to its focal distribution in urbanized areas and limited dispersal. The use of organophosphate products versus a combination of pyrethrins and pyrethroids increased the magnitude of reduction estimated for Cx. pipiens while the difference by broad insecticide class was not significant for Cx. tarsalis. The effects of aerial sprays on Cx. tarsalis populations were likely moderated by the broad dispersal ability of this species, its large population sizes and the vast expanses of productive larval habitat in the study area. Therefore, the best control of Cx. tarsalis would appear to be achieved in areas with isolated or highly spatially segmented populations. For both species, aerial spraying reduced abundance at high spatial coverage while reductions were also estimated at lower spatial coverage, at albeit greatly reduced magnitudes, indicating that aerial sprays had some impacts beyond the target zone. There was also evidence for population rebounds at periods of 2–4 weeks post-spraying. Our modeling approach allowed us to utilize observational data to isolate aerial treatment effects while taking into account contextual factors, such as spatio-temporal relationships, weather and habitat, that contribute to stochastic variation in nightly trap counts. This is an important advance that complements experimental trials and expands upon conventional observational approaches that summarize population changes following aerial treatments at individual time points. Further work should expand upon these methods to estimate the change in WNV transmission potential and resulting human infections following aerial spray events.

Availability of data and materials

The data that support the findings of this study were obtained from the Sacramento-Yolo Mosquito and Vector Control District. These data were used with permission for the current study and are not publicly available. Data are however available from the authors upon reasonable request and with permission from the Sacramento-Yolo Mosquito and Vector Control District.

Abbreviations

GAM:

Generalized additive model

SYMVCD:

Sacramento-Yolo Mosquito and Vector Control District

WNV:

West Nile virus

References

  1. Hayes EB, Sejvar JJ, Zaki SR, Lanciotti RS, Bode AV, Campbell GL. Virology, pathology, and clinical manifestations of West Nile virus disease. Emerg Infect Dis. 2005;11(8):1174–9.

    Article  PubMed  PubMed Central  Google Scholar 

  2. McLean RG, Ubico SR, Docherty DE, Hansen WR, Sileo L, McNamara TS. West Nile virus transmission and ecology in birds. Ann N Y Acad Sci. 2001;951:54–7.

    Article  CAS  PubMed  Google Scholar 

  3. Kilpatrick AM, LaDeau SL, Marra PP. Ecology of West Nile virus transmission and its impact on birds in the western hemisphere. Auk. 2007;124(4):1121–36.

  4. Turell MJ, Dohm DJ, Sardelis MR, O Guinn ML, Andreadis TG, Blow JA. An update on the potential of North American mosquitoes (Diptera : Culicidae) to transmit West Nile virus. J Med Entomol. 2005;42(1):57–62.

    Article  PubMed  Google Scholar 

  5. Kramer LD, Styer LM, Ebel GD. A global perspective on the epidemiology of West Nile virus. Annu Rev Entomol. 2008;53:61–81.

    Article  CAS  PubMed  Google Scholar 

  6. Goddard LB, Roth AE, Reisen WK, Scott TW. Vector competence of California mosquitoes for West Nile virus. Emerg Infect Dis. 2002;8(12):1385–91.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Reisen WK, Lothrop HD, Chiles RE, Madon MB, Cossen C, Woods L, et al. West Nile virus in California. Emerg Infect Dis. 2004;10:1369–78.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Mostashari F, Bunning ML, Kitsutani PT, Singer DA, Nash D, Cooper MJ, et al. Epidemic West Nile encephalitis, New York, 1999: results of a household-based seroepidemiological survey. Lancet. 2001;358(9278):261–4.

    Article  CAS  PubMed  Google Scholar 

  9. Hughes JM, Wilson ME, Sejvar JJ. The long-term outcomes of human West Nile virus infection. Clin Infect Dis. 2007;44(12):1617–24.

    Article  Google Scholar 

  10. California Department of Public Health. Human West Nile virus activity, California, 2003–2019. 2019. http://westnile.ca.gov/. Accessed 1 Feb 2021.

  11. Gubler DJ, Campbell GL, Nasci R, Komar N, Petersen L, Roehrig JT. West Nile virus in the United States: guidelines for detection, prevention, and control. Viral Immunol. 2000;13(4):469–75.

    Article  CAS  PubMed  Google Scholar 

  12. California Department of Public Health, Mosquito and Vector Control Association of California, University of California. California mosquito-borne virus surveillance & response plan. 2019. https://westnile.ca.gov/resources_reports.php?resource_category_id=9. Accessed 1 Feb 2021.

  13. Rose RI. Pesticides and public health: integrated methods of mosquito management. Emerg Infect Dis. 2001;7(1):17–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. California Department of Health Services Vector-Borne Disease Section. Overview of mosquito control practices in California. 2005. https://www.cdph.ca.gov/Programs/CID/DCDC/CDPH%20Document%20Library/OverviewofMosquitoControlinCA.pdf. Accessed 1 Feb 2021.

  15. Schleier JJ, Peterson RKD. Pyrethrins and pyrethroid insecticides. In: Lopez O, Fernfmdez-Bolafios JG, editors. Green trends in insect control. London: Royal Society of Chemisty; 2011. p. 94–131.

    Google Scholar 

  16. Reisen WK. Using “Mulla’s Formula” to estimate percent control. In: Atkinson PW, editor. Vector biology, ecology, and control. New York: Springer SBM; 2010. p. 127–38.

    Chapter  Google Scholar 

  17. Mulla MS, Norland RL, Fanara DM, Darwezeh HA, McKeen DW. Control of chironomid midges in recreational lakes. J Econ Entomol. 1971;64:300–7.

    Article  CAS  Google Scholar 

  18. Nielsen CF, Reisen WK, Armijos V, Wheeler S, Kelley K, Brown D. Impact of climate variation and adult mosquito control on the West Nile virus epidemic in Davis, California during 2006. Proc Pap Mosq Vector Control Assoc Calif. 2007;75:125–30.

    Google Scholar 

  19. Macedo PA, Nielsen CF, Reed M, Kelley K, Reisen WK, Goodman GW, et al. An evaluation of the aerial spraying conducted in reponse to West Nile virus activity in Yolo county. Proc Pap Mosq Vector Control Assoc Calif. 2007;75:107–14.

    Google Scholar 

  20. Lothrop H, Lothrop B, Palmer M, Wheeler S, Gutierrez A, Lothrop H, et al. Evaluation of pyrethrin aerial ultra-low volume applications for adult Culex tarsalis control in the desert environments of the Coachella Valley, Riverside county California. J Am Mosq Control Assoc. 2007;23(4):405–19.

    Article  CAS  PubMed  Google Scholar 

  21. Elnaiem DA, Kelley K, Wright S, Laffey R, Yoshimura G, Reed M, et al. Impact of aerial spraying of pyrethrin insecticide on Culex pipiens and Culex tarsalis (Diptera: Culicidae) abundance and West Nile virus infection rates in an urban/suburban area of Sacramento County California. J Med Entomol. 2008;45(4):751–7.

    Article  PubMed  Google Scholar 

  22. Reisen WK, Milby MM, Reeves WC, Eberle MW, Meyer RP, Schaefer CH, et al. Aerial adulticiding for the suppression of Culex tarsalis in Kern County, California, using low volume propoxur: 2. Impact on natural populations in foothill and valley habitats. J Am Mosq Control Assoc. 1985;1(2):154–63.

    CAS  PubMed  Google Scholar 

  23. U.S. Census Bureau. State and County QuickFacts. 2019. https://www.census.gov/quickfacts. Accessed 14 Oct 2019.

  24. United State Census Bureau. Percent urban and rural in 2010 by state and county. https://www.census.gov/geo/reference/ua/urban-rural-2010.html. Accessed 14 Oct 2019.

  25. Cline G, Neigher A, Bellinder A. Climate of Sacramento, California. Sacramento: National Weather Service Office. Sacramento; 2010.

    Google Scholar 

  26. Sacramento-Yolo Mosquito & Vector Control District. Mission & vision. Elk Grove: Sacramento-Yolo Mosquito & Vector Control District; 2018. http://www.fightthebite.net/about/about-us/. Accessed 5 Oct 2019.

  27. Newhouse VF, Chamberlain RW, Johnson JG, Sudia WD. Use of dry ice to increase mosquito catches of the CDC miniature light trap. Mosq News. 1966;26(1):30–5.

    Google Scholar 

  28. California Vectorborne Disease Surveillance System (CalSurv). 2018. https://vectorsurv.org/. Accessed 1 Feb 2021.

  29. R Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2020. https://www.r-project.org/.

  30. Bivand R, Keitt T, Rowlingson B. rgdal: bindings for the “geospatial” data abstraction library; 2019. https://cran.r-project.org/package=rgdal.

  31. PRISM Climate Group, Oregon State University. Daily mean and monthly average temperature datasets. 2018. http://prism.oregonstate.edu. Accessed 12 May 2018.

  32. Ciota AT, Matacchiero AC, Kilpatrick AM, Kramer LD. The effect of temperature on life history traits of Culex mosquitoes. J Med Ent. 2014;51(1):55–62.

    Article  Google Scholar 

  33. Reisen WK. Effect of temperature on Culex tarsalis (Diptera: Culicidae) from the Coachella and San Joaquin valleys of California. J Med Entomol. 1995;32(5):636–45.

    Article  CAS  PubMed  Google Scholar 

  34. Multi-Resolution Land Characteristics Consortium. NLCD Land Cover (CONUS) database. 2011. https://www.mrlc.gov/data?f%5B0%5D=category%3Alandcover&f%5B1%5D=year%3A2011 . Accessed 20 May 2018.

  35. Reisen WK, Reeves WC. Bionomics and ecology of Culex tarsalis and other potential mosquito vector species. In: Reeves WC, editor. Epidemiology and control of mosquito-borne arboviruses in California, 1943–1987. Sacramento: California Mosquito and Vector Control Association; 1990. p. 254–329.

    Google Scholar 

  36. Bailey SF, Eliason DA, Hoffmann BL. Flight and dispersal of mosquito Culex tarsalis Coquillett in Sacramento Valley of California. Hilgardia. 1965;37(3):73–113.

    Article  Google Scholar 

  37. Reisen WK, Milby MM, Meyer RP, Pfuntner AR, Spoehel J, Hazelrigg JE, et al. Mark-release-recapture studies with Culex mosquitoes (Diptera: Culicidae) in southern California. J Med Entomol. 1991;28(3):357–71.

    Article  CAS  PubMed  Google Scholar 

  38. Dow RP, Reeves WC, Bellamy RE. Dispersal of female Culex tarsalis into a larvicided area. Am J Trop Med Hyg. 1965;14(4):656–70.

    Article  Google Scholar 

  39. Reisen WK, Milby MM, Meyer RP. Population dynamics of adult Culex mosquitoes (Diptera: Culicidae) along the Kern River, Kern county, California, in 1990. J Med Ent. 1992;29(3):531–43.

    Article  CAS  Google Scholar 

  40. Hastie T, Tibshirani R. Varying-coefficient models. J R Stat Soc Ser B-Methodol. 1993;55(4):757–96.

    Google Scholar 

  41. Hastie T, Tibshirani RJ. Generalized additive models. Stat Sci. 1986;1(3):297–318.

    Google Scholar 

  42. Wood SN. Thin plate regression splines. J R Stat Soc Ser B. 2003;65(1):95–114.

    Article  Google Scholar 

  43. Wood SN. Generalized additive models: an introduction with R. Boca Raton: Chapman & Hall/CRC; 2006.

    Book  Google Scholar 

  44. Wood SN. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J R Stat Soc Ser B-Stat Methodol. 2011;73:3–36.

    Article  Google Scholar 

  45. Aike HAI. A new look at the statistical model identification. IEEE Trans Automat Contr. 1974;19(6):716–23.

    Article  Google Scholar 

  46. Wood SN. Fast stable direct fitting and smoothness selection for generalized additive models. J R Stat Soc Ser B-Stat Methodol. 2008;70:495–518.

    Article  Google Scholar 

  47. Hopkins CC, Hollinger FB, Johnson RF, Dewlett HJ, Newhouse VF, Chamberlain RW. The epidemiology of St. Louis encephalitis in Dallas, Texas, 1966. Am J Epidemiol. 1975;102(1):1–15.

    Article  CAS  PubMed  Google Scholar 

  48. Wekesa JW, Yuval B, Washino RK. Spatial distribution of adult mosquitoes (Diptera:Culicidae) in habitats associated with the rice agroecosystem of northern California. J Med Entomol. 1996;33(3):344–50.

    Article  CAS  PubMed  Google Scholar 

  49. Mitchell CJ, Kilpatrick JW, Hayes RO, Curry HW. Effects of ultra-low volume applications of malathion in Hale county, Texas II: mosquito populations in treated and untreated areas. J Med Ent. 1970;7(1):85–91.

    Article  CAS  Google Scholar 

  50. Reisen WK, Yoshimura G, Reeves WC, Milby MM, Meyer RP. The impact of aerial applications of ultra-low volume adulticides on Culex tarsalis populations (Diptera: Culicidae) in Kern County, California, USA, 1982. J Med Entomol. 1984;21(5):573–85.

    Article  CAS  PubMed  Google Scholar 

  51. Mutebi JP, Delorey MJ, Jones RC, Plate DK, Gerber SI, Gibbs KP, et al. The impact of adulticide applications on mosquito density in Chicago, 2005. J Am Mosq Control Assoc. 2011;27(1):69–76.

    Article  PubMed  Google Scholar 

  52. Clifton ME, Xamplas CP, Nasci RS, Harbison J. Gravid Culex pipiens exhibit a reduced susceptibility to ultra-low volume adult control treatments under field conditions. J Am Mosq Control Assoc. 2019;35(4):267–78.

    Article  PubMed  Google Scholar 

  53. Reddy MR, Spielman A, Lepore TJ, Henley D, Kiszewski AE, Reiter P. Efficacy of resmethrin aerosols applied from the road for suppressing Culex vectors of West Nile virus. Vector-Borne Zoonotic Dis. 2006;6(2):117–27.

    Article  PubMed  Google Scholar 

  54. Holmes PR. A study of population changes in adult Culex quinquefasciatus Say (Diptera: Culicidae) during a mosquito control programme in Dubai, United Arab Emirates. Ann Trop Med Parasitol. 1986;80(1):107–16.

    Article  CAS  PubMed  Google Scholar 

  55. Green RH. Sampling design and statistical methods for environmental biologists. Chichester: Wiley Interscience; 1979.

    Google Scholar 

  56. Underwood AJ. Beyond BACI: the detection of environmental impacts on populations in the real, but variable, world. J Exp Mar Bio Ecol. 1992;161:145–78.

    Article  Google Scholar 

  57. Stewart-Oaten A, Murdoch WW, Parker KR. Environmental impact assessment: “pseudoreplication” in time? Ecology. 1986;67(4):929–40.

    Article  Google Scholar 

  58. Stewart-Oaten A, Bence JR. Temporal and spatial variation in environmental impact assessment. Ecol Monogr. 2001;71(2):305–39.

    Article  Google Scholar 

  59. Underwood AJ. On beyond BACI: sampling designs that might reliably detect environmental disturbances. Ecol Appl. 1994;4(1):3–15.

    Article  Google Scholar 

  60. Rochlin I, Iwanejko T, Dempsey ME, Ninivaggi DV. Geostatistical evaluation of integrated marsh management impact on mosquito vectors using before-after-control-impact (BACI) design. Int J Health Geogr. 2009;8(1):1–20.

    Article  Google Scholar 

  61. Wolfram G, Wenzl P, Jerrentrup H. A multi-year study following BACI design reveals no short-term impact of Bti on chironomids (Diptera) in a floodplain in Eastern Austria. Environ Monit Assess. 2018;190(12):1–7.

    Article  CAS  Google Scholar 

  62. Smokorowski KE, Randall RG. Cautions on using the Before-After-Control-Impact design in environmental effects monitoring programs. Facets. 2017;2(1):212–32.

    Article  Google Scholar 

  63. Stewart-Oaten A, Bence JR, Osenberg CW. Assessing effects of unreplicated perturbations: no simple solutions. Ecology. 1992;73(4):1396–404.

    Article  Google Scholar 

  64. Lewis D. Causation. J Philos. 1973;70:556–67.

    Article  Google Scholar 

  65. Rubin DB. Comment: Neyman (1923) and causal inference in experiments and observational studies. Stat Sci. 1990;5:472–80.

    Article  Google Scholar 

  66. Greenland S, Rothman KJ, Lash TL. Measures of effect and measures of association. In: Rothman KJ, Greenland S, Lash TL, editors. Modern epidemiology. 3rd ed. Philadelphia: Lippincott Williams and Wilkins; 2008. p. 51–70.

    Google Scholar 

  67. Reisen WK, Meyer RP, Shields J, Arbolante C. Population ecology of preimaginal Culex tarsalis (Diptera: Culicidae) in Kern County, California. J Med Entomol. 1989;26(1):10–22.

    Article  CAS  PubMed  Google Scholar 

  68. Loomis EC, Meyers EG. California encephalitis surveillance program: mosquito population measurement and ecologic considerations. Am J Epidemiol. 1960;71(3):378–88.

    Article  CAS  Google Scholar 

  69. Barker CM, Eldridge BF, Reisen WK. Seasonal abundance of Culex tarsalis and Culex pipiens complex mosquitoes (Diptera: Culicidae) in California. J Med Ent. 2010;47(5):759–68.

    Article  Google Scholar 

  70. Loomis EC, Eide RN, Caton JR, Merritt DA. Physical and cultural control: Mosquito and fly problems in dairy waste-water systems. Calif Agric. 1980;34(3):37–8.

    Google Scholar 

  71. Su TY, Webb JR, Meyer RR, Mulla MS. Spatial and temporal distribution of mosquitoes in underground storm drain systems in Orange County, California. J Vector Ecol. 2003;28(1):79–89.

    PubMed  Google Scholar 

  72. Reisen WK. The contrasting bionomics of Culex mosquitoes in Western North America. J Am Mosq Control Assoc. 2012;28(4):82–91.

    Article  PubMed  Google Scholar 

  73. California Department of Public Health Vector-Borne Disease Section. California mosquito pesticide resistance summary. 2005. http://westnile.ca.gov/downloads.php?download_id=3137&filename=Pesticide_Resistance.pdf. Accessed 1 Feb 2021.

  74. Reed M, Macedo PA, Brown D. Increased tolerance to permethrin in Culex pipiens complex population from Sacramento County, California. Proc Pap Mosq Vector Control Assoc Calif. 2012;80:56–8.

    Google Scholar 

  75. Reed M, Macedo PA, Goodman G, Brown D. Tough mosquitoes—why they should be everyone’s problem. Proc Pap Mosq Vector Control Assoc Calif. 2013;81:75–9.

    Google Scholar 

  76. Sorensen MA, Stevenson JA. An adulticide resistance “heat map” for Culex pipiens and Culex tarsalis (Diptera: Culididae) in Placer County. Proc Pap Mosq Vector Control Assoc Calif. 2015;83:8–11.

    Google Scholar 

  77. Dethier VG, Browne LB, Smith CW. The designation of chemicals in terms of the responses they elicit from insects. J Econ Entomol. 1960;53(1):134–6.

    Article  CAS  Google Scholar 

  78. Roberts DR, Chareonviriyaphap T, Harlan HH, Hshieh P. Methods ot testing and analyzing excito-repellency responses of malaria vectors to insecticides. J Am Mosq Control Assoc. 1997;3(1):13–7.

    Google Scholar 

  79. Miller JR, Siegert PY, Amimo FA, Walker ED. Designation of chemicals in terms of the locomotor responses they elicit from insects: an update of Dethier et al. (1960). J Econ Entomol. 2009;102(6):2056–60.

    Article  CAS  PubMed  Google Scholar 

  80. Lindsay SW, Adiamah JH, Miller JE, Armstrong JRM. Pyrethroid-treated bednet effects on mosquitoes of the Anopheles gambiae complex in The Gambia. Med Vet Entomol. 1991;5:477–83.

    Article  CAS  PubMed  Google Scholar 

  81. Kennedy JS. The excitant and repellent effects on mosquitos of sub-lethal contacts with DDT. Bull Entomol Res. 1945;37(4):593–607.

    Article  Google Scholar 

  82. Davidson G. Experiments on the effect of residual insecticides in houses against Anopheles gambiae and A. funestus. Bull Entomol Res. 1953;44:231–54.

    Article  CAS  Google Scholar 

  83. Boonyuan W, Bangs MJ, Grieco JP, Tiawsirisup S, Prabaripai A, Chareonviriyaphap T. Excito-repellent responses between Culex quinquefasciatus permethrin susceptible and resistant mosquitoes. J Insect Behav. 2016;29:415–31.

    Article  Google Scholar 

  84. Sathantriphop S, Ketavan C, Prabaripai A, Visetson S, Bangs MJ, Akratanakul P, et al. Susceptibility and avoidance behavior by Culex quinquefasciatus Say to three classes of residual insecticides. J Vector Ecol. 2006;31(2):266–74.

    Article  PubMed  Google Scholar 

  85. Macedo PA, Schleier JJ, Reed M, Kelley K, Goodman GW, Brown DA, et al. Evaulation of efficacy and human health risk of aerial ultra-low volume applications of pyrethrins and piperonyl butoxide for adult mosquito management in response to West Nile virus activity in Sacramento county. California J Am Mosq Control Assoc. 2010;26(1):57–66.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We would like to thank Ruben Rosas, Marcia Reed, Sarah Wheeler, Samer Elkashef and Gary Goodman from the Sacramento-Yolo Mosquito & Vector Control District for providing trapping and aerial treatment data used in this study and for their insight on mosquito control practices. We also thank the two anonymous reviewers for their helpful comments which improved the manuscript.

Funding

KH acknowledges funding support from the Floyd & Mary Schwall Fellowship in Medical Research at UC Davis, and KH and CB acknowledge support from the Pacific Southwest Center of Excellence in Vector-Borne Diseases funded by the U.S. Centers for Disease Control and Prevention (Cooperative Agreement 1U01CK000516). KH acknowledges support by the National Center for Advancing Translational Sciences, National Institutes of Health, through grant number UL1 TR001860 and linked award TL1 TR001861.

Author information

Authors and Affiliations

  1. Department of Pathology, Microbiology, and Immunology, School of Veterinary Medicine, University of California, Davis, CA, 95616, USA

    Karen M. Holcomb & Christopher M. Barker

  2. Institute for Health Metrics and Evaluation, University of Washington, Seattle, WA, 98121, USA

    Robert C. Reiner

Authors
  1. Karen M. Holcomb
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Robert C. Reiner
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Christopher M. Barker
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

KH performed the analysis and wrote the initial draft of the manuscript. CB conceived the project idea, supervised the research, and contributed to revisions of the manuscript. RR provided statistical guidance. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Christopher M. Barker.

Ethics declarations

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

: Figure S1. Location of CO2-baited mosquito trapping events and aerial spray events stratified by year (2006–2017). A random jitter of ≤ 1 km was applied to trapping locations for visualization of repeated events at the same site. Each spray event polygon represents the area targeted during a single aerial spray application.

Additional file 2

: Figure S2. Location of CO2-baited mosquito trapping events and aerial spray events stratified by season. Season defined into three-month intervals. A random jitter of ≤ 1 km was applied to trapping locations for visualization of repeated events at the same site. Each spray event polygon represents the area targeted during a single aerial spray application.

Additional file 3

: Text S1. R script outlining our workflow of covariate development, GAM model fitting and estimated change in abundance.

Additional file 4

: Figure S3. Spatial surface for Cx. pipiens at midpoint of the typical WNV season (week of 1 Aug). Surface presented for each year (2006–2017) reflects the relative abundance of the species and is a slice from the three-dimensional spatio-temporal smoothed function. Contours applied for visualization of estimates.

Additional file 5

: Figure S4. Spatial surface for Cx. tarsalis at midpoint of the typical WNV season (week of 1 Aug). Surface presented for each year (2006–2017) reflects the relative abundance of the species and is a slice from the three-dimensional spatio-temporal smoothed function. Contours applied for visualization of estimates.

Additional file 6

: Table S1. Smooth functions included in the final GAMs for both Cx. tarsalis and Cx. pipiens. Table outlining the construction of smooth function used in the final models, including spline type and basis dimensions chosen.

Additional file 7

: Table S2. Change (%) in nightly abundance from expected for collections preceding aerial spraying. Table indicating the final model estimates for change in expected trap-counts for collections in the 1 to 4 weeks preceding an aerial spray event.

Additional file 8

: Figure S5. Spatial and temporal distribution of model deviance residuals for Cx. pipiens for 2006–2017. Residuals presented spatially at the associated trapping location (random jitter of ≤ 1 km was applied for visualization of repeated events).

Additional file 9

: Figure S6. Spatial and temporal distribution of model deviance residuals for Cx. tarsalis for 2006–2017. Residuals presented spatially at the associated trapping location (random jitter of ≤ 1 km was applied for visualization of repeated events).

Additional file 10

: Figure S7. Estimated percentage change in Cx. pipiens and Cx. tarsalis populations with Mulla’s formula. Change estimated for the 36 aerial sprays in Sacramento and Yolo counties, California (2006–2017) with associated trap collections within the targeted zone (treated) and an adjacent 5-km buffer (control) within 1 week before and 1 week following spraying.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Holcomb, K.M., Reiner, R.C. & Barker, C.M. Spatio-temporal impacts of aerial adulticide applications on populations of West Nile virus vector mosquitoes. Parasites Vectors 14, 120 (2021). https://doi.org/10.1186/s13071-021-04616-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13071-021-04616-6

Keywords

Parasites & Vectors

ISSN: 1756-3305

Contact us