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 (Jul mean temperature: 25.8 °C, May-Sep mean total rainfall: 3.18 cm) and mild, rainy winters (Jan mean temperature: 9.6 °C, Oct-Apr mean total rainfall: 47.30 cm) (25) and extensive irrigated agriculture, especially rice and row crops such as tomatoes. Sacramento-Yolo Mosquito and Vector Control District (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).
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 (1,021 unique nights of spraying during 930 spray events).
Mosquito collection records for CDC CO2-baited EVS traps (28) from SYMVCD for the years 2006-2017 (Fig 1c) were obtained with permission through the CalSurv Gateway (29), an online database hosting data from California vector control agencies. Any records that indicated trap malfunctions or which unfeasibly ran longer than one 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. Distribution of traps and spray events by year (Additional file 1) and season (Additional file 2) are presented in the Supplementary Information.
Any records that indicated trap malfunctions or which were operated for more than one night were excluded. For each species separately, we removed the collection reports corresponding to the top 5% of trap counts 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; (30,31)) for all data processing and analysis.
Covariate Development
To isolate the effects of aerial insecticide treatments within out 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 the absence of control, 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 (32). 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 (33,34), we also calculated the average temperature during the two-week period immediately preceding the trap collection to capture 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 30x30m gridded land cover data from the 2011 National Land Cover database (27). We used the classification of all pixels within a 5km 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 the 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 which are predominantly rice 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 four weeks preceding the collection. Using a conservative estimate on Culex flight distance (36–39) and to account for insecticide drift during application, we used a 5km radius collection area for both species. The temporal sequence of overlapping sprays during the four 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 > 5km 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 two-fold approach, namely with the factor corresponding to the sequence of weeks when spraying overlapped during the preceding four 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 that 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 four 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 spatial-temporal trend over the twelve-year study period. Using the 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 generalized additive models (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, two-week average temperature, and nightly deviations from average temperature during trapping, the presence of a spray event in the following one to four weeks (high risk period), and pesticide class used in the aerial spray. Either a smooth function or a factor were used in fitting the covariates with the choice between these forms, along with 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 (aka between Dec 31 and Jan 1). 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) (30,44) with restricted estimation maximum likelihood (REML) as the smoothing parameter estimation method. We choses a negative binomial function to account for the over-dispersed nature of trap count data. We used backward selection guided by AIC (45) to reach our final model. In each model, we included an on 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. Estimates of efficacy from the model were compared with those derived from Mulla’s formula (16,17).