Modeling the association between Aedes aegypti ovitrap egg counts, multi-scale remotely sensed environmental data and arboviral cases at Puntarenas, Costa Rica (2017–2018)

Problems with vector surveillance are a major barrier for the effective control of vector-borne disease transmission through Latin America. Here, we present results from a 80-week longitudinal study where Aedes aegypti (L.) (Diptera: Culicidae) ovitraps were monitored weekly at 92 locations in Puntarenas, a coastal city in Costa Rica with syndemic Zika, chikungunya and dengue transmission. We used separate models to investigate the association of either Ae. aegypti-borne arboviral cases or Ae. aegypti egg counts with remotely sensed environmental variables. We also evaluated whether Ae. aegypti-borne arboviral cases were associated with Ae. aegypti egg counts. Using cross-correlation and time series modeling, we found that arboviral cases were not significantly associated with Ae. aegypti egg counts. Through model selection we found that cases had a non-linear response to multi-scale (1-km and 30-m resolution) measurements of temperature standard deviation (SD) with a lag of up to 4 weeks, while simultaneously increasing with finely-grained NDVI (30-m resolution). Meanwhile, median ovitrap Ae. aegypti egg counts increased, and respectively decreased, with temperature SD (1-km resolution) and EVI (30-m resolution) with a lag of 6 weeks. A synchrony analysis showed that egg counts had a travelling wave pattern, with synchrony showing cyclic changes with distance, a pattern not observed in remotely sensed data with 30-m and 10-m resolution. Spatially, using generalized additive models, we found that eggs were more abundant at locations with higher temperatures and where EVI was leptokurtic during the study period. Our results suggest that, in Puntarenas, remotely sensed environmental variables are associated with both Ae. aegypti-borne arbovirus transmission and Ae. aegypti egg counts from ovitraps.

In Costa Rica, dengue, chikungunya and Zika are the most common arboviruses affecting humans, and account for the highest vector-borne disease burden in the country, given their high case number, which is in the magnitude of thousands of annual cases (Soto-Garita et al., 2016;Vigilancia de la Salud, 2019). This case load is well above what is observed for other vector-borne diseases, such as cutaneous leishmaniasis, where annual cases amount to a few hundred (Chaves et al., 2008), or malaria, which is on the brink of elimination (Chaves et al., 2020a). In Costa Rica, dengue transmission has been traditionally located at low altitudes, on coastal regions from the Pacific and Caribbean basins (Mena et al., 2011). On the Pacific basin, the county of Puntarenas has experienced high transmission of Ae. aegypti-borne pathogens over recent decades. For example, a serological survey found over 90% of the people exposed to dengue virus (Lee-Lui et al., 2008). Moreover, Puntarenas county is regularly among the top 10% of counties reporting dengue cases in Costa Rica (Vigilancia de la Salud, 2019).
Aedes aegypti has a long history in Puntarenas and was documented to be present by 1920 (Alfaro, 1921;Serre, 1921) and the late part of the 1930s (Kumm et al., 1940). In 1961, the country was certified as free from Ae. aegypti as part of yellow fever eradication efforts (Soper, 1967), a status held until at least 1988 (Gubler, 1989). However, in 1993 local dengue transmission started, implying that Ae. aegypti was again present in Puntarenas (S aenz et al., 1999). Recent studies have identified that key larval habitats of Ae. aegypti in Puntarenas include artificial discarded containers, but also water holding tanks, and other containers with active domestic use Marín Rodríguez & Díaz Ríos, 2012). A couple of studies have described insecticide resistance patterns in Ae. aegypti from Puntarenas. Bisset et al. (2013) reported resistance to the pyrethroid deltamethrin and to the organophosphate temephos used for larval control. Zardkoohi et al. (2020) confirmed the resistance to deltamethrin, and the emergence of resistance to the pyrethroid cypermethrin, linking pyrethroid resistance to the co-occurrence of the V1016I and F1534C kdr mutations in the voltage-gated sodium channel gene.
Several studies have linked weather changes with both Ae. aegypti abundance and arboviral pathogen transmission (Scott et al., 2000;Johansson et al., 2009;Chaves, 2017a;Reinhold et al., 2018). However, no attempt has been made at linking temporal changes in Ae. aegypti-borne arboviral transmission, and Ae. aegypti abundance, with weather variables at Puntarenas, which is, and has historically been, a major focus for dengue transmission in Costa Rica. This potential association is particularly important because a global meta-analysis found limited evidence for dengue transmission prediction based on data from entomological surveys of Ae. aegypti aquatic stages (Bowman et al., 2014). One major problem to test if dengue transmission was associated with entomological indicators was the lack of a spatially standardized sampling (Bowman et al., 2014), which is more feasible with ovitraps (CEN-APRECE, 2015). In Costa Rica, ovitraps have been successfully used to survey mosquito biodiversity (Chaverri et al., 2018) and to study mosquito population dynamics (Romero et al., 2019). Moreover, the Costa Rican National programme for integrated vector management started to use ovitraps to guide vector control in 2017, as done elsewhere in Latin America (Regis et al., 2008;CENAPRECE, 2015;PAHO, 2017). All these conditions make Puntarenas an ideal site to test if arboviral cases can be associated with Ae. aegypti egg counts from ovitraps.
Here, we present the results of an 80-week longitudinal study where weekly observations on the number of eggs at 92 ovitrap locations were recorded between 2017 and 2018. We tested whether remotely sensed environmental data and Ae. aegypti egg counts from ovitraps can be associated with arboviral transmission in Puntarenas, a city on the Pacific Coast of Costa Rica. Additionally, we also tested whether egg counts were associated with remotely sensed environmental variables.
For our analysis, we employed satellite images with different resolution. Specifically, we included Sentinel 2 (10 m resolution), Landsat 8 (30 m) and MODIS (250 m and 1 km). We used several temporal, spatiotemporal and spatial modelling techniques. We started by (i) performing an ancillary analysis comparing land surface temperature estimates obtained with MODIS images for the Puntarenas city area, where cases were recorded, with estimates based on Landsat 8 images for the Puntarenas peninsula area, the subsection of Puntarenas city where most of the transmission is assumed to occur (Marín Rodríguez & Díaz Ríos, 2012).
We then developed (ii) time series models to assess the association between the remotely sensed environmental variables with Ae. aegypti egg counts and Ae. aegypti-borne arboviral cases. We specifically employed seasonal autoregressive (SAR) models given their ease for fitting and interpretation (Chaves & Pascual, 2007). Among other things, these models allow testing for the significance of environmental variables at time lags whose association with the studied time series is not an artifact of a similar, but unrelated, seasonality (Priestley, 1988). These models allow to evaluate if the arboviral cases were associated with egg counts and remotely sensed environmental variables at different time lags. This analytic framework was also used to assess the association between egg counts and remotely sensed environmental variables at different lags.
We also performed (iii) a synchrony analysis of Ae. aegypti egg counts and the remotely sensed environmental variables based on Landsat 8 and Sentinel 2 images. This is an innovative analysis to study spatio-temporal patterns in mosquito oviposition, which we have previously used to study mosquito larval (Chaves et al., 2020b) and adult (Chaves, 2017b) abundance patterns. Briefly, a synchrony analysis shows the degree of concerted fluctuation between values of a studied variable, egg counts in this study, as a function of the distance separating the sampling locations (Chaves et al., 2013). This analysis allows to test if the correlations between egg counts are synchronous, i.e. not changing with distance. Alternatively, the correlations can decay with distance as expected under dispersal limitation, or have a fluctuating sign, changing from positive to negative, with distance, as expected when dispersal occur in traveling waves (Ranta et al., 2006). Thus, the synchrony analysis provides a description of the type of spatial heterogeneity in oviposition, measured through time. The obtained oviposition synchrony can also be compared with the synchrony in environmental data to establish a dependence of the former on the latter (Chaves et al., 2013).
We finished by performing a (iv) spatial analysis of the mean number of Ae. aegypti eggs at each sampled location with average values for the remotely sensed environmental variables, based on Landsat 8 and Sentinel 2 images, at the pixels where traps were located. For this analysis, we employed generalized additive models. This modelling framework enables to easily incorporate nonlinearities in the association between variables (Faraway, 2006). This analysis allowed to test if average egg counts were associated with average environmental conditions estimated from finely-grained (30 m and 10 m) remotely sensed data.

Study site
Our study was carried out in Puntarenas city (9 58 0 34.50 00 N, 84 50 0 18.10 00 W), the largest city on the Pacific coast of Costa Rica (Fig. 1), capital of Puntarenas county and Puntarenas province. The urban area of Puntarenas city is formed by Puntarenas, Chacarita, Barrancas, El Roble and Pitahaya districts (Fig. 1). The city has a population of 81,187 inhabitants, according to the last census (2011), and 193.98 km 2 of area, for a population density of 418 inhabitants/km 2 (INEC, 2012). Puntarenas peninsula, the peninsular part of Puntarenas district (Fig. 1), concentrates government, education and commercial activities of Puntarenas city (Villanueva, 2009). Bus and ferry terminals that serve as transportation hubs for Puntarenas city, Puntarenas county and the Pacific basin of Costa Rica are also located in the Puntarenas peninsula .
Puntarenas city has a marked seasonal climatic pattern with a dry season from December to April, and a rainy season over the rest of the year (World Meteorological Organization, 2020). The average annual total cumulative rainfall is 1,600 mm, while temperature has an annual average of 32.8 C, which is about 2.1 C hotter during the dry season than the rainy season (World Meteorological Organization, 2020).

Mosquito sampling
Mosquitoes were sampled using ovitraps made with PVC plastic. The ovitraps have a cylindrical shape (diameter of 11 cm and height of 14 cm), with a closed and an open end ( Fig. 2A) and were donated by the Mexican Government through the Center for Infectious Disease Control and Prevention (CENAPRECE), a branch of Mexico's Secretary of Health. The ovitraps have two 5 mm in diameter openings 8 cm above the closed end, thus having a water holding volume of 750 ml. Ovitraps were set using Scott® paper towels (Kimberly Clark Co., Neenah, WI, USA), which were placed around the internal side of the ovitrap. Following the recommendations from CENAPRECE (2015) we filled the ovitraps with 750 ml of tap water. This water came from the city water distribution network, where chlorine concentrations are strictly monitored and fluctuate between 0.3 and 0.6 mg/l (Marín Mena, 2007). The paper towels  A Ovitrap label indicating its use for dengue vector surveillance and providing contact information for the Programa de Manejo Integrado de Vectores of Costa Ricaʼs Ministry of Health. B Field deployment of an ovitrap, next to a plant pot, a common larval habitat for Aedes aegypti in Puntarenas, Costa Rica were added to sample Ae. aegypti eggs, given the oviposition behavior of this species which places eggs on surfaces above the water line (Chadee & Corbet, 1987). Traps were labeled ( Fig. 2A) and deployed (Fig. 2B) following a protocol developed by the CENAPRECE (2015).
A total of 92 traps were deployed on April 17th, 2017 (epidemiological week 16) through the Puntarenas peninsula (Fig. 3). Both the water and the paper towels were replaced every Monday from epidemiological week 17 of 2017 (April 24th 2017) to epidemiological week 44 of 2018 (October 29th 2018). This was done to avoid biases in mosquito sampling related to ovitrap water age and/or conspecific presence (Chadee, 2009). Ovitraps were placed inside the houses (n ¼ 90) or outside the houses but within the household space, e.g. the backyard, specifically an area outside the houses where a sink for washing clothes ("pila" in the Spanish of Costa Rica) is traditionally located (n ¼ 2). Ovitraps were placed next to the walls defining the residential premises, i.e. within approximately within one meter from the wall (within 1 m), in places selected by the homeowners, but restricted to the first floor of multi-story buildings. Traps were deployed following CENAPRECE (2015) recommendations, i.e. a maximum of 4 ovitraps were deployed by residential block, "manzana" in Spanish, trying to ensure a minimal distance between any pair of traps was 20 m, something done because of the finely-grained nature of Ae. aegypti oviposition . Homeowners were asked for consent to place the ovitraps in their households as part of efforts for arbovirus transmission surveillance.
In light of insights about dengue transmission in public spaces with a high turnover of visitants and human movement (Stoddard et al., 2009;Arunachalam et al., 2010;Smith et al., 2014), ovitraps were located within the Puntarenas peninsula. This decision was made given that Puntarenas peninsula serves as a hub for human movement within Puntarenas city  and because previous studies have shown persistent mosquito infestations in the area (Marín Rodríguez & Díaz Ríos, 2012) suggesting that the peninsula is a place where most Ae. aegypti-borne arboviruses are transmitted among inhabitants of the larger Puntarenas city (Marín Rodríguez & Díaz Ríos, 2012).
Ovitraps were georeferenced using the recorded address for each ovitrap with Google maps. Paper towels were individually transported inside identified medium sized (17.7 Â 18.8 cm) freezer Ziploc® bags (SC Johnson Co., San Mateo Otzacatipan, M exico) to the Coordinaci on Regional of Programa Nacional de Manejo Integrado de Vectores in Puntarenas, Costa Rica, where eggs were counted under a 10Â dissecting microscope. The counting was carried out under a dissecting microscope to ensure that eggs of Limatus spp., if present, were separated from Aedes spp. eggs, given the striking morphological differences between these two mosquito genera and the potential colonization of ovitraps by Limatus durhami Theobald, a species present in the study area (Calder on-Arguedas et al., 2009), and considering that Limatus spp. also oviposit on surfaces above the water line (Santos Neto & Azevedo Marques, 1996;Ceretti-Júnior et al., 2014;Carvalho et al., 2017). Briefly, Limatus spp. eggs have a shape similar to a boomerang (or elongated rhombus), as opposed to the elongated elliptical shape of Aedes spp. eggs (Santos Neto & Azevedo Marques, 1996). Other species in the studied area (Vargas, 1961;Calder on-Arguedas et al., 2009) that could potentially colonize ovitraps included Culex spp. and Uranotaenia spp., both of which oviposit in rafts and whose eggs are projected into the water surface (Harbach & Knight, 1978;Day, 2016), and Toxorhynchites spp., which oviposit single eggs, that are significantly larger than eggs from other mosquito genera and that are also projected into the water surface . For further mosquito identification quality control, each week up to 100 eggs from 3 randomly chosen traps were allowed to hatch to verify the identity of the species collected. Throughout the study, eggs from each ovitrap location were checked at least once; all samples collected belonging to Ae. aegypti. In addition, Puntarenas peninsula is an area without current or historical records for the presence of Aedes albopictus (Skuse) (Futami et al., 2015). During this study, we did not find Limatus spp. eggs while examining and counting sampled eggs under the dissecting microscope.

Arboviral disease data
Weekly counts for Zika, dengue and chikungunya cases from Puntarenas county, for all the weeks of 2017 and 2018, were obtained from the weekly epidemiological bulletins published by the epidemic surveillance division (Vigilancia de la Salud) from Costa Rica's Ministry of Health which are publicly available at: https://www.ministeriodesalud. go.cr/index.php/vigilancia-de-la-salud/boletines/enfermedades-de-tra nsmision-vectorial-2017. These case counts include all clinically diagnosed dengue cases, of which a variable proportion, ensuring all presumed transmission clusters are sampled, are confirmed by PCR and/or serological methods by the Costa Rican National Reference Laboratory for viral infections (Gonz alez Elizondo, 2018). All Zika and chikungunya cases are confirmed by laboratory methods (PCR and/or serology), before being consolidated and reported by the epidemic surveillance division at Costa Rica's Ministry of Health. Independently of the moment when the final diagnostic is performed, cases are recorded for the week individuals attend health facilities with symptoms of an arboviral infection. We downloaded daily images for land surface temperature and emissivity (LST&E) with a spatial resolution of 1 km (MOD1101, Version 6) based on Moderate Resolution Imaging Spectroradiometer (MODIS) images (Wan et al., 2015). The images were re-scaled by multiplying each pixel by a factor of 0.02, and the resulting temperatures were transformed from K to C by subtracting 273 (Wan et al., 2015). Downloaded MODIS images can be seen in the supplementary online Video S1. We generated Landsat 8 land surface temperature images with a spatial resolution of 30 m using a script for online use with Google Earth Engine (Ermida et al., 2020). Downloaded Landsat 8 images can be seen in the supplementary online Video S2 and were available every 8 days.
Supplementary videos related to this article can be found at https ://doi.org/10.1016/j.crpvbd.2021. 100014 We also downloaded images for the Normalized Difference Vegetation Index (NDVI) and Enhanced Vegetation Index (EVI), both of which are considered proxies for vegetation growth (Pettorelli et al., 2005) and have been extensively used to study the ecology of mosquitoes (Hurtado et al., 2018b;Chaves et al., 2019;Poh et al., 2019;Rigg et al., 2019;Nguyen et al., 2020;Chaves & Friberg, 2021) and other insect vectors (Kitron et al., 1996). The images for vegetation indices were downloaded for the multiple spatial scales encompassed by MODIS, Landsat 8 and Sentinel 2. The MODIS NDVI and EVI images were processed every 16 days, have a spatial resolution of 250 m (MOD13Q1, version 6 and MYD13Q1, version 6) and are collected by Terra and Aqua satellites (Didan, 2015a, b). These images were re-scaled by dividing the values of each pixel by 10,000, thus obtaining normalized values that range between À1 and 1 (Pettorelli et al., 2005). We used data from Terra (MOD13Q1) and Aqua (MYD13Q1) satellites, to increase the number of observations used to estimate environmental time series. This also allowed us to increase the frequency of observations to every 8 days provided the delay between the images from each satellite for the study site.
We also used information from the red (RED, band 4) and near infrared (NIR, band 5) Landsat 8 bands, with a resolution of 30 m, downloaded with the Landsat 8 land surface temperature to estimate NDVI using the following equation (Vermote et al., 2016): Adding information from the Landsat 8 blue band (BLUE, band 2), we estimated EVI using the following equation (Vermote et al., 2016): The resulting images were used to estimate Landsat 8-based NDVI and EVI. Equations (1) and (2) were also used to estimate NDVI and EVI using top of the atmosphere Sentinel 2 images, which have a resolution of 10 m and where BLUE is band 2, RED is band 4 and NIR is band 8 (Drusch et al., 2012). Sentinel 2 images for the study period were processed using the Google Earth Engine to remove pixels covered by clouds or whose values exceeded the theoretical bounding values of À1 to 1 (ESA, 2015). Downloaded Sentinel 2 images can be seen in the supplementary online Video S3, and were available every 3.4 AE 2.0 days.
Supplementary videos related to this article can be found at https ://doi.org/10.1016/j.crpvbd.2021.100014 To generate MODIS-based temperature, NDVI and EVI time series, each image was cropped using a polygon from the joint area of the 5 districts conforming the urban area of Puntarenas ( Fig. 1) using the package "raster" in R (Brunsdon & Comber, 2015) and mean values were estimated for the selected pixels containing information. We repeated this process for the Landsat 8-based temperature images, and Landsat 8and Sentinel 2-based NDVI and EVI images, but cropping only the area corresponding to the Puntarenas peninsula (Fig. 1). The difference in the cropped areas was due to three main reasons: (i) over 90% of the arboviral cases occur in the Puntarenas city area of Puntarenas county, with over 80% of the cases linked with activities in the Puntarenas peninsula area (Vigilancia de la Salud, 2019) and such area is more easily covered by the lower resolution MODIS images; (ii) traps were deployed in the Puntarenas peninsula area, at distances where the spatial resolution of Landsat 8 and Sentinel 2 images is more likely to detect differences in the local environments were traps were located (which had a mean AE standard deviation (SD) of 495 AE 311 m, Supplementary Fig. S1); and (iii) having multi-scale spatial images allowed us to test the hypothesis that both local and regional environmental factors might be involved in arbovirus transmission and Ae. aegypti oviposition. In addition, for MODIS and Landsat 8 temperature images we estimated time series for the SD and kurtosis, respectively, in the Puntarenas city and peninsula area, given the possibility that arboviral cases and egg counts are sensitive to patterns of variability in the environment, as predicted by Schmalhausen's law, the biological principle stating that organisms follow both the mean and higher order moments, i.e. SD and kurtosis, of environmental variation (Chaves, 2017a).
The resulting time series were then smoothed using the loess algorithm for local polynomial regression fitting (Cleveland & Devlin, 1988) by performing a regression of the estimated temperature, NDVI or EVI as a function of the day when the images were acquired by the satellites. For the loess, we employed a neighborhood size of 10% of the data and second-degree polynomials.
The resulting smoothed time series were then used to obtain weekly estimates of NDVI and EVI, extrapolating values to the Monday of each week in the case of NDVI and EVI from MODIS, Landsat 8 and Sentinel 2 images, and also for land surface temperatures, including mean, SD and kurtosis from Landsat 8 images. However, for MODIS-based temperatures, we estimated weekly average, SD and kurtosis from Tuesday to Monday, using the daily temperature time series.
All MODIS and Landsat 8 images were courtesy of the NASA Land Processes Distributed Active Archive Center (LP DAAC) and the United States Geological Survey (USGS)/Earth Resources Observation and Science (EROS) Center (Sioux Falls, South Dakota). Copernicus Sentinel 2 data were processed, and provided free of charge, by the European Space Agency (ESA). All MODIS images were downloaded from the LP DAAC server (NASA, 2018) using the package MODIStsp in R (Busetto & Ranghetti, 2016). Landsat 8 and Sentinel 2 images were retrieved from the Google Earth Engine Dataset Catalog (Google, 2020).

Statistical analysis
Our statistical analysis had four main components.

Comparison of land surface temperature (LST) time series estimates with MODIS and Landsat 8 images
We examined LST time series estimated with MODIS ( Supplementary  Fig. S2A) and Landsat 8 ( Supplementary Fig. S2B) using the Pearson's correlation coefficient (Zar, 1998). We also compared MODIS and Landsat 8 data examining their coefficient of variation (CV). The CV, formally defined as the ratio between the SD and the mean of a variable (Stearns, 1981), was estimated for each observation in the MODIS ( Supplementary Fig. S2C) and Landsat 8 ( Supplementary Fig. S2D) time series. Then we assessed any potential bias in the estimation of time series observations related to the number of available pixels by performing a linear regression of the LST CV as a function of the number of pixels used to estimate LST. We used MODIS images for the Puntarenas city area ( Supplementary Fig. S2E) and Landsat 8 images ( Supplementary  Fig. S2F) for the Puntarenas peninsula area when estimating the regressions.

Time series modeling
We studied the association between Ae. aegypti-borne arboviruses (Zika, chikungunya and dengue) and Ae. aegypti egg counts and remotely sensed environmental variables using the general framework of seasonal autoregressive (SAR) models (Shumway & Stoffer, 2011). SAR models have the general form: where X t is the focal variable, i.e., arboviral cases or mean number of eggs in this study, observed at time t, X t is a function of itself at time t-j, and at time t-s, where j and s are the time lags for the "autoregressive" and "seasonal" parameters (φ j and φ s ; respectively) that account for autoregressive and cyclical patterns in the time series of the focal variable studied, and Cov(i) indicates covariates, which are related to the focal variable under study according to parameters β i . The SAR models can have any k number of covariates, as long as this number is shorter than the time series length, minus the number of seasonal and autoregressive parameters plus one, μ is the mean of the time series, and ε t is the error, which is assumed to be normal, independent and identically distributed (Shumway & Stoffer, 2011). To develop SAR models, we started by the identification of the maximum time lags for the autoregressive and seasonal components of the model. For lag identification we employed autocorrelation (ACF) and partial autocorrelation (PACF) functions which show the correlation of time series observations at different time lags considering, respectively, the full time series or only consecutive lags (Hoshi et al., 2014). Based on information from these functions, a null SAR model was fitted and used to pre-whiten time series of covariates whose association is then evaluated by examining cross-correlation functions (CCFs) between the covariate time series and the focal time series. This analysis allows the identification of lags at which covariates are correlated with the studied variables. Pre-whitening was used to avoid the spurious identification of significant lags of association that emerge from time series having similar autocorrelation structures (Hoshi et al., 2017). Models for the weekly combined case count of dengue, chikungunya and Zika considered as potential covariates the median, the SD and kurtosis of weekly collected eggs across all ovitraps, the mean of remotely sensed temperature, NDVI and EVI at the different spatial resolutions considered in this study. We also studied the median number of weekly collected eggs as a function of the mean of remotely sensed temperature, NDVI and EVI at the different spatial resolutions considered in this study. In both cases we considered the MODIS-based time series estimated for the whole Puntarenas city and the Landsat 8-and Sentinel 2-based time series estimated for the Puntarenas peninsula. This was done to test if environmental conditions specific to the area where most transmission is assumed to occur, i.e., the Puntarenas peninsula, were associated with transmission recorded for the larger Puntarenas city area in models for arboviral cases. Similarly, in models for Ae. aegypti egg counts we considered the MODIS-based environmental time series for the Puntarenas city area to test if environmental phenomena occurring at larger spatial scales could influence Ae. aegypti oviposition dynamics recorded in the Puntarenas peninsula.
When fitting the models, all selected covariates were demeaned to ease their interpretation in terms of changes above or below their means . We selected the best models through the minimization of the Akaike Information Criterion (AIC), a metric that trades-off model goodness-of-fit and parameter number, employing a strategy of backward elimination, where models are simplified by comparing "nested" models with the same number of parameters that are not significantly different from the "parent" model that considered all variables that are left out in each of the simplified "nested" models (Chaves, 2016;, or until there were no significant differences through a Chi-square likelihood ratio test between the model minimizing AIC and simplified versions of such model .

Synchrony analysis
We performed a synchrony analysis for Ae. aegypti oviposition and environmental variables estimated with Landsat 8 and Sentinel 2 images. This analysis was carried out to understand if concerted fluctuations in eggs/ovitrap during the study period followed similar changes in temperature and the vegetation indices that we studied. For the analysis, we employed a Mantel correlogram with inference based on a Monte Carlo test, whose null hypothesis is that synchrony across the range of distances separating sampling locations is equal to a global average, unless more extreme (Chaves, 2017b;Chaves et al., 2020b). For the analysis, we only considered up to 50% of the spatial extent of the data to avoid spurious results at long distances due to small sample size (Gouhier & Guichard, 2014), but estimated regional synchrony considering samples from all the georeferenced ovitraps that had, at most, one missing observation. To estimate the synchrony, we used the time series from the 44 ovitrap locations that did not have (n ¼ 27) or had at most one missing observation (n ¼ 17), which were inputted by taking the average of the two immediate observations, i.e. one time step before and after the missing observation, in order to increase the power for statistical inference (Gouhier & Guichard, 2014). We also estimated synchrony for environmental variables at the pixels that contained the ovitraps considered in the oviposition synchrony analysis. For this we extracted Land Surface Temperature (LST), NDVI and EVI data from Landsat 8 images, and NDVI and EVI data from Sentinel 2 images, at the pixels containing the 44 ovitraps used for the Ae. aegypti oviposition synchrony analysis. We inputted missing values for the resulting environmental time series following the procedure described in the time series modeling section. The 44 ovitraps used for this analysis did not include any of the two ovitraps placed outside the houses.

Spatial analysis
We tested whether mean number of Ae. aegypti egg counts at georeferenced ovitraps (n ¼ 81, after excluding two locations whose remotely sense data was defective and nine that were impossible to be georeferenced, among which were the two placed outside the housing premises) were correlated with Landsat 8 (LST, NDVI, EVI) and Sentinel 2 (NDVI and EVI) remotely sensed data. We employed a linear regression described by the following equation: where Y i is the mean number of Ae. aegypti eggs, through the study period, at ovitrap location i, Cov(j) indicates the mean of remotely sensed environmental covariates, which are related to the mean number of Ae.
aegypti eggs according to parameters β j , μ is the intercept, and ε i is the error, which is assumed to be normal, independent and identically distributed (Faraway, 2004). Given that observations used to estimate the mean number of Ae. aegypti eggs at the ovitrap locations were not homogeneous during the study period, as a few traps were vandalized or contents accidently removed by household residents, we used error weights proportional to the number of trap-weeks used to estimate the means to accurately represent sampling effort during parameter estimation (Faraway, 2004). We then proceeded with model selection based on the minimization of the AIC, following the same steps already described for the time series models. We tested the best model for nonlinearities in the association with the selected covariates, and if any of the associations was non-linear, we proceeded with fitting a generalized additive model (Venables & Ripley, 2002) that used a smoothed function for the covariates whose association with the mean number of Ae. aegypti eggs was nonlinear. For the generalized additive model, we also used error weights proportional to the sampling effort at each ovitrap location.
Finally, we performed a Moran's I test on the best model residuals, a test with a null hypothesis that a variable is spatially independent (Brunsdon & Comber, 2015). For the Moran's I test we generated a spatial weights matrix linking ovitrap sampling locations within a distance of 145 m, the largest minimum distance between any two ovitrap sampling locations. Once neighbors were identified, weights were made proportional to the number of neighbors for each sampling location (Brunsdon & Comber, 2015).

Data summary statistics
A total of 381 cases of Ae. aegypti-borne arbovirus cases were recorded in Puntarenas county during the ovitrap monitoring period, i.e. from epidemiological week 17 of 2017 to epidemiological week 44 of 2018 (Fig. 4A). Of these, 11 were Zika infections, 19 were chikungunya infections and the remaining 351 cases were due to dengue infections (Fig. 4A). During the study period we counted a total of 291,369 Ae. aegypti eggs, over a total sampling effort of 7,360 ovitrap weeks. A total of 278 times (3.78% of the total sampling effort), we were unable to collect the paper towels with eggs from ovitraps due to diverse reasons, mainly accidental disturbances from residents, and very few instances of vandalism (only 5 ovitraps needed to be replaced during the study duration). Ovitraps had no eggs 39.56% of the times they were sampled, i.e. 2,912 times, with a weekly average (AESD) of ovitraps without eggs of 36.40 AE 10.91, weekly ranging from 14 to 62 ovitraps without eggs across the studied period. The average number of eggs (AESD) per ovitrap/week was 41.14 AE 75.94 with weekly trap counts ranging from 0 to 1,000 eggs (Fig. 4B). Raw data for individual ovitraps are presented in Supplementary Fig. S3. In general, the mean number of eggs for ovitraps was larger than the median, indicating oviposition was clustered at few ovitrap locations during the study period (Fig. 4B). The SD (Fig. 4C) and kurtosis (Fig. 4D) of eggs/ovitrap/week had regular peaks during the study period, which did not seem to coincide with peaks in arbovirus transmission (Fig. 4A) or egg counts (Fig. 4B).

Comparison of land surface temperature (LST) time series estimates with MODIS and Landsat 8 images and remotely sensed data seasonality
MODIS and Landsat 8-based Land Surface Temperature showed a clear seasonal pattern where temperature peaks were observed during the dry season of Puntarenas, weeks 1 to 16 of 2018 (Fig. 5A, raw MODIS and Landsat 8 data presented, respectively, in Supplementary Figs. S2A and S2B). These time series had a low CV (around 7%) for both MODIS (Supplementary Fig. S2C) and Landsat 8 ( Supplementary Fig. S2D), meaning that pixel had a variability of up to 0.7 C for each 10 C of mean temperature in the areas over which the MODIS (Puntarenas city) and Landsat 8 (Puntarenas peninsula) time series were estimated. The estimates in the time series also had a very low sensitivity to the number of pixels considered for area estimations. Estimates, on average, only decreased on the order of 10 À5 for MODIS ( Supplementary Fig. S2E) and 10 À6 for Landsat 8 (Supplementary Fig. S2F) for each pixel not considered on the area estimations. Interestingly, both time series, land surface temperature from MODIS and Landsat 8, were highly correlated (Pearson's b r ¼ 0.734) despite the difference in resolution and areas over where these temperature time series were estimated. The MODIS-based temperature SD time series ( Supplementary Fig. S4A) had a seasonal pattern different to its mean value, as the SD peaked during the dry season of Puntarenas. However, no clear patterns were observed for the MODISbased temperature kurtosis time series (Supplementary Fig. S4B) which seemed to increase randomly, becoming more leptokurtic, i.e. with high Seasonality in NDVI (Fig. 5B) and EVI (Fig. 5C) were opposite to the seasonality observed in temperature, as both vegetation indices decreased when temperature reached peak values. For the study period, MODIS-based NDVI values were higher than those observed for Landsat 8 and Sentinel 2 (Fig. 5B). Meanwhile, MODIS-and Landsat 8-based EVI values were very similar, and higher than values generated with Sentinel 2 images (Fig. 5C). Raw data for MODIS NDVI (Supplementary Fig. S5A) and EVI ( Supplementary Fig. S5B), Landsat 8 NDVI ( Supplementary  Fig. S5C) and EVI ( Supplementary Fig. S5D), Sentinel 2 NDVI (Supplementary Fig. S5E) and EVI ( Supplementary Fig. S5F) also showed a marked seasonal pattern, where both indices decreased during weeks 1 to 16 of 2018, suggesting the seasonality was not an artifact of smoothing the time series.

Time series analysis
The inspection of the ACFs (Fig. 6A) and PACFs (Fig. 6B) suggested that Ae aegypti-borne arboviral cases and Ae. aegypti median egg counts followed first order seasonal autoregressive processes with a 3-week period. This means that observations, for both cases and egg counts, at any given time were associated with observations with a lag of one and three weeks. Using this information for these two variables, a null model with the following general form was fitted: This was then used to pre-whiten the covariate time series in order to select lags at which MODIS-based (Fig. 6C) and Landsat 8-based (Fig. 6D) variables were associated with arbovirus cases according to their CCFs. We used this information to build a full time series model for arboviral cases, which after a process of model selection by backward elimination (Table 1) resulted in a first-order seasonal autoregressive model. The selected model showed that number of cases had a nonlinear response to changes in MODIS-based temperature SD increasing 1 case for each 2 units of SD increase with a 4-week lag, while decreasing almost 1 case for each 2 units of MODIS-based SD increase with a 3-and 4-week lag and decreasing up to five cases for each Landsat 8 temperature SD unit increase with a 1-week lag. Arboviral cases increased around 4 cases for each 0.1 units of Landsat 8-based NDVI increase during the same week cases were recorded (Table 2). Arboviral cases were not significantly associated with the median, SD or kurtosis of Ae. aegypti egg counts ( Supplementary Fig. S6A) nor with Sentinel 2-based NDVI and EVI ( Supplementary Fig. S6B).
The CCFs of median number of Ae. aegypti eggs per ovitrap were not significantly associated with Sentinel 2-based NDVI and EVI (Supplementary Fig. S6C). In contrast, the median number of Ae. aegypti eggs per ovitrap was significantly associated with several MODIS-based (Fig. 6E) and Landsat 8-based (Fig. 6F) environmental variables presented in the full model in Table 3. After the process of model selection, the seasonal component was not significant (Table 3) and the median of eggs increased by approximately three units for each unit of MODIS-based temperature SD, while decreasing approximately 7 units for each 0.1 increase of Landsat 8-based NDVI, both variables having their impact on median egg number with a lag of 6 weeks ( Table 4). The ACF (Supplementary Fig. S7A), PACF (Supplementary Fig. S7B) and CCFs with MODIS ( Supplementary Fig. S7C) and Landsat 8 ( Supplementary  Fig. S7D) of the mean number of Ae. aegypti per ovitrap were similar to the ones observed for the median number of Ae. aegypti egg counts per ovitrap.

Synchrony
Aedes aegypti oviposition synchrony patterns are presented in Fig. 7A. Overall, oviposition synchrony was very low, with a regional synchrony of 0.04, which although statistically significant was very close to 0, i.e. a total lack of synchrony. The observed synchrony pattern is also Fig. 5 Remotely sensed environmental time series. A Weekly mean land surface temperature derived from satellite images. B Weekly NDVI derived from satellite images. C Weekly EVI derived from satellite images. Temperature data originate from daily MODIS images, and Landsat 8 images collected every 8 days. NDVI and EVI time series also included Sentinel 2 data, which were collected every 3.4 AE 2.0 days (mean AE SD). For guidance, please refer to the inset legend of panel A. Raw and weekly smoothed data used to estimate these time series are presented in Supplementary Fig. S2 for MODIS and Landsat 8 land surface temperature, and raw data for the NDVI and EVI, derived from MODIS, Landsat 8 and Sentinel 2 are provided in Supplementary Fig. S5 suggestive of a travelling wave, as suggested by the 7 point estimates shown in Fig. 7A, where synchrony changes from positive to negative, around 400 m (where synchrony is significantly different from the regional estimate, as indicated by the black dot in Fig. 7A), suggesting patterns observed in Fig. 4B reflect a heterogeneous, or with dynamic spatiotemporal clusters, oviposition through time and space, also suggesting dynamic changes in the entomological exposure to Ae. aegypti in Puntarenas city. Landsat 8-based temperature (Fig. 7B) NDVI (Fig. 7C) and EVI (Fig. 7D) were significantly more synchronous than Ae. aegypti oviposition, with synchrony below 100 m being significantly higher than the regional values, suggesting these variables were not associated with Ae. aegypti oviposition in Puntarenas peninsula. Sentinel 2-based NDVI (Fig. 7E) synchrony significantly decreased at distances around 170 m when compared with the regional synchrony, while Sentinel 2-based EVI (Fig. 7F) had a pattern similar to the one observed for Landsat 8-based environmental variables with synchrony below 100 m being significantly higher than the regional values. Sentinel 2-based vegetation indices had low regional synchrony values similar to the ones observed for oviposition, suggesting oviposition synchrony might be driven by finely-grained environmental phenomena. For example, NDVI (Fig. 7E) decreased below the regional synchrony at distances of 170 m, the pattern most similar to what was observed for oviposition synchrony (Fig. 7A).

Spatial models
Selection for the best spatial model (Table 5) allowed us to identify Landsat 8-based EVI kurtosis and temperature as variables significantly associated with the mean number of Ae. aegypti eggs per ovitrap. The relationship with Landsat 8-based EVI kurtosis (LEVIK) was linear, with the number increasing by 3 eggs for each unit of kurtosis increase ( Table 6), suggesting that leptokurtic vegetation changes were associated Fig. 6 Correlation functions. A Autocorrelation (ACF) and B partial autocorrelation function (PACF) for arboviral cases and median number, per ovitrap, of Aedes aegypti egg counts weekly time series. Lines are identified in the inset legend of panel B. C Cross-correlation function (CCF) between weekly arboviral cases and MODIS-based environmental variables, including land surface temperature, NDVI and EVI mean values, and also the standard deviation (SD) and kurtosis (K) of land surface temperature. D CCF between weekly arboviral cases and Landsat 8-based environmental variables, including land surface temperature, NDVI and EVI mean values, and also the SD and kurtosis of land surface temperature. E CCF between median Ae. aegypti egg counts and MODIS-based environmental variables, including land surface temperature, NDVI and EVI mean values, and also the SD and kurtosis (K) of land surface temperature. F CCF between median Ae. aegypti egg counts and Landsat 8-based environmental variables, including land surface temperature, NDVI and EVI mean values, and also the SD and kurtosis (K) of land surface temperature. In all panels the dot-dashed 95% confidence interval lines (see panel A for reference) indicate that correlations within the area are expected by chance. This means that peaks outside the 95% confidence interval lines are the time lags at which the number of arboviral cases (or median egg counts) was associated with itself in the ACF and PACF plots, or with the covariates in the CCF plots with the oviposition of more eggs. However, the association with temperature (T) was nonlinear and we fitted a Gaussian generalized additive model described by the following equation: where s( ) denotes a smoothed function and parameters and assumptions were presented when defining Equation (4). The fitted surface is shown in Fig. 8, where it can be observed that the mean number of eggs increases with both T and LEVIK. The model surface shows that number of eggs increase at a rate of approximately 10 additional eggs for each C above 36 C, with temperature having a minimal impact on the mean number of Ae. aegypti eggs when below 36 C. Finally, lack of significance for the Moran's I statistic (Table 6) supports that assumptions of spatial independence for the error were met, and therefore inferences for the model presented in Table 6 are valid.

Discussion
Our data analysis suggests that ovitrap counts are not useful to predict Ae. aegypti-borne arboviral cases at the studied site. This is interesting, given the general assumption that mosquito, and more generally vector abundance measurements, are useful to understand the risk for pathogen transmission (Smith et al., 2014), and evidence about associations between vector abundance and transmission at different spatial and temporal scales for dengue (Barrera et al., 2011;Cromwell et al., 2017;Parra et al., 2018;Li et al., 2019) and other vector-borne diseases (Chaves et al., , 2014aPoh et al., 2019). Among the diverse factors that may have influenced our results we think some might be related with Ae. aegypti ecology, as it is related to pathogen transmission, and some might be related to the predictability of events that are rare, or that occur at different spatial scales (Levin, 1992;Levins et al., 1994). To start, one of the reasons why egg counts might not be informative about disease transmission is related to the highly likely nonlinear association between egg abundance and adult mosquito abundance. The latter stage is responsible for pathogen transmission, but adult abundance might actually decrease with egg counts because of density-dependence, especially if the relationship is similar to what has been observed between adult and larval abundance (Wilson et al., 1990;Barrera et al., 2011;Chaves et al., 2012). Similarly, the productivity of Ae. aegypti changes across diverse container types (Schneider et al., 2004) and how adult mosquito productivity, i.e. the number of emerging adult mosquitoes, might ultimately reflect trade-offs in oviposition habitat selection by Ae.  (3)

431.29
Notes: AIC indicates the Akaike information criterion for each model. AIC is minimized by the best model, which is presented in bold. Lag indicates the time lag (in epidemiological weeks) for the correlation between arbovirus case number and the environmental variables considered. Abbreviation: SE, standard error. a Statistically significant (P <0.05).  (3)  aegypti Wong et al., 2012). The hypothesis of nonlinear abundance changes through the mosquito life-cycle has been supported by data fits to mathematical models looking at the abundance of Ae. aegypti (Chaves et al., 2012;Lana et al., 2014). More generally, ontogenetic niche shifts, such as those of mosquitoes that alternate between being aquatic and terrestrial organisms during their life-cycle, filter signals of density dependence that might be useful to predict population dynamics (Ripa & Ranta, 2007;Ranta et al., 2006;Ratikainen et al., 2008). Thus, not surprisingly, our results follow a pattern similar to what has been observed in data from other Ae. aegypti aquatic stages, where entomological indices based on the presence, or abundance, of either larvae or pupae, rarely correlate with dengue transmission (Bowman et al., 2014). Mosquito eggs are not conspicuous life stages, which increases the ability to bias their counts. For example, accidently losing eggs during the manipulation of samples, might be a factor that could bias abundance estimates (Lenhart et al., 2005). The accidental loss of eggs is likely to occur when a large number of ovitraps are monitored by a mosquito surveillance crew with multiple responsibilities related to the management of vector-borne disease transmission (CENAPRECE, 2015;Tovar--Zamora et al., 2019). There is also the possibility that egg counts from ovitraps deployed over the larger area of Puntarenas could have been a better predictor for arboviral cases. However, within the Puntarenas peninsula it was clear that egg counts are highly variable at finely-grained spatial scales, and previous studies have shown this area to be the one with most persistent mosquito infestations (Marín Rodríguez & Díaz Ríos, 2012).
Nevertheless, the fact that temperature SD, a weather variable, had a significant impact on egg counts further supports that ovitraps are sensitive devices to estimate Ae. aegypti adult abundance, and most definitely to detect the presence of this important mosquito vector, as suggested by trials comparing multiple trap types (Resende et al., 2013;Codeço et al., 2015). For example, ovitraps had a higher sensitivity to detect Ae. aegypti than larval surveys (Nascimento et al., 2020). In this sense, the patterns observed in Costa Rica are similar to what has been observed for ovitraps in locations as diverse as Texas (Martin et al., 2019), Puerto Rico (Barrera et al., 2014), Brazil (Lana et al., 2014), Trinidad and Tobago (Chadee, 2009), M exico (Tovar-Zamora et al., 2019 and Argentina (Gimenez et al., 2020). These patterns also echo population dynamics patterns observed for adult mosquitoes, which are also sensitive to weather factors (Barrera et al., 2011;Ng et al., 2018;Scott et al., 2000). The relatively long delay of 6 weeks to see the impact of temperature and EVI might emerge via resonance, a phenomenon where natural populations cycles become amplified through a few generations in the presence of the right environmental conditions (Chaves et al., 2014b). For example, 6 weeks is an exact harmonic of the 3-week  (Nisbet & Gurney, 1976). The best model explaining egg counts included EVI, a remotely sensed vegetation index that is associated with water in the environment (Pettorelli et al., 2005). EVI association with median egg counts was negative, probably reflecting that Ae. aegypti is primarily an artificial container species outside its native range in Africa (Dye, 1984;Wallace et al., 2018), needing water for the creation of its man-made larval habitats (Barrera et al., 1993;Predescu et al., 2006). However, Ae. aegypti larval habitats can be destroyed by an excess of rain, which can lead to a flushing of aquatic populations (Koenraadt & Harrington, 2008), ultimately reducing Ae. aegypti abundance (Scott et al., 2000). Still, the association is likely dynamic, as a previous, purely spatial study using high resolution satellite images for Puntarenas did not find an association between vegetation and Ae. aegypti abundance, which was more associated with the built urban environment (Fuller et al., 2010).
The travelling wave pattern of synchrony, where synchrony decreases with distance, and changes its sign before increasing back to its regional value (Ranta et al., 2006), to the best of our knowledge has been only reported for adults of Armigeres subalbatus (Coquillet), a common urban Aedini in Asia (Chaves, 2017b). This is an interesting pattern since it suggests that as much as human movement is an important factor for dengue transmission (Stoddard et al., 2009), vector movement might play a crucial role on Ae. aegypti-borne arboviral transmission. For example, vector control and Ae. aegypti dispersal ecology are prone to be asynchronous (Juarez et al., 2020). This asynchrony might ensure Ae. aegypti persistence in urban landscapes, as has been reported for both adult and larval Ae albopictus in Japan (Chaves, 2017b;Chaves et al., 2020b). Interestingly, although the pattern was not exactly the same, the synchrony of Sentinel 2-based NDVI also showed a significant decrease in synchrony at a shorter distance (170 m) than the one observed for Ae. aegypti eggs at 400 m. This pattern is suggestive of resonance in a two-dimensional expansion wave, a phenomenon, that to the best of our knowledge, has not been described in mosquitoes or other animal populations, but that has been described in neuron populations with hybrid synapses (Sun et al., 2013). Thus, the traveling wave pattern in Ae. aegypti egg count synchrony might emerge because of constrains in Ae. aegypti dispersal that might be conditioned, for example, by the availability of oviposition habitats (Edman et al., 1998), which in turn is linked with weather patterns (Juarez et al., 2020). This might reflect a process similar to the way neuron synapses are regulated in hybrid networks, those undergoing both chemical and electrical synapses (Sun et al., 2013). Interestingly, the evidence suggesting the emergence of this spatial pattern is further supported by the positive impact of high kurtosis, or leptokurtic, EVI on Ae. aegypti mean egg counts at the study site. This observation suggests that Ae. aegypti populations are likely to prefer less stable environments, as observed for treehole Aedini mosquitoes (Sota et al., 1994). This result also implies a nonlinear response to weather conditions, as mean egg numbers spatially increased with temperature in a fashion similar to what has been observed, temporally, for adult Ae. aegypti populations elsewhere (Chaves et al., 2012). For example, in Thailand, Ae. aegypti outbreaks, i.e. sudden changes in mosquito abundance, have been linked with the canalization of high temperatures into life history traits (Chaves et al., 2012(Chaves et al., , 2014b. More specifically, models have shown that prolonged high temperatures can reduce density-dependent mortality while simultaneously increasing fecundity and mosquito productivity (Chaves et al., 2014b).
Like mosquito egg counts, Ae. aegypti-borne arboviral case counts were associated with temperature SD, at multiple spatial scales. This result is in accordance with global observations for dengue, chikungunya and Zika (Wallace et al., 2018) where temperature has been shown as an important factor to explain time series of dengue cases across the globe (Stoddard et al., 2014;Chuang et al., 2017;Siraj et al., 2017;Oidtman et al., 2019). However, at the reduced spatial scale of this study we also found an impact for NDVI, a factor that has been reported as significant for dengue transmission dynamics at a country-wide scale in Costa Rica  and Vietnam (Nguyen et al., 2020). This result might be related to the fact that vegetation growth is associated with several weather variables, like rainfall and temperature (Pettorelli et al., 2005). Finally, the fact that both egg and human cases were associated with temperature SD highlights the dependence of the system not only on average environmental conditions, but also on higher-order moments of environmental variability. This could emerge from Schmalhausen's law, the biological principle stating that organisms are sensitive to both the mean and the variability of the environment (Chaves, 2017a). Overall, our results also highlight the value of environmental information derived from satellite images for the surveillance of Ae. aegypti and the arboviruses it transmits, as documented in Southeast Asia (Nguyen et al., 2020). This pattern has also been observed for Culex pipiens L. and West Nile virus transmission in the USA Poh et al., 2019); and Anopheles albimanus (Wiedemann) abundance and malaria infection in Mesoamerica (Hurtado et al., 2018a,b;Rejmankova et al., 1993;Rigg et al., 2019). Thus, our results are particularly important in the era of big data and unlimited computational  Full model: LST, LNDVI, LEVI, SNDVI, SEVI, LSTSD, LNDVISD,  LEVISD, SNDVISD, SEVISD, LSTK, LNDVIK, LEVIK, SNDVIK, SEVIK   566.91   LST, LNDVI, LEVI, SNDVI, SEVI, LSTSD, LNDVISD, LEVISD Notes: AIC indicates the Akaike information criterion for each model. AIC is minimized by the best model, which is presented in bold. The following covariates were considered: Landsat 8-based temperature (LST); Landsat 8-based NDVI (LNDVI); Landsat 8-based EVI (LEVI); Sentinel 2-based NDVI (SNDVI); Sentinel 2-based EVI (SEVI); and also the SD and kurtosis of each of those variables which are indicated by adding, respectively, SD or K as suffix to the environmental variables. power. Nowadays, satellite data can be streamlined into algorithms for vector-borne disease transmission forecasting and dynamic risk assessments. Finally, ovitrap-based entomological surveillance might also help to guide vector control activities in sites with spatially heterogeneous and asynchronous changes in vector abundance, like these observed for Ae. aegypti in the present study. Major advantages of ovitrap use for Ae. aegypti surveillance include their low cost and ease for systematic field deployment and monitoring (CENAPRECE, 2015;Chaverri et al., 2018;Romero et al., 2019) and the potential involvement of local residents in vector surveillance efforts (Hamer et al., 2018;Tarter et al., 2019;Sousa et al., 2020). This study was carried out during a low transmission period for dengue in Costa Rica and the rest of the Americas in 2017 and 2018 (Perez et al., 2019). This low transmission prevented the use of mosquito control during our study, because mosquito control with pesticides in Costa Rica is only used for human disease outbreak management, and based on human case numbers (Vigilancia de la Salud, 2013). However, the Costa Rican national programme for integrated vector control management is currently evaluating the use of ovitraps to guide vector control in areas prone to dengue epidemics.

Ethical approval
This study was carried out in accordance with Article 7 from Law 9234 for biomedical research which grants the Epidemic Surveillance Division (Vigilancia de la Salud) of Costa Ricaʼs Ministry of Health (Ministerio de Salud) the ability to perform activities with the dual goal of surveillance and research, which are exempted from the approval of an Institutional Research Board, as these efforts are deemed essential for health policy planning and decision making, and do not release individually identifiable data.

Data availability
Data are available upon reasonable request. The request needs to be accompanied by an appropriate ethical clearance in accordance with law 9234 for Biomedical Research in Costa Rica.

Declaration of competing interests
The authors confirm that there are no known conflicts of interest associated with this publication and that there has been no significant financial support for this work that could have influenced its outcome.

Funding
This study was funded by Costa Rica's Ministry of Health.
CRediT author statement JAVC, EM, JMGA, GD, LFC and RMR designed the study. JAVC and RMR coordinated field sampling. GD, LMR, LFC, JMGA, EM, georeferenced ovitrap locations. GD and CAA verified all entomological data. LFC performed all tasks related with the remote sensing data processing and analysis. LFC made all maps. JAVC prepared Fig. 2. MRR, CAA and LMR compiled arboviral case data. LFC and GD analyzed the data. LFC, GD, MRR and LMR wrote the original draft of the manuscript. All authors read, edited and approved the final manuscript.
Puntarenas that allowed us to place ovitraps in their households. This study benefited from discussions at PAHO sponsored entomo-virological surveillance meetings and from the generosity of the Mexican Government that donated the ovitraps and provided training for ovitrap use through the CENAPRECE. Finally, we thank Dr Sofia L. Ermida, Instituto Português do Mar e da Atmosfera (Lisbon, Portugal), for her advice to use scripts for processing Landsat 8 images using the Google Earth Engine, and Dr Mariel Dalmi Friberg Aponte, NASA Goddard Spatial Flight Center (Greenbelt, MD, USA), for her comments about satellite imagery and remote sensing.