Insights and Challenges for the Development of Mosquito Control Action Thresholds Using Historical Mosquito Surveillance and Climate Datasets

Mosquito control has evolved over the last century to incorporate a large tool kit of tactics based on mechanical control, chemical control, biological control, and the current practice of integrated vector management (IVM) (Benelli et al. 2016, Patterson 2016). Integrated vector management is designed to combine tactics from all 3 areas while emphasizing the importance of understanding the biology and ecology of mosquitoes, including differences between life stages and species type (Beier et al. 2008, Lizzi et al. 2014, Marcos-Marcos et al. 2018). Mosquito control is typically managed at the local administrative level. Thus, with the large scope of available control strategies paired with area-dependent environmental, disease, and financial priorities, vector control interventions and implementation strategies vary widely between regional agencies (NACCHO 2017). The operational timing and initiation of control is commonly based on action thresholds. An action threshold, as defined by the United States Environmental Protection Agency (EPA) under the broad umbrella of integrated pest management (IPM), is “a point at which pest populations or environmental conditions indicate that pest control action must be taken” (EPA 2022). Integrated pest management is primarily guided by economic threats of agriculture, whereas IVM is specifically motivated by risks of nuisance or transmission of vector-borne disease (VBD) from arthropod vectors (Lizzi et al. 2014) such as mosquitoes, ticks, sandflies, midges, and fleas (EFSA 2018). In this article, VBD specifically refers to those diseases vectored by mosquitoes. Vector-borne diseases are a global public health concern, since mosquitoes vector a multitude of arboviral and parasitic agents worldwide (Dahmana and Mediannikov 2020). It is estimated 2.5–9.3% of the ∼3,500 mosquito species are associated with a human disease, with 76% of all known and/or potential vectors belonging to the three genera of Anopheles, Culex, and Aedes (Yee et al. 2022). It can be argued that action thresholds for mosquito control are also a response to economic threats, since VBDs lead to financial and overall societal costs to the detriment of households (lost wages, treatment costs) to nations (reduced labor force, slower economic growth, reduced GDP) (WHO 2017). Historically to the current day, vector control remains the primary way to prevent the spread of VBDs, especially for those neglected tropical diseases that lack an accessible vaccine or effective clinical treatment (Wilson et al. 2020).

Successful IVM programs rely on active surveillance of mosquito populations to inform control decisions (AMCA 2017). However, alongside the differences in chosen control strategies, action thresholds and overall surveillance activities are rarely standardized between local agencies due to variation in both operational capabilities and public health priorities (Aryaprema et al. 2022). Action thresholds are triggered by exceedance of predetermined values assessed via passive or active surveillance activities. Assessed surveillance factors most commonly include human or animal disease case reports, presence or abundance of adult mosquitoes, mosquito larval indices, mosquito landing rate counts, arbovirus detection, sentinel animal seroconversion, and service requests from community members (AMCA 2017). The thresholds enforced by mosquito control programs to initiate control can include district level decisions on top of (but still compliant with) codified legislative justifications to account for the nuances of control methods enacted by different programs. This is in addition to the nuanced necessity for implementing a certain method versus others which requires weighing decisions based on information such as the seasonality, perceived issue (e.g., excess pest biting, VBD case reports), and extent of geographic coverage required (e.g., aerial spraying over several subdivisions versus ground ULV spraying on one street). Most action thresholds are inherently reactive and are more than likely triggered timeline-wise on a date that temporally surpasses the temporal period when initiation of control activities would have more effectively prevented the start of nuisance issues or disease transmission.

In a recent survey of global mosquito control programs, some programs in the United States, Asia, Africa, and Australia reported that they do not have a designated action threshold or follow any regional or military guidelines, even if they collect the necessary basic surveillance data to develop thresholds (Aryaprema et al. 2022). International guidelines do exist with action thresholds for yellow fever and dengue that rely on Stegomyia indices (WHO 1971, PAHO 1994, WHO 2009). Numerous studies have utilized these thresholds as a target for operational control or a baseline to compare study results (Aryaprema et al. 2023), but the overall generalizability and effectiveness of such indices is under question (Focks et al. 2004, Bowman et al. 2014). There is a lack of standardized guidelines for programs to develop their own action thresholds with or without utilization of historical surveillance data. A recent literature review of action thresholds found a limited number of studies from the last decade that used empirical methods to develop an objective threshold value or a threshold model (Aryaprema et al. 2023). Those threshold studies, although focused primarily on control of Aedes species and dengue, described numerous methods to develop a proactive operational threshold using common surveillance parameters. Importantly, most studies incorporated mosquito surveillance counts into development and/or the final threshold index value. Only a few instead created thresholds based on local climate measurements such as temperature, rainfall, and relative humidity. Previous literature has made evident that these weather variables can impact the development, behavior, reproduction, and mortality of mosquitoes (Gage et al. 2008), albeit to various degrees dependent on species, study location, study design, etc. (Aryaprema et al. 2023).

Overall, there is an increasing need for proactive and evidence-based action thresholds that are sustainable and realistic to implement at local levels. It is also important to create more generalizable action thresholds that do not rely on the typical surveillance of mosquito counts or disease cases, especially for entities that do not have the resources for sustainable surveillance. Action thresholds based on local climate measurements would be a good solution to address this issue. The primary motivator for this study is that many pest managers deployed in US military operation zones abroad lack both sufficient vector control guidelines and routine mosquito surveillance systems. In the event of a pest issue, the Armed Forces Pest Management Board (AFPMB) provides support to global military installations through deployment of a limited number of specialized personnel. Integrated pest management technical guidelines exist (AFPMB 2013) but do not specify direct action thresholds to trigger proactive control activities. For an already problematic situation, there will be a delay in arrival of specialized personnel, assessment through temporary surveillance activities, and initiation of appropriate control measures. Therefore, there is crucial time lost and risk accrued when an issue does arise because of the reactive response protocol. An optimal action threshold would entail 1) passively collected input parameters that do not require resource-heavy data collection protocols, 2) generalizable values that account for different ecological niches across genera and environment type, and 3) alert for upcoming rather than current nuisance/disease transmission issues with several weeks of lead-time.

The Anastasia Mosquito Control District (AMCD) in St. Johns County (SJC), Florida, utilizes an IVM approach with routine mosquito and arbovirus surveillance activities followed by a combination of adult and larval control techniques. Adult control by the AMCD is prioritized following one of the few legislative action thresholds set forth by the rule 5E-13.036 of the Florida Administrative Code and Florida Administrative Register that states the minimum threshold is met if “adult mosquito populations build to levels exceeding 25 mosquitoes per trap-night or 5 mosquitoes per trap hour during crepuscular periods” (Florida Department of State 2006). AMCD’s threshold for larval control is ≥ 1 larva per dip into a water body. In addition to weekly mosquito trap threshold, other surveillance parameters that factor into operation control decisions include a positive sentinel chicken seroconversion (West Nile virus, St. Louis encephalitis virus, and eastern equine encephalitis virus), service requests, and the knowledge and intuition of operator personnel. Data from all surveillance programs and operational activities have been recorded and digitally archived since 2004.

The main objectives of this study were to 1) explore the analytic capacity of AMCD historical mosquito surveillance records, 2) determine key variables for predicting optimal initiation of mosquito control using historical records in conjunction with public climate data, and 3) optimize proactive and generalizable climate -based action threshold for controlling mosquito genera of public health concern. The results published here significantly address the first two objectives while yielding guidance for understanding the needs and limitations of accomplishing the third objective, which was unsuccessful for reasons discussed here.

The two trap systems for mosquito surveillance at AMCD are 1) BG traps with BG lure and CO₂ (dry ice) attractants (BG) and 2) CDC light traps with octenol lure stick (CDC LT). Centers for Disease Control and Prevention LT data has been collected since 2004 while BGs have been in use since 2014. Both trap systems are a weekly responsibility during the active surveillance season. Traps are set out in designated locations on a consistent weekday morning and then picked up the following morning, with approximately 24 h of total field deployment. BGs are placed only in the coastal downtown city of St. Augustine for approximately 48 wk of the year; meanwhile, CDC LT are set out throughout the entire county for a shorter active season, March/May–December (Fig. 1). After pickup, the trap collections are transported to AMCD for identification to species level using appropriate taxonomic keys (Darsie and Ward 2005). Total species counts (female only) for each individual trap site are recorded both in an Excel surveillance record template and AMCD’s online database.

For this study, the AMCD database was mined and curated to best format the historical records for analysis using a mixed bag of software tools including Excel (Microsoft; Redmond, WA, USA), ArcMap v.10.8 (Esri; Redlands, CA, USA), and R v 4.1.0 (R Foundation for Statistical Computing, Vienna, Austria) alongside manual examination of weekly records to find events of trap inactivity. As seen in the original surveillance record sheets, each year commonly had instances that a trap site was clearly not an accurate measure of mosquito counts. These instances were due to either trap dysfunction (dead battery, full of water, stolen) or lack of deployment (holiday, bad weather conditions). These instances were important to find and mark as data points with null values “NA” instead of a value of “0” to avoid skewing of statistical analyses as much as possible. A simple summary of all trap locations, surveillance season duration, and geographic coordinates was updated for each year for future analyses using these historical data. The number and geographic distribution of CDC LT varied during the sixteen-year timespan (2004–2019) with the minimum of active seasonal traps in 2015 (n = 30) and a maximum in 2006–2007 (n = 52) (Fig. 1). Comparatively, BG locations remained consistent from 2014–2019, aside from the increase in trap locations between 2014 and 2016 (n = 6) to 2017–2019 (n = 12).

A timeline of epidemiological weeks (epiweek) was developed to match collection days across trap systems and years. Each epiweek was set as Monday–Sunday to match the start of AMCD’s workweek and operational surveillance and control activities. Daily readings of total rainfall, minimum temperature, average temperature, and maximum temperature were retrieved via the Climate Data Online Search tool (National Oceanic and Atmospheric Administration) from a weather station located in a township of St. Johns County, FL (Hastings 4 NE, FL US GHCND: USC00083874; 29.7652°, −81.4697°) (Fig. 1). Daily relative humidity data were retrieved from the Florida Automated Weather Network (FAWN) (University of Florida Institute of Food and Agriculture Sciences Extension), which sourced from a nearby station in the same township (FAWN ID: 270; 29.69332°, 81.44485°). Daily climate values were averaged (temperature, relative humidity) or summed (rainfall) across each epiweek to calculate the final climate parameters of minimum temperature (tmin), maximum temperature (tmax), average temperature (tave), total rainfall (prcp), and relative humidity (rh) per epiweek.

Out of approximately 40 mosquito species present in St. Johns County, 8 were chosen as the focus of modeling efforts due to their higher abundance in the county and known/suspected status as vectors of infectious human disease agents including West Nile virus, Chikungunya virus, Zika virus, dengue virus, Yellow Fever virus, Mayaro virus, eastern equine encephalitis virus, Japanese encephalitis virus, La Crosse encephalitis virus, St. Louis encephalitis virus, malaria Plasmodium, and Keystone virus. Species counts of Aedes taeniorhynchus (Wiedemann), Ae. sollicitans (Walker), Ae. atlanticus (Dyar and Knab), Ae. infirmatus Dyar and Knab, Culex nigripalpus Theobald, Cx quinquefasciatus Say, and Anopheles crucians were taken from CDC LT records, whereas Ae. aegypti (L.) and Ae. albopictus (Skuse) were taken from BG trap records. This inclusion of species records from CDC LT versus BG was based on the differential attraction of those species to each trapping system with preference given to CDC LT due to the larger sampling size. Species were grouped in separate models according to genera and basic ecological habitat, hereby referred to as “species group” (Table 1). Another category named Total Mosquitoes was created with the total number of all mosquitoes collected, disregarding species level, at each CDC LT site for each active surveillance week. The value of “mosquito trap-night average” was calculated for each individual epiweek across all years using the total number of mosquitoes (by species group) divided by the total number of active (and presumed functional) traps deployed that week.

The corresponding β coefficients for final models include the intercept (β0), epiweek by individual week (β_{ith epiweek}), and β of significant climate variables. Only the climate β coefficients are multiplicative in the model equation. The parameter β0 is static and β_{ith epiweek} is dynamic (1 coefficient per epiweek 1–51), but neither are affected or effect the inclusion of novel climate data on output model predictions. Backwards elimination entailed the initial inclusion of all the top lagged climate variables (see equation for potential final model structure) and one-by-one removal of nonsignificant variables until only significant parameters remained.

Of note, 3 years of mosquito collection records were excluded when earliest iterations of Models #1 failed to be generated, albeit with a different statistical analysis method (ANCOVA). Exploratory data analysis revealed that the years 2006, 2010, and 2017 were outliers of annual precipitation and were noticeably drier (2006, 2010) or wetter (2017) years compared to the rest of the 2004–2019 period. Final training datasets for CDC LT species groups covered the years 2004–2019 (except the three outlier years), whereas BG covered all 2014–2019, including 2017, due to the already reduced sample size of years.

Models #1 results prompted a more thorough investigation of outlier points in 2004–2021 mosquito and climate datasets. We distinguished weekly outliers by year and conducted exploratory data analyses using comparison with climate. Results from this analysis guided decisions for Models #2.

The second iteration models were based off results from the first using a new paradigm. The same timeline for training and testing data was used, except the initial 3 excluded years were re-included, and epiweeks of influential (not outlier) trap-night averages were removed according to species model.

Climate parameters again included tmin, tave, tmax, prcp, and rh. Index threshold values were chosen for the number of “hot,” “wet,” and “humid” days based on exploratory data and outlier analyses for each parameter. The rainfall parameter prcp had one potential threshold value to test (> 0 inches) as “wet,” while three threshold degree measures were separately tested as “hot” for tmin (60°F, 65°F, 70°F), tave (70°F, 75°F, 80°F), tmax (80°F, 85°F, 90°F), and as “humid” for rh (75%, 80%, 85%). A binary classification was ordered each day as over “1” or under “0” the threshold value of a certain parameter. The total number of days over threshold was summed for each epiweek and combination of epiweeks. To retain a lag period and additionally investigate the statistical power of a variable spanning more than one week, fifteen separate temporal variables for each climate parameter were tested (Fig. 3). These temporal variables included the total of each individual epiweek individually of 1 to 6 wk previously (n = 6), and combination lag totals (n = 9) of 2 to 4 epiweeks starting from the epiweek lagged 1, 2, or 3 wk previously. Spearman’s correlation was used to distinguish the temporal variable of greatest correlation coefficient (ρ) magnitude at the individual lagged epiweek and each 2, 3, 4 combination lag epiweek for all three potential thresholds of tmin, tave, tmax, and rh. Spearman was used in favor of Pearson’s correlation for this iteration due to rank ordering of ordinal climate values (i.e., number of days). The optimal threshold index value for each temperature and relative humidity parameter moved forward to modeling, alongside the top 4 temporal variables of prcp with its 1 index threshold.

By species group, univariable models were built for each top temporal variable (n = 16). The most optimal univariable climate model for each species group (n = 1-3) was isolated using criteria of Akaike information criterion (AIC) and significant p-value alongside minimization of Pearson Chi Square/Degrees of Freedom (DF). The optimal univariate parameters were then combined and tested as a multivariate climate model until only significant model parameters remained.

The same process of examination for testing and validating the first iteration of models was used for the second iteration.

The selected climate variables differed in lag time, significance, and correlation strength between the species groups (Fig. 4). Four of the species groups significantly correlated to at least one lag variable for all 5 climate parameters, while relative humidity did not reach significance for either Saltmarsh or Container groups. Overall, prcp at a 2-wk lag (prcp2) was the most consistent climate variable and was selected for inclusion in 5 of 6 species group models (See Fig. 2 to distinguish lagged climate variable names).

Final parameters in species group models were reduced to 1—2 climate variables, alongside the temporal class variable of epiweek (Table 2). Predictors and estimated effect differed between species group GLMs. The variable prcp2 was the most common final predictor in five of six species group models (excluding Swamp), and it is the only significant climate variable for two models (Standing, Total Mosquitoes). Besides prcp2, a single temperature parameter remained significant in the final GLMs for Container (tmin1), Floodwater (tmax3), and Saltmarsh (tmin3). Swamp stood out because rh1 was the only significant climate parameter, whereas rh failed to make the final cut in the other GLMs. The parameter tave was eliminated from all models due to its extreme collinearity with other temperature parameters. The β coefficients revealed that temperature parameters had a reductive effect, while prcp and rh effects were instead augmentative to the predicted trap-night counts. The epiweek predictor had a significant effect in all models; however, individual β coefficients by epiweek were only significant for a portion of the season for Total Mosquitoes (wk 21–43) and Container (wk 17–42), or no weeks at all for the other species groups. A least-squared means (LSMEANS) estimate of trap-night averages was calculated for each species group across each epiweek in a typical surveillance season (CDC LT: 10-51; BG: 1-51) (Fig. 5). The estimates reflected the observed frequency of epiweek events in 2004–2019 (CDC LT) and 2014–19 (BG), and thus future likelihood, that the trap-night average of a certain species group surpassed either 25 mosquitoes (Total Mosquitoes: n = 180, Floodwater: n = 80, Container: n = 25, Swamp: n = 23, Standing: n = 9, Saltmarsh: n = 3) or 8.33 mosquitoes (Total Mosquitoes: n = 340, Floodwater: n = 158, Container: n = 85, Swamp: n = 152, Standing n = 37, Saltmarsh n = 13).

Calibration (Table 2) was based on available diagnostics on model fit and residual diagnostics from the GLIMMIX procedure (plots=residual panel). One avenue of model development included attempts to pair variables from the same lag time. Half of the final GLMs (Floodwater, Saltmarsh, and Container) had two climate parameters from different lagged weeks. The models were rerun after switching out the lag time of one climate parameter (prcp) to match the lag week of the other climate parameter (tmin/tmax), and vice versa, and gauged model fit criterion alongside the significance of the replacing variables. No attempts to match lag weeks succeeded in minimizing AIC or Pearson Chi-Square/DF.

Diagnostics of residuals of the training GLMs with quantile-quantile plots and histograms revealed that error was relatively normally distributed, but there was heteroscedasticity and increased spread as observed mosquito trap-night averages increased. Outliers were examined but ultimately not removed from training datasets (explained below).

Two GLMs (Total Mosquitoes, Floodwater) were run with training GLM parameters on climate testing data accurately reflected significance of all model predictors. However, the prcp2 β coefficients were dissimilar and exceeded the training 95% CI (Table 3). Meanwhile, the prcp2 β coefficients for Saltmarsh and Container were within the training 95% CI, but the predictors themselves did not reach significance on testing data. The parameter tmax3 for the Floodwater GLM was the only variable to remain significant and land in the CI.

Testing GLMs for four species groups did not significantly fit with any or all training GLM parameters. Epiweek (but not climate) was significant for Container, climate (but not epiweek) was significant for Standing, and neither climate nor epiweek remained significant for Swamp and Saltmarsh. The observation that epiweek alone did not reach significance for 2020–2021 Swamp, Saltmarsh, or Standing data was unexpected. The RMSE paralleled the problematic variation of residuals in training model calibration but was much lower in the models that failed to be significant in validation. Total Mosquitoes and the most numerous species group, Floodwater, had the largest error. Standing, Swamp, then Saltmarsh had the least error. Residual plots for all models had similar behavior of residuals; normal distributions, but heteroscedasticity and greater spread as observed mosquito trap-night averages increased.

The sensitivity, specificity, accuracy, and kappa statistic differed between each species group and the established (25) and trap-subjective (8.33) mosquito threshold (Table 4). The kappa statistic landed in the range of 0.20–0.40 between models, meaning the models achieve a fair agreement with the ground truth that is 20–40% better than chance alone (Landis and Koch 1977).

The Total Mosquito model with the prcp2 variable displayed the highest sensitivity, but the lowest specificity, of all species group models with both 25 and 8.33 mosquito thresholds. Meanwhile, the models Floodwater and Container were more specific than sensitive at the 25-mosquito threshold, but more sensitive than specific at the 8.33-mosquito threshold. Neither Standing nor Saltmarsh surpassed either threshold during 2020 or 2021, which was correctly predicted by models of both species groups with climate testing data. Thus, both threshold values for the two species groups were completely accurate and specific but sensitivity could not be evaluated. Contrastingly, Swamp also did not have observations surpassing the mosquito threshold, but its model incorrectly predicted thirteen weeks to surpass the 8.33-threshold which resulted in a lowered specificity and overall accuracy.

The selected climate variables slightly differed in optimal lag and combination lag variables, threshold temperature indices, and correlation strength between the species groups (Table 5). Excluding many Saltmarsh or Container versus rh variables, most Spearman correlations with climate variables during initial parameter selection were significant at p < 0.0001 including all top combination variables listed in the table. All correlation relationships were positive (See Fig. 3 to distinguish lagged climate variable names).

Univariate GLM models with climate variables were tested and the variables at optimal threshold indices and lag times were combined into multivariate models (Table 2). In this first exclusion stage, humid failed to reach any significance with Saltmarsh or Container and was not included in those multivariable model attempts. Total Mosquitoes and Standing were the only two groups with all three climate parameters (hot, wet, humid) included in the final GLM equation. Floodwater, Container, and Swamp each included only two climate parameters in the multivariate model, while Saltmarsh remained the singular univariate climate model. Overall, the variables of wet and hot were included in five final species group models with humid only in three. Temperature threshold values designating days as hot varied between the four species groups; only Standing and Floodwater shared the same optimal definition of “hot” as exceeding daily tmax 80°F. Lag times of hot were also inconsistent; however, hot thresholds based on tave (hot1, hot123) relied on more recent time lags than those based on tmax (hot34, hot345). All wet thresholds were the same from the start (0 inches), and all final lag combinations (wet23, wet234, wet2345) were very similar in timing. The optimal humid variable for the three including models has a consistent threshold index value of 80% rh and consistent timing (humid12, humid1234).

Epiweek was significant in all models. β coefficients by individual week were again significant only for Total Mosquitoes (wk 16-45) and Container (wk 14-45), a slightly broader range compared to Models #1. The candidate GLMs for Models #2 for each species group (n = 2–5) were culled to one optimal model using calibration diagnostics. Final models and climate parameters are highlighted in Table 2.

Diagnostics of residuals of the training GLMs with quantile-quantile plots and histograms revealed that error was relatively normally distributed, but there was heteroscedasticity and increased spread as observed mosquito trap-night averages increased. The residuals were more densely concentrated in the bottom half, implying that the model tended to overestimate the average mosquito trap count.

A LSMEANS estimate of trap-night averages was again calculated across each epiweek in a typical surveillance season (CDC LT: 10–51; BG: 1–51) (Fig. 6). The LSMEANS estimates between the Models #1 and Models #2 iterations had similar patterns; Saltmarsh and Standing exceeded neither 25- nor 8.33-mosquito thresholds, Swamp exceeded only 8.33, while Total Mosquitoes exceeded both thresholds. Contrastingly, Container and Floodwater LSMEANS did not surpass the 25-mosquito threshold this time. The LSMEANS estimates were also examined for the other initial models (not shown). Notably, GLMs built on the number of hot days versus the number of wet days varied within a species group. For example, Models #2 candidate models of Total Mosquitoes and Floodwater involving only temperature (tmin, tave, tmax) often had higher weekly estimates and more weeks exceeding the 25-mosquito threshold compared to models with only prcp. Meanwhile, models with only prcp had estimates exceeding the 8.33-mosquito threshold earlier in the surveillance season. Those built on only rh had similar early seasonal exceedances of 8.33-mosquito threshold like prcp models, however with greater estimates that more often exceeded values of 8.33 or 25, dependent on species.

All GLMs built on testing data lacked the significance of one or more training predictors, and the majority of β coefficients were dissimilar and exceeded the training 95% CI (Table 3). The prcp predictor (wet) remained significant in four of five testing GLMs, rh predictor (humid) remained significant in one of three, while none of the temperature (hot) variables remained significant. Like the Models #1 iterations, the epiweek again remained significant in Total Mosquitoes, Floodwater, and Container and not significant with Standing. No GLM could be built using Saltmarsh or Swamp testing data. β coefficients of humid for Total Mosquitoes (humid1234) and wet for Container (wet23) were the only ones within the 95% CI while all other β coefficients exceeded the 95% CI from training GLMs.

The RMSE was highest in Total Mosquitoes, approximately double the next highest error measure seen with Floodwater. The order of species group by RMSE matched that of Models #1 iterations, except for a switch between Container and Swamp, with Saltmarsh and Standing again having the lowest error overall. Diagnostics of residual plots matched those in testing of Models #1, with increasing but equally distributed spread as observed mosquito trap-night averages increased.

All sensitivity, specificity, accuracy, and kappa statistics results were the same or slightly improved compared to Models #1 iterations, aside from the decreased specificity and accuracy of the Swamp and Standing GLM at the 8.33 mosquito threshold (Table 4). The kappa statistics landed in the range of 0.27–0.49 between models, meaning the models achieved a fair to moderate agreement with the ground truth (Landis and Koch 1977).

There were outliers in the testing and training mosquito datasets that likely influenced the model accuracy; however, removing natural ecological phenomena to smooth the datasets was considered contradictory to the aims of the study. In mosquito control programs, population peaks are vital knowledge for operational personnel to prepare for and manage as best possible. A secondary exploratory data analysis on outliers’ post-completion of Model #1 revealed patterns to follow up on. Annual climate and annual mosquito collections were also tested against year in case of temporal patterns.

Each season had frequent outliers in mosquito trap-night averages (in context of statistical distribution of values in timeframe) with 0–6 outliers dependent on year and species. The Floodwater group had on average the most outliers per season (3.72, n = 68, followed by Saltmarsh (3.11, n = 56), Container (2.88, n = 23), Standing (2.39, n = 43), Total Mosquitoes (2.11, n = 38), and lastly Swamp (1.67, n = 30). Species models differed in the average, minimum, and maximum values of trap-night averages that were identified as an outlier during the training and testing seasons (Table 6). Species groups commonly displayed repeating outlier events over 2–4 consecutive weeks including Floodwater (n = 15 years), Saltmarsh (n = 7 years), Standing (n = 8 years), Container (n = 3 years), but less so Swamp (n = 2 years). These consecutive outliers were not always simultaneous across species; Saltmarsh, Floodwater, and Standing sometimes showed distinct outlier epiweek clumps. The seasonal timing of outlier events differed percentage wise by species group (Table 7), with the majority (50–71%) of outliers for each species group occurring in the 3rd quarter of the year (July–Sept; epiweeks 27–39).

Timing of outlier climate events is likely important to the magnitude of effect imparted on mosquito populations, particularly the synergy of temperature and rainfall. There was no clear consistent order of species group outliers during a surveillance season; however, most outlier rainfall events did lead to an outlier of one or more mosquito species groups 1–4 weeks later. Seasonal temperature showed interesting changes over the 16-year time span. The variables tmax and tmin had outliers that only started appearing in 2014 and continued to pop up 1–3 times per following season. Of note, all these outliers were at temperatures lower than the normal distribution. In addition, correlation of annual climate trends versus year (2004–2021) revealed a significant positive relationship with the tmax, tmin, and tave of average seasonal temperature (Table 8) with annual temperatures increasing approximately 1-3°F (with intra-annual cycling) over the study period. The parameters for prcp and rh did not show a distinct positive or negative trend. However, prcp outliers were common during the study period with 1-5 outliers in 10 of 13 training years (average 1.72 outliers). The driest year, 2006, had the most outliers (n = 5) but these outliers occurred mostly during the epiweeks peripheral to the active season (6, 8, 24, 51, 52). The wettest season, 2017, had the next most prcp outliers (n = 4) but these outliers occurred closer to the middle of the active season (24, 37, 39, 40). Prcp outliers were typically smaller in dry versus wetter years. Approximately 10% of prcp outliers occurred in the 1st quarter (Jan–March), 32% in the 2nd (April–June), 39% in the 3rd (July–Sept), 19% in the 4th (October–December). There was no correlation between the trap-night averages of species groups and climate outliers by year.

The value of outlier trap-night averages in 2004–2021 of each species group was broken down into percentages of their individual species (Table 9). Proportions between combined species for species groups Floodwater (Ae. atlanticus 84.21%, Ae. infirmatus 84.21%), Saltmarsh (Ae. taeniorhynchus 93.43%, Ae. sollicitans 6.57%), Standing (Cx. nigripalpus 87.08%, Cx. quinquefasciatus 12.92%), and Container (Ae. aegypti 39.05%, Ae. albopictus 60.95%) closely reflect the percentages by trap type (see Table 1), aside from the more closely matched Ae. aegypti and Ae. albopictus. Unsurprisingly, the dominant species in each species group by annual quantity was also the dominant driver of outliers.

Correlation tests for annual climate versus annual mosquito trap-night averages were run for individual species. Total Mosquitoes (r = 0.7718, p = 0.0002), Ae. atlanticus (r = 0.6903, P = 0.0015), Ae. infirmatus (r = 0.7185, P = 0.0008), Cx. nigripalpus (r = 0.5977, P = 0.008), and An. crucians (r = 0.6402, P = 0.0042), had a positive correlation with total annual prcp. Aedes aegypti had a positive relationship to tmax (r = 0.7447, P = 0.0340), while the other container species Ae. albopictus was negatively correlated to all three temperature parameters: tmax (r = −0.8521, P = 0.0072), tmin (r = −0.8521, P = 0.0072), and tave (r = −0.9019, P = 0.0022). No correlations were found with the rh parameter with Saltmarsh species Ae. taeniorhynchus or Ae. sollicitans. Additionally, both Ae. albopictus (r = −0.7767, P = 0.0234) and, though barely significant, Cx. quinquefasciatus (r = −0.4772, P = 0.0452) had a negative relationship with year itself.

The current literature includes numerous statistical analyses conducted to understand the primary factors driving the dynamics of mosquito species’ populations and VBDs’ transmission. Fewer studies have constructed empirical action thresholds with a direct value or categorical risk model of mosquitoes, climate, or disease cases that can optimally trigger mosquito control activities (Aryaprema et al. 2023). This here is the first known analysis to use historical control program data to create dynamic lagged climate models for habitat genera groups and to validate model effectiveness using current operational thresholds.

Climate predictors have been historically associated with dynamics of mosquito populations and disease transmission. Importantly, we found in this study that rainfall, temperature, relative humidity varied in significance in estimating average trap-night counts across separate species groups. Rainfall was the leading predictor across our models as determined by the number of inclusions in final model equations. Temperature was a close 2nd in significant impact while relative humidity showed none to minimal impact in most species groups, with clear contrast to the 2 final models of Swamp, i.e., An. crucians. Observing such a range in temperature predictors cannot be unexpected because mosquito species typically inhabit specialized ecological niches and can therefore be differentially impacted by daily and seasonal climate patterns. This observation was supported by tracking population dynamics of species groups and individual species. The abundance levels between species of a group were clearly unproportionate, and thus models based off the combination total of differentially active species were of initial concern. A tangential analysis of seasonal time points of reaching 5%, 50%, and 95% population numbers within a season (not shown) revealed close dynamic trends between species combined to a species group, and we felt the decision to model this way was justified. Species that purportedly share the same ecological habitat can be expected to react similarly to routine environmental cycles or challenges. Biological differences may also account for variation in population dynamics, as seasonal weather can theoretically have delayed effects on the populations of following years through overwintering behaviors and egg desiccation.

This study also demonstrated that lag time was critical to the significance of each climate parameter and that climate threshold indices can be optimized for separate species groups. Previous threshold models and generated climate thresholds rely on equivocal climate measurements (Aryaprema et al., 2023). Studies have differentially incorporated temperature in the form of maximum, minimum, average, or degree days, or compiled rainfall data into total or averaged values. In this analysis we expected the correlation between the three temperature indices, but also saw a pronounced difference in the range of daily and seasonal variation between minimum and maximum temperature indices. Presumably, it would be preferential for a certain ecological zone to rely on the index that provides the most information about changing climate. A more distinguished variation would have a better chance of predicting variation in the dependent variable, i.e., mosquito abundance. Our exploratory analysis of climate baseline and outliers led to decisions to modify the format of climate parameters. The Model #2 iteration was a specific attempt to understand if the number of rainy days was a more accurate predictor than the total amount of rain. We additionally translated different threshold values of temperature and humidity that correspond to approximate physiological ranges associated with mosquitoes to define optimal hot and humid days. The lag time was stretched further to potentially capture longer delayed effects, and it was observed that impact of minimum temperature and average temperature had offset timings. Still, the predictors between the 2 model iterations were closely matched, particularly with the significant impact and timing of rainfall, and the no to low impact of humidity on all groups except Swamp.

It is worthwhile to focus on the differences in those models with iterations containing both rainfall and temperature variables. The effect of average weekly temperature matched or exceeded the effect of total weekly precipitation in Models #1; meanwhile the number of wet days had a far greater impact than the number of hot days in Models #2. Relative humidity was significant for more species groups in Models #2 than in Models #1 although with minimal effect. Overall, it is difficult to objectively state which model iteration performed better. Only moderate accuracy was achieved during validation tests by either iteration, albeit slightly higher in Models #2, and the residual error steadily increased alongside the observed trap-night averages. These models are not intended to predict an exact count of mosquitoes but to predict whether that mosquito count will exceed a set threshold value. There was a clear tradeoff between sensitivity and specificity using the two validation thresholds tested. An 8.33-mosquitoes-per-trap-night average was more sensitive and good at correctly predicting when the observed collections would exceed threshold (true positives), while the 25-mosquitoespertrap-night average had better specificity and was good at predicting when the observed collection numbers did not exceed threshold (rue negative). A subsequent step would be to create receiver-operating characteristic curves (ROC) to optimize the threshold for model performance and strike the balance between sensitivity and specificity. We did not perform an ROC analysis in this study due to the dominant limitations of ongoing operational control activities that could not be properly addressed in the models.

The third takeaway is that seasonality (aka epiweek) was an important modifier to the effect of climate. The significance of the epiweek β estimate did change from week 1 to week 51 and was more reliably significant in the middle of the season, aka the summer months, and only for the most abundant species. This is likely an artifact of the limitations of AMCD’s data collection timeline limitations and the unproportionate abundance of mosquitoes in different species groups. Exploratory data analysis of temporal species dynamics revealed jagged population curves for each species and all surveillance seasons. Many of these peaks were distinguished as outliers, likely due to favorable climate conditions or weather events preceding the population spike. In the Model #1 iterations we did not exclude outlier values of mosquito counts because of the assumed inherent tie to climate predictors that were the focus of our analysis. These models had generally poor predictive capability particularly when trap-night averages surpassed 10–20 mosquitoes. There are no standardized guidelines to dealing with outlier values of mosquito count data. We ultimately removed only influential data points (extreme outliers) instead of all statistically “unusual” trap-night averages. Mosquito model publications may state the presence of outliers but might fail to clearly mention what happens to such data points whether modified, excluded, or included in analysis (Lindahl et al. 2012, Bravo-Barriga 2017, Reiskind et al. 2017, Giordano et al. 2021). Some studies do specifically state how outliers and/or extreme outliers were defined and separately ran models with and without (Okanga et al. 2013, Karki et al. 2016) with Karki et al (2016) noting a more defined trend line and statistical correlation between trap type collections using datasets that excluded the outliers. We performed an in-depth exploratory analysis on the occurrence of the statistical outliers of mosquito counts in relation to outlier values of weekly climate measurements. Individual species varied by the quantity and frequency of outlier events in weekly trap averages during the timeline. The sudden population reductions in weeks following outlier population booms are most likely due to targeted mosquito abatement activities in prioritized zones. It is likely that if the trap site areas were left uncontrolled, population trends would be far more smoothed or continuing exponential curves and would encompass our outlier points into the normal distribution of seasonal trap counts. The type and lag time of important predictors are further impacted by the degree to which a species is impacted by the activity of other species, such as from intraspecies competition or action threshold exceedance in the same operational zone. This idea was discussed by Steck et al. (2021) based on a comparison of seasonal population trends of Ae. atlanticus versus Ae. infirmatus in SJC. Control operations intended for one species, especially ultra-low volume fogging missions, theoretically affect unproblematic species before they reach their out peaks. The low values of outlier events in Standing, Saltmarsh, and Swamp are likely the result of premeditated management. Overall, it is taken as a sign of the effectiveness of control operations that these species do not experience outlier events or reach greater events of mosquito counts. In SJC, Floodwater species, specifically Ae. atlanticus, are widespread throughout the county and frequently exceed the 25-mosquitoes per trap-night average action threshold value and are thus a leading trigger species for control. Anopheles crucians can reach similar abundance levels sufficient to trigger a fog mission, but contrastingly, high trap counts go largely ignored and do not lead to a large-scale adulticide application.

Another related question is how the amount of total rainfall potentially switches from augmentative to reductive toward mosquito populations, and the delayed effect of extreme and thus outlier weather events. For example, Hurricane Irma hit SJC on September 12, 2017 (epiweek 37). Trap collections the following 3 weeks were 37,51, 5,804, 228, then 15,589, respectively. Meanwhile, epiweek 33 had the 2nd highest collection of the season at 13,002 with no major rainfall events but consistently rainy weeks in the previous 6 wk (0.56–3.26 inches). An overabundance of rain has been shown to lead to a flushing out effect of larval habitats and reductions in disease transmission (Paaijmans et al. 2007, Dieng et al. 2012, Seidahmed and Eltahir 2016, Benedum et al. 2018). Another Florida mosquito control district, overseeing Collier County in the southwest of the state, associated a significant population crash in their seasonal Ae. taeniorhynchus trap counts with the impact of Hurricane Irma the season before (Lucas et al. 2019). Extreme weather events like hurricanes, monsoons, flooding, or simply above average rainfall may also result in novel or expanded mosquito habitats that are associated with immediate or delayed disease results (Barrera et al. 2019, Caillouët and Robertson 2020, Coalson et al. 2021).

It should be noted that there were numerous limitations to our analysis. The methodology for model development here incorporated the testing of various iterations to resolve potential issues of masked or overexposed statistical relationships due to the fickleness of real-time ecological data. Statistical power issues were encountered as is expected in IPM datasets that are guided by operational and logistic decisions. The number and geographic placement of AMCD trap sites, specifically CDC LT, greatly differed across surveillance seasons. The use of county-wide trap-night averages for mosquito counts reduced the spatial resolution and is a major limitation of the analysis. We could not account for either spatial autocorrelation or the climatic and environmental variation within the county. A related sub-limitation is any temporal or spatial population changes in individual mosquito species at certain trap sites. St. Johns County encompasses a 608-square-mile area between the St. Johns River (approximately 14–19 miles inland) and the Atlantic coast, with a parallel intercoastal waterway. Climate parameters measured at one inland station are not the best representation of all trap sites in the county, likely only representing general seasonal trends, and certainly do not capture microhabitat conditions. Southward in Palm Beach County, FL, temperature and relative humidity significantly differed with distances 0–15 km inland from the Atlantic coast in a study of the changing population dynamics of Ae. aegypti versus Ae. albopictus (Hopperstad and Reiskind 2016). Logistically, any control program would rely on one or very few weather stations to indicate overall trends, even if at a coarse resolution. Temperature, rainfall, and relative humidity are fundamental through the literature covering nuisance and disease vector species, particularly from the perspective of understanding consequences due to climate change (Reiter 2001, Caldwell et al. 2021). This analysis showed statistical evidence of a gradual increase in annual SJC temperature indices over the 18-year timespan (2004–2021). The environment including land use land cover (LULC) and LULC change is also not accounted for. According to land-cover proportions calculated with National Land Cover Database data (Multi-Resolution Land Characteristics Consortium [MRLC]) approximately one-third (30%) of the land-water area of SJC has changed by LULC category sometime between 2004 and 2019 (MRLC 2021). There were noticeable increases in area with developed land (open space, low intensity, medium intensity, high intensity) (22–149%) and decreases in shrub/scrub (47%), hay/pasture (31%), and herbaceous (32%). A previous analysis using AMCD surveillance records found evidence of a shift in the population abundance and distribution of Ae. taeniorhynchus mosquitoes in permanent trapping locations following LULC changes (Qualls et al. 2021). Intraspecies competition is another gray area likely to confound results, since its impacts at the county level are unclear, while significant relationships have been distinguished at a trap-site level (Qualls et al., unpublished data). The significant negative trend of Ae. albopictus and Cx. quinquefasciatus versus year found during correlation analyses is an interesting observation and should be examined further in context of LULC change and climate change. In the context of implementing these exact models for Anopheles, Culex, and Aedes species, it is important to consider that validation with AMCD data showed a reduced predictive capability as mosquito populations increased. Florida and the AMCD already have a legislative mosquito action threshold in place. The models here tended to overestimate and lead to false positives during validation; thus if AMCD realistically moved by these models it could lead to wasted time, money, and labor resources. A follow-up analysis is planned both to investigate the quantitative impact of AMCD ULV adulticide missions on mosquito population numbers and to understand the overall pattern of theoretically triggered and/or completed missions. Another future direction incorporating control data is to investigate the effectiveness of control activities on directly targeted species (>25 trap count) and indirectly targeted species (<25 trap count) to understand if a lowered threshold number effectively reduces populations sooner and for longer periods. Otherwise, we could rely on model equations to back calculate using 25 and 8.33 as the response variable to understand the threshold precipitation and temperature values of concern.

The final takeaway is that these models are useful despite the limitations and necessary assumptions. It is meaningful that statistically significant models were successfully developed for ecologically distinct habitat genera groups, albeit with subjective predictive capacity and operational capabilities. Importantly, no actual thresholds were optimized due to the analytical limitations described and the logistical hurdles to investigate more robust models given the dataset constraints. In the context of implementation of these models or the development of models using the same methodology, it is crucial for any program to account for available resources (personnel, time, equipment) and field knowledge (chemicals, mosquito biology/ecology) in the context of regional priorities (political, economic, disease transmission). Questions for any novel control initiative to address include (1) presence of which nuisance and/or vector species, 2) potential of disease transmission in the area, 3) physical and financial ability to conduct surveillance activities, 4) physical and financial ability to conduct surveillance activities, and 5) typical climate patterns. Programs lacking historical data or rigorous entomological training can further consider utilization of online databank platforms that enable access to mosquito and mosquito-borne disease data (Brown et al. 2021). Leverage of existing datasets can 1) provide knowledge of local species and disease incidence, and 2) lend local or geographically similar data for threshold model development. Leverage of climate forecast models may lend sufficient lead time to provide climate predictor estimates (Lowe et al. 2018). These models for SJC do not provide a silver bullet answer to each specific context worldwide but do lend thought to how to optimize models for distinct regions and resource-dependent contexts.

Precipitation was the most important predictor for SJC, potentially because of its subtropical climate, whereas temperature and especially relative humidity might be of larger significance in other eco-regions due to larger seasonal variations. The variables of rainfall amount versus number of rainy days have subtle differences in statistical power and biological significance, and attention should be paid to the occurrence and lagged effect of outlier rainfall events. Lag time and temperature indices were important when constructing models and seasonality is an important modifier. Sensitivity and specificity should be balanced during validation of any threshold before implementation. Any control activities against a target species are going to impact nontarget species sharing the same ecological niche. Optimizing models for priority species rather than all species can help to fulfill reduction of disease transmission in conjunction with nuisance nonvector species. Operational thresholds need to be proactive and geared toward specific control activities, since timing and strategy strongly depend on target species, life stage, and habitat. Leverage of existing mosquitoes, mosquito-borne disease, and climate data during model development, threshold optimization, or ongoing surveillance can increase operational efficiency.

This research project was funded by the Department of the Army, U.S. Army Contracting Command, Aberdeen Providing Ground, Edgewood Contracting Division, Ft. Detrick, MD, under Deployed War Fighter Protection (DWFP) program grant W911QY20110004 awarded to the AMCD. The funders had no role in data study design, data collection, and analysis.