Study area
The study was conducted in the WVMVCD of southwestern San Bernardino County, encompassing six cities—Chino, Chino Hills, Ontario, Upland, Montclair and Rancho Cucamonga—with a total estimated population of 600,000 (Fig. 1). The area is characterized by a Mediterranean climate, experiencing hot, dry summers and mild, wet winters [26]. Seasonal precipitation occurs predominantly between November and April, with occasional summer thunderstorms. Summer temperatures in the WVMVCD frequently surpass 38 °C, with occasional spikes exceeding 42 °C (Supplementary Figures S1 and S2). The native vegetation includes chaparral, oak woodlands, and riparian habitats, although urbanization has led to significant alterations in the landscape [27]. Furthermore, the region has experienced rapid population growth in recent decades, accompanied by extensive urban development, including residential, commercial, and industrial structures.
Adult mosquito trapping
Mosquito surveillance was carried out by the WVMVCD, a special district responsible for vector surveillance and control in the region. Mosquito collections took place from January 2017 to December 2023 using BG-Sentinel traps (Biogents AG, Regensburg, Germany). BG-Sentinel traps are recognized for their effectiveness in monitoring Ae. aegypti across spatial and temporal trends [28]. Trapping sites were selected based on a combination of criteria: locations from which mosquito-related service requests had been recently received by the vector control district, and sites with a history of high Aedes abundance based on prior surveillance records. In addition, exploratory trapping was conducted in areas with no historical data to capture current mosquito dynamics. This targeted approach ensured optimal allocation of surveillance resources to areas of greatest public health concern and potential mosquito activity. Traps were deployed in shaded, wind-sheltered outdoor locations in close proximity to human habitations or areas of known mosquito activity, while avoiding direct sunlight and potential competing attractants. Each trap was equipped with a BG-Lure, a proprietary attractant that mimics human skin odors, placed in the designated lure compartment to enhance mosquito attraction. To further increase trap efficacy, an additional attractant, delivered via dry ice placed in insulated containers positioned adjacent to the trap, was included. Traps were powered by 12 V rechargeable batteries and operated continuously for a 24 h period. The following day, collection nets containing the captured mosquitoes were carefully removed from the traps, labeled, and transported in cool, insulated containers to preserve specimen integrity. In the laboratory, adult mosquitoes were knocked down by dry ice exposure and subsequently identified morphologically to species level under a microscope. Relevant data including Global Positioning System (GPS) coordinates were recorded for each sampling event to support interpretation of trap yields and spatio-temporal trends.
Trapping was systematically conducted throughout the district and set in residential areas. With the homeowner’s consent, traps were left overnight. Mosquitoes were collected the next morning and then counted and sorted by species. Each trap included a geographical location (latitude and longitude) and a timestamp indicating when it was set. Trap sites were strategically chosen based on historical mosquito hotspots, service requests from residents (e.g., complaints about mosquito biting), and an exploratory strategy to monitor mosquito activity across different areas of the district. All traps were brought to the WVMVCD Lab for mosquito counting and species identification by trained technicians. The district maintains records of adult mosquito counts, categorized by location and date, from 2017 to 2023.
The number of trapping sites varied from week to week depending on the number of Aedes mosquitoes collected at each location. Only traps that captured more than 20 Ae. aegypti individuals during the previous collection were eligible for redeployment, unless they were designated as routine sites. Additionally, mosquito-related service requests from residents informed trap placement, contributing to spatial variability across the study period. Notwithstanding these changes, eight routine sites were monitored weekly and consistently from 2017 to 2023 (see Supplementary Figure S10).
Environmental and meteorological variables
Environmental and meteorological data were obtained through the Google Earth Engine (Table 1) from a variety of Earth observation data sources [29]. We extracted data on the built environment, elevation, precipitation, Normalized Difference Water Index (NDWI, a measure of surface water or moisture), and average ambient temperature—all chosen for their hypothesized or known impact on mosquito abundance (listed in detail in Table 1).
These data were downloaded and merged into the trapping dataset based on location and trap date. Buffers 150 m in size were drawn around each of the trap locations, and environmental and meteorological variables were attributed to the traps based on these buffers. This buffer size was chosen based on the assumption that Ae. aegypti dispersal is generally limited to within 100 m of their breeding sites [30]. In our exploratory analyses, we also tested the use of 250-m buffers, based on a study from central California that suggested this species may have broader dispersal in nearby areas [31].
We extracted precipitation, surface water, and average temperature with 14- and 28-day lag periods from the mosquito collection date to account for the time needed for these environmental factors to influence mosquito development and population dynamics [32].
Statistical analysis
Temporal patterns in trapped Aedes
We used descriptive statistics to summarize Ae. aegypti trap data from 2017 to 2023, including annual totals, means, standard deviations, medians, interquartile ranges, and the proportion of traps with zero mosquitoes. We investigated counts per trap by year to assess long-term trends in abundance, and by month across all years to detect seasonal patterns (results in Table 2 and Supplemental Table S1). These descriptive analyses helped contextualize environmental and spatial trends explored in later modeling.
Spatial patterns in trapped Aedes
We mapped annual Ae. aegypti counts per trap from 2017 to 2023 to examine spatial patterns and changes in mosquito abundance over time. We then employed two main approaches to examine the potential spatial patterns of Ae. aegypti abundance from 2017 to 2023. First, we created spatial correlograms for each year to investigate how spatial autocorrelation in mosquito abundance changed with increasing distance between trap locations. These correlograms were constructed by calculating Moran’s I at increasing distance intervals up to 3 km, using 500-m bins. To account for sampling effort variations, we used the mean weekly number of mosquitoes per trap location within each year for these calculations. This approach provided insights into the scale of spatial dependence and its annual variations.
Following this global assessment, we employed the Getis-Ord Gi* statistic to identify significant local clusters of high (hot spots) and low (cold spots) Ae. aegypti abundance. We applied a fixed-distance spatial weights matrix with a 500-m threshold, selected based on results from the correlogram analysis. Locations with statistically significant high positive Z-scores (> 1.96, P < 0.05) were identified as hot spots, while those with statistically significant low negative Z-scores (< −1.96, P < 0.05) were classified as cold spots.
To ensure robustness, we conducted a sensitivity analysis by repeating the correlogram analysis and Getis-Ord Gi* using 250-m bins/threshold, comparing results with our primary analysis (Supplementary Figures S4 and S5). We generated annual hotspot maps and created a combined visualization of annual spatial correlograms.
This multi-approach method allowed us to track changes in spatial clustering patterns over time. Our exploratory spatial analyses provided a foundation for investigations into environmental correlates of Ae. aegypti abundance.
Bivariate analysis between environmental and meteorological covariates and Ae. aegypti abundance
We hypothesized that the spatial and temporal distribution of Ae. aegypti in the WVMVCD would be positively associated with precipitation, surface water, and the built environment. We also hypothesized that higher ambient temperatures would be associated with higher mosquito abundance. To visualize potential associations between Ae. aegypti abundance and environmental factors, mosquito abundance data from traps were plotted against continuous environmental variables. This approach allowed for a preliminary examination of potential associations between each factor and Ae. aegypti abundance in our study area.
Generalized additive model for associations between environmental and meteorological covariates and Ae. aegypti abundance
We then used a negative binomial generalized additive model (GAM) to model counts of trapped Ae. aegypti mosquitoes, using the environmental variables as covariates in the model and accounting for trap location and time, for the duration of the 2017–2023 study period. The GAM allowed us to account for the complex and potentially non-linear associations between environmental variables and Ae. aegypti abundance through the use of smoothing spline functions [33].
To account for potential lagged effects (i.e., precipitation might have a delayed effect on mosquito abundance), we explored different lags (14 and 28 days) for the precipitation, surface water, and average temperature. We began by developing a model for all data from 2017 through 2023, with all covariates specified as thin plate regression spline functions. All covariates were centered on their means and standardized by their standard deviation to facilitate the interpretation and comparison of effect sizes.
Final variable selection and model specification was based on a combination of scientific hypotheses, model goodness-of-fit tests (deviance explained and the Akaike information criterion (AIC)] (Supplementary Tables S2 and S3), and interpretability—referring to the inclusion of variables with meaningful, biologically relevant associations that could inform vector control strategies. Through comparison of model fit using deviance explained and the AIC, we determined that a 14-day lag for these variables provided the best model fit. The final negative binomial GAM included the following variables: NDWI, precipitation, built environment, elevation, average temperature, geographical location, season, and year. It also incorporated interaction terms between season and average temperature, as the effect of temperature shifted over time.
The results of our regression analysis are presented as graphs showing the estimated smoothed spline functions for each predictor variable (Fig. 6), along with a summary table indicating statistical significance (Supplementary Table S2).
Software
All statistical analyses were performed using the R statistical software (version 4.3.3). The spatial autocorrelation analyses were run using the spdep and tmap packages. The GAM was created using the mgcv package. All maps were created using QGIS (version 3.34.9).