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 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–37]. Aerial insecticide treatments typically target only a small fraction of the total available habitat for Cx. tarsalis, often 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 for Cx. tarsalis populations would be achieved in isolated areas, as has been reported previously . 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 abundance of Culex mosquitoes show some similarity to previous published estimates of aerial spray events from Sacramento and Yolo counties (Table 1), our methodology also accounted for contextual factors, resulting in more robust estimates of the average effect of spraying. Previous estimates exhibited spatial heterogeneity . 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.  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 use differing time interval lengths to estimate mosquito abundance before and after spray events. The heterogeneity in these time intervals combined with the heterogeneity of resulting estimates (Table 1) highlights 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 timeframe 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 have 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 vs. 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 seven days apart . An average 65.3% reduction in Cx. pipiens/restuans populations was observed within 24 hours of truck-mounted applications of a pyrethroid . 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 . Up to 75% reduction in the two-day counts of female Cx. quinquefasciatus was reported during a month-long period with truck-mounted pyrethroid sprayed five days a week in Dubai, United Arab Emirates . Without untreated comparison locations, it is hard to directly compare these results to those of our study.
Our model structure most closely compares to the study design used by Elnaiem et al.  where a timescale of one 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 twelve-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 upon earlier contributions from the fields of environmental science and epidemiology. Mulla’s formula can be considered an extension of the Before/After and Control/Impact (BACI) analysis framework. Originally defined by Green  and extended and applied by others [56–59], BACI 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–58]. An ANOVA test is used to detect a significant difference in the trajectories before vs. after the disturbance in the impact area as compared to the control area. Location and timing of sampling in each is chosen to ensure each are 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 entire spatio-temporal context as a way 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 timeframe chosen . Additionally, our methods do not depend on the ANOVA assumptions of independence and homoscedasticity of samples , 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-66]. Counterfactual theory as a basis for causal inference is premised on the idea 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 are trapping locations, and we are seeking 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 spray. 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: (1) it allows for use of rich observational data sets that exist already and capture pesticide usage in real operational contexts, as opposed to experimental settings that are often closer to ideal conditions and (2) 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-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 (aka 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 generally had lower Cx. pipiens abundance across the season, reflecting the general lack of high quality suitable larval habits in these areas.
Temperature anomalies during trapping and the two-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 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 abundance while extreme deviations result 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 above, 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 localize 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 have 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–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 where abundance was estimated to be higher than expected two to four (Cx. pipiens under pyrethrin and pyrethroid sprays) and three to four (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; contact excitation (increased activity upon contact) and non-contact spatial repellency [77–79]. These non-toxic behavioral impacts of pesticides were first identified in Anopheles mosquitoes in response to DDT and later with insecticide-treated bed nets and indoor residual spraying [80–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 < 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, and the impacts of short-term (night) and longer-term (two 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 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 collections 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 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 . Our use of a continuous variable to measure spray coverage within the 5km-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 (one vs two vs three 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.