Modelling distributions of Aedes aegypti and Aedes albopictus using climate, host density and interspecies competition

Florida faces the challenge of repeated introduction and autochthonous transmission of arboviruses transmitted by Aedes aegypti and Aedes albopictus. Empirically-based predictive models of the spatial distribution of these species would aid surveillance and vector control efforts. To predict the occurrence and abundance of these species, we fit a mixed-effects zero-inflated negative binomial regression to a mosquito surveillance dataset with records from more than 200,000 trap days, representative of 53% of the land area and ranging from 2004 to 2018 in Florida. We found an asymmetrical competitive interaction between adult populations of Aedes aegypti and Aedes albopictus for the sampled sites. Wind speed was negatively associated with the occurrence and abundance of both vectors. Our model predictions show high accuracy (72.9% to 94.5%) in validation tests leaving out a random 10% subset of sites and data since 2017, suggesting a potential for predicting the distribution of the two Aedes vectors.


Introduction
Aedes mosquitoes, in particular, Aedes aegypti (Linnaeus) and Aedes albopictus (Skuse), are the primary vectors of multiple arboviruses including dengue virus (DENV), Zika virus (ZIKV), yellow fever virus, and chikungunya virus (CHIKV) [1][2][3]. The incidence of these viruses in humans is driven, in part, by the close overlapping habitats of humans and these vectors [4]. In the absence of effective vaccines, reducing contact between mosquitoes and humans through targeted mosquito control is regarded as the most effective approach to reducing the risk of mosquito-borne arbovirus transmission. There have been several efforts to create large-scale estimates of the spatial presence and abundance of these vectors using a variety of collection methods and data from literature reports and entomological surveys of mosquito occurrence [5,6]. Global maps have been generated using climate and socio-economic variables, relying on a strong dependence of mosquito populations to temperature and rainfall [7][8][9]. These efforts have uncertainty associated with publication bias and variability of collection methods. Large-scale data collected by standardized surveillance methods could improve the certainty and precision of occurrence and abundance maps.
Florida has suffered from the introduction and autochthonous transmission of DENV [10,11], CHIKV [12] and ZIKV [13,14] and remains at high risk of transmission due to repeated pathogen introductions, high densities of Ae. aegypti and Ae. albopictus [7] and favourable meteorological conditions [13,15]. Studies have shown a positive relationship between human Zika and dengue cases and larger Ae. aegypti populations in urban areas [13,16]. Therefore, characterizing the population size of the two Aedes species over time and space could aid in examining the risk of local arbovirus transmission and spread in Florida and inform more effective and targeted mosquito control efforts.
Although coexistence of the two Aedes vectors is reported [17], declining populations and displaced habitats of Ae. aegypti have been observed in several places, including Florida [3,[18][19][20]. In particular, the habitats of Ae. aegypti were restricted to urban areas while those of Ae. albopictus were found to increase in suburban and rural areas in Florida [21]. The proposed mechanisms for the displacement of Ae. aegypti include species interactions such as the superiority of Ae. albopictus to compete for resources at the larval stage and asymmetric sterilization at the adult stage after interspecific mating, which favours Ae. albopictus [1,3,22]. Previous studies modelled the current spatial distribution of Ae. aegypti and Ae. albopictus by applying boosted regression trees to a comprehensive global database of Aedes occurrence [5,7] and characterized the spatial and temporal abundance of the two Aedes species in a local southern Florida county [23][24][25][26][27][28]. In this study, we build on these previous findings by incorporating longitudinal data collected from a standardized format, providing information on both occurrence and absence with the temporal component. Additionally, a recent systematic review [29] reported inconsistent findings on the associations between the species interactions between Ae. aegypti and Ae. albopictus and meteorological factors, we therefore consider these factors in our model.
The objective of this study was to simultaneously characterize the occurrence and abundance of the Ae. aegypti and Ae. albopictus mosquitoes using routine mosquito surveillance data in Florida. To estimate if mosquitoes were present or not and, if present, the number of adults in each trap location, a mixed-effects zero-inflated negative binomial (ZINB) regression was performed. Various predictors were examined, like climate and human population density covariates, and their potential impact on Ae. aegypti and Ae. albopictus spatial and temporal abundance. In order to evaluate to what extent the model can provide accurate predictions we assessed the performance of the models by validating them against independent data withheld from the model fitting process, especially without the benefits of pre-existing knowledge on abundance and localized spatial variations.

Mosquito surveillance data
Statewide surveillance data on 16 Aedes species were obtained by networking with Florida's mosquito control districts, Clarke Scientific, the Florida Department of Agriculture Consumer Services, and the Florida Department of Health. Each control district is required to trap mosquitoes prior to conducting their control efforts by Florida Statutes 388 and 482. The traps were placed to acquire a representative sampling of the district including baseline traps placed in the same location annually, at risk areas due to environmental factors like increased standing water, locations within areas of known arbovirus transmission, and frequent areas of complaint. Information collected from these traps includes the species-identified count of the trapped adult mosquitoes (total number were calculated where male and female were recorded seperately), date and duration of collection, type of trap (i.e. BG-Sentinels, light traps and other types; details in S1 Table) used, and coordinates of the trap sites. The collected mosquitoes were identified to species level according to standardized mosquito keys [30]. For missing data, the duration of collection was assumed to be one day, according to the common trapping practices, and coordinates were extracted from Google Maps based on the address of the site. The full dataset was aggregated to include data on adult Ae. aegypti and Ae. albopictus, two vectors that transmit arboviruses, on a weekly basis. The longitudinal training dataset for the ZINB regression model was extracted from the full dataset and included only data collected from sites with at least four consecutive weeks of surveillance and no missing explanatory variables (S1 Fig).

Abiotic variables
To examine the potential effects of meteorological factors on the trap rate of the two Aedes species, temperature (˚C), wind speed (meter per second) and relative humidity (%) were included in the model. We obtained the daily meteorological data for Florida from the regional data of NASA Prediction of Worldwide Energy Resources [31] and applied the inverse distance weighting method [32] to interpolate the daily weather raster of Florida with a 5 km×5 km resolution. To examine the effect of data source, interpolation and spatial resolution on our results, we also conducted sensitivity analyses using meteorological data from National Oceanic and Atmospheric Administration (NOAA, 5 km×5 km) and Daily Surface Weather and Climatological Summaries (Daymet, 1 km×1 km) [33,34]. Data from NOAA was available for stations that were disproportionately distributed in densely populated places, and we interpolated these data using the above-mentioned inverse distance weighting. We downloaded the Daymet data using the coordinates and "daymetr" package [35] without interpolation. The weekly average of weather conditions was calculated as the mean of the weather conditions on the days the traps were collected. To account for the collinearity of the maximum and minimum temperature, we used the residuals of the linear regression of maximum temperature on minimum temperature as a proxy of the maximum temperature in the model, which was calculated as ΔT max = T max -(α+βT min ), where T max and T min denoting the observed maximum and minimum temperature, respectively, while α and β were estimated from the linear regression. We used human population density as a proxy for urbanization. We used data on population density obtained from the Center for International Earth Science Information Network with a 5 km×5 km resolution for the year 2015 [36]. If the value was missing for a site, we extracted the corresponding environmental variables based on its coordinate and used the average drawn from a 5km buffer around the site.

Statistical methods
We applied a ZINB regression model to the weekly abundance of Ae. aegypti and Ae. albopictus from the longitudinal training dataset, respectively, to account for the excessive zeros in the abundance data and the over-dispersed count of trapped mosquitoes, simultaneously. The ZINB model comprises a binary component (corresponding to the absence/presence of mosquitoes), and a negative binomial component (corresponding to the abundance of mosquitoes). The estimates from the binary component (presented as odds ratio, OR) and the negative binomial component (presented as incidence rate ratio, IRR) represent the associations between the covariates and the occurrence and abundance of these Aedes vectors, respectively. The potential factors included in the ZINB model for both species are: the previous abundance of Ae. aegypti and Ae. albopictus up to three weeks prior (in log-scale), weekly site-specific meteorological factors (i.e. wind speed, minimum temperature, the residual of maximum and relative humidity), human population density (in log-scale) and type of mosquito traps (i.e. BG-Sentinels, light traps and other types). We examined the potential interaction between Ae. aegypti and Ae. albopictus by examining the relationship between the current abundance of one species with the previous abundance of the other species. We used counts of each species detected in recent weeks to predict future weeks. To do this, we only considered records when data was available for four consecutive weeks prior. Trap type was included as an explanatory covariate as each of the traps used has a different effectiveness in trapping each species. We also included the random effects at both site level and county level, which were modelled for both components of the ZINB model simultaneously. The detailed equations used for the ZINB model are provided in the S1 Text. Parameters were estimated by maximizing the likelihood using "glmmTMB" package [37] in R version 3.5.0 (R Foundation for Statistical Computing, Vienna, Austria).

Goodness of fit of the model
We assessed the goodness of fit of the model by comparing the observations with the predictions of occurrence and abundance from the longitudinal dataset. We assessed the spatial pattern by calculating the site-specific mean of differences between observations and predictions. The absence and presence were assigned as 0 and 1 respectively for calculation purposes. Moran's I was calculated to assess the spatial autocorrelation of Aedes distribution [38]. We examined the temporal pattern of the model fitting by assessing the monthly 2.5% and 97.5% percentile of the difference between the predicted and observed abundances for the two Aedes species.

Cross-validations
We also performed cross-validation of the model both spatially and temporally. Prediction of a test set was based on a model fit from a training set and comparing the predicted and observed occurrence and abundance. In the spatial prediction, we randomly selected records from 127 (around 10% of total) sites to be the spatial testing set and used the records from the remainder of the sites as a spatial validation training set (S1 Fig). In the temporal prediction, we used data up to the year of 2017 as the temporal validation training set to predict data after 2017 (S1 Fig). The area under the receiver operating characteristic (AUC) was used to measure the performance of prediction on the mosquito occurrence, using the predicted probability of the occurrence as the predictor. In order to assess the accuracy of predictions on abundance, we divided the observed and predicted trap rate (r, per trap day) into four categories, i.e. r = 0, 0<r<10, 10�r<100 and r�100. We define the predictions on abundance as correct if the predicted abundance category is the same with observed abundance category. The proportion of correct predictions were calculated by dividing the number of correct predictions on abundance by the number of traps where the mosquito was found and predicted to occur (excluding the impact of misclassifications of presence/absence).

Utility of model prediction
In order to assess the model utility in predictions, we first compared the observations of occurrence and abundance from the longitudinal dataset with predictions using models that incorporate 1) both random effects and prior abundance information, 2) random effects only ("no abundance model"), 3) prior abundance information only and 4) none of random effects and prior abundance information. All models accounted for climate factors, human population density and trap types. We evaluated whether the proportion of the predicted presence and abundance categories were consistent with observations. In addition, to test the model's utility for real-time prediction, we used the "no abundance model" fit using longitudinal data to predict an external "no abundance testing dataset" (S1 Fig) and compared the predictions with observations. The no abundance testing dataset consists of the surveillance records in the full dataset that failed on the four consecutive four-week criteria (S1 and S2 Figs). We also assessed the performance of predictions calculated by removing the random effects estimates of the "no abundance model" in order to illustrate the mean trap rate in Florida.

Results
We integrated a full dataset from counties in Florida that contains 180,242 weekly records and representative of around 102,000 km 2 (73% of Florida) between 2004 and 2018 (S1 and S3 Figs). From the full dataset, we extracted a longitudinal training dataset that only included data collected from sites where at least four consecutive weeks of surveillance were available. The longitudinal training dataset included 132,088 weekly records from 1,246 unique sites for Ae. aegypti and Ae. albopictus, covering 33 out of 67 counties from 2004 to 2018. The land area covered by the counties that we have data for represents 53% of the land area in Florida (Table 1 and S1 and S3 Figs and S1 Video). Traps were typically set for one day but a minority of collaborators reported counts from a trap that was set for multiple days (7.4%). Approximately 87.4% and 84.8% of trap episodes reported no adults collected for Ae. aegypti or Ae. albopictus, respectively. The majority (81.4%) of traps used were light traps, and the remaining 7.3% and 11.3% of traps used were BG Sentinel traps or other mosquito traps (Tables 1 and  S1), respectively. A wider distribution range and higher trap rate was reported for Ae. albopictus compared to Ae. aegypti in Florida, and, as expected from previous studies, most Ae. aegypti were reported in central and southern Florida (Figs 1 and S3). Both Ae. aegypti and Ae. albopictus were trapped more often between May to October (S3 Fig and S1 Video).
The median human population density at the locations where the traps were set is 480.8 persons per km 2 (Interquartile range (IQR), 112.5 to 1,165.

Presence and abundance of Aedes
The results from ZINB regression suggested the probability of presence of Ae. aegypti and Ae. albopictus in the current week was positively associated with their respective abundance in the previous week. ( Table 2). The abundance of both Aedes species was more likely to be higher if a higher abundance was reported for its own species (e.g. Incidence rate ratio (IRR) 1.03 and 1.02 for one week prior for the two vectors, respectively) ( Table 2). The abundance of Ae. aegypti was negatively associated with the abundance of Ae. albopictus in the last three weeks (IRR: 0.992, 0.994 and 0.990 for one, two and three weeks earlier, respectively), while the abundance of Ae. albopictus was not associated with the previous abundance of Ae. aegypti (Table 2).
We found both the presence (Odds ratio (OR): 0.98, 95% confidence interval (CI) 0.95 to 1.01 and 0.97, 95% CI, 0.95 to 0.99, respectively) and abundance (IRR: 0.97, 95% CI, 0.96 to 0.99 and 0.97, 95% CI, 0.95 to 0.98, respectively) of Ae. aegypti and Ae. albopictus were negatively associated with the average wind speed of the week, although such association was not significant for the presence of Ae. aegypti. Minimum temperature was positively associated with the occurrence (OR: 1.08, 95% CI, 1.07 to 1.09) for Ae. albopictus and the abundance (IRR: 1.13, 95% CI, 1.12 to 1.14 for Ae. aegypti and 1.09, 95% CI, 1.08 to 1.10 for Ae.  Table). Both the probability of presence (OR: 0.95, 95% CI, 0.94 to 0.97) and abundance (IRR: 0.98, 95% CI, 0.97 to 1.00) of Ae. albopictus were negatively associated with a higher human population density, while the probability of the presence of Ae. aegypti was positively associated with human population density (OR: 1.05, 95% CI, 1.03 to 1.07). The probability of presence and abundance of Ae. agypti and Ae. albopictus was lower when using light traps and other trap types, compared to using BG traps ( Table 2). We also found substantial heterogeneities of presence and abundance of these two Aedes species across trap sites and counties ( Table 2). The heterogeneity was greater at the county level (random effects (RE): 12.27 for Ae. aegypti and 6.59 for Ae. albopictus) compared to the site level for both species.

Goodness of fit of the model
We compared the predictions from the main ZINB model with observed presence and abundance from the longitudinal training datasets (Figs 2 and S6 and S2 Video). Overall, our model fits well with both the occurrence and abundance estimates for Ae. aegypti and Ae. albopictus (Figs 2 and 3). Places where inconsistent predictions and observations on presence/absence were observed are places with higher trap rates of the two species (Fig 2). We observed that 91.1% (95% CI, 91.0% to 91.3%) and 84.9% (95% CI, 84.7% to 85.1%) of the predicted presence/absence was consistent with the observations of Ae. aegypti and Ae. albopictus, respectively (S4 Table). Similarly, 78.7% (Ae. aegypti, 95% CI, 77.9% to 79.4%) and 84.9% (Ae. albopictus, 95% CI, 84.3% to 85.5%) of the predicted abundance was consistent with the observations among traps where the mosquito was captured and predicted to be present (S4 Table). The values of Moran's I are 0.47 (p < 0.01) and 0.08 (p = 0.02) for Ae. aegypti and Ae. albopictus, respectively, and is -0.03 (p = 0.81) for Ae. aegypti after removing data from Miami-Dade (S5 Table). Temporal differences and average trap rates were relatively higher between May and September for both species (S6 Fig).

Cross-validations
We performed the cross-validations both spatially (i.e. holding out 10% of sites) and temporally (i.e. holding out data after 2017). For spatial ( Fig 3A) and temporal (Fig 3D) cross-validations, the model predictions are highly consistent with the observed presence of both Ae.  Table).

Utility of model prediction
In order to assess the minimum information needed for the model to provide accurate predictions, we fit the different models incorporating various combinations of random effects and prior abundance information. We found greater reductions in goodness of fit when only random effects were removed than when only removing prior abundance information (Fig 4 and  S4 Table). The goodness of fit was reduced by 15.0% (Ae. aegypti) and 26.1% (Ae. albopictus) for the model simultaneously when removing random effects and prior abundance compared to the model that included both. Using the model that included random effects but not prior abundance information ("no abundance model" hereafter), we predicted the trap rate of Ae. aegypti and Ae. albopictus at all points in the state in 2018 with BG traps (Figs 5 and S7 and S7 Table). Results in Fig 5 show the predictions incorporating random effects and represent systematic differences in trap rates by county (Fig 5A  and 5B). The predictions from fixed effects, which represent the mean trap rate, can capture more temporal trends than spatial heterogeneity in Florida (Fig 5C and 5D). We further conducted validation of the above predictions with the no abundance testing dataset (S1 Fig), which is an external dataset that was not used for model training. The performance of model incorporating random effects was better (AUC: 0.90 for Ae. aegypti and 0.85 for Ae. albopictus) than predictions of only fixed effects (AUC: 0.81 for Ae. aegypti and 0.75 for Ae. albopictus) (S7 Fig and S7 Table).

Discussion
We built models using more than 132,000 routine mosquito surveillance records from 33 counties in Florida collected from 2004 to 2018 to characterize and predict the occurrence and

PLOS NEGLECTED TROPICAL DISEASES
abundance of Ae. aegypti and Ae. albopictus. Our analysis is set apart by the extensive localized data set that was collated from counties throughout Florida, which strengthens the predictive power of our model. Our model performed well, particularly considering the stochastic nature of mosquito populations, trap efficiency and small-scale trap locations. We modelled random effects across sites and counties to account for unmeasured abiotic factors, inconsistencies, and randomness and found the highest random effect was for the probability of presence at the county level, suggesting great heterogeneity of occurrence across counties possibly down to differences in varying micro-scale environment, surveillance, and domestic mosquito control across counties.
Our results suggest a broad distribution of Ae. albopictus in Florida, while Ae. aegypti was more likely to be found in counties in southern Florida, a pattern similar to reports during the past two decades [17]. Such results were consistent when using different sources of climate data. This is also consistent with previous observations about the declining population of Ae. aegypti after the invasion of Ae. albopictus in the southern United States [21,39]. However, there is some evidence to suggest limited local recoveries of Ae. aegypti in relation to Ae. albopictus, in part, attributable to evolution of resistance to satyrization [1,17,40,41]. Our findings on the positive association between the probability of presence of adult mosquitoes of the two Aedes species suggest their niches have some overlap, particularly in urban areas [42]. This is supported by the observed coexistence of Ae. aegypti and Ae. albopictus in Florida [3,23,39]  and the similar breeding behavior of the two species [43]. The abundance of Ae. aegypti was negatively associated with the previous abundance of Ae. albopictus, with the greatest effective size observed for the abundance of Ae. albopictus during the previous three-week period. These results suggest competitive interactions between the two species. A previous study revealed the breeding preference of Ae. aegypti in habitats without Ae. albopictus [43]. Our findings support the hypothesis that the two Aedes species can coexist but the abundance of adult Ae. aegypti are suppressed due to its failure to outcompete at the larval stage and/or the impact of interspecific mating [3,22,40]. Evolution of resistance to interspecific mating (i.e., satyrization-resistance) in Ae. aegypti populations is likely to promote coexistence [1]. Future control efforts targeting the Aedes species, especially Ae. albopictus, need to consider the risk of resurgence of Ae. aegypti, which has been documented in Brazil [2], and could be possible in Florida considering recent reports of the rapid evolution of satyrization-resistant aegypti [1], coupled with an observed increased in insecticide resistance as compared to Ae. albopictus [44].
We found the presence and abundance of Ae. albopictus are negatively associated with human population density, while the presence of Ae. aegypti was positively associated with the human population density, which matches with reports that anthropophilic Ae. aegypti are more likely to be found in urban areas and Ae. albopictus has wider range of habitats including peri-urban, vegetated and rural areas [3,45], mostly due to its wide range of host preference and a greater adaptation to different climates [39]. Land cover status, which is an important predictor of distribution of these species by other reports [7,46], was not included in our main analysis as it may be associated with the human population density and the vast majority of the Ae. aegypti were collected in developed areas with a large human presence, consistent with other studies [47,48]. The observed positive association between human and Ae. aegypti densities has practical implications for targeted mosquito control because these areas represent the greatest risk for arboviral infections (e.g., dengue [49]).
Elevated abundance of both species between May and October has also been reported. This timeframe corresponds to Florida's rainy season, associated availability of breeding sites, and a warmer temperature [23]. Related to these abiotic factors, the negative association between wind speed and the presence and abundance of both species can be explained by high wind speed hindering the effective trapping of the mosquitoes [50,51]. Traps are more likely to under-catch mosquitoes on windy days. Likewise, mosquito host-seeking activity has been shown to be impacted by higher wind speeds, which is presumably due to the wind's influence on flight distance, pattern, and poor dispersal of the CO 2 plume over both short and long distances.
Furthermore, our results suggest positive associations between the minimum temperature and the observed abundance of adult Ae. aegypti and Ae. albopictus when using the NOAA data. However, inconsistent findings on the association between the residuals of maximum temperature and the abundance of Ae. aegypti was found when using the NASA data. One study suggested higher tolerance of low temperatures in adult Ae. aegypti compared to Ae. albopictus leading to a relatively lower mortality of adult Ae. aegypti at low temperatures and a milder effect of temperature on the presence of Ae. aegypti [52]. Alternatively, another study observed that Ae. albopictus prefer to live in cooler areas in Florida [39]. However, different local adaptations by these Aedes species to climatic changes were reported both in and out of Florida [3,53].To investigate this, we modeled the residuals of maximum temperature to avoid its collinearity with minimum temperature. This might, however, hinder the interpretation of the residuals of maximum temperature due to the non-linear association between temperature and mosquito survival. Further, the same residuals of maximum temperature may impact the mosquito differently conditional on different minimum temperatures. We performed a sensitivity analysis using mean temperature only, and the associations identified in this study seem to be unaffected by the choice of temperature (S8 Table). Despite these discussions pertaining to the relation between mortality of the two Aedes vectors and temperature, seasonality can be used to predict the patterns of presence and abundance of these two Aedes species and the incidence of diseases transmitted by the these mosquito vectors [15,23,54].
We find a negative correlation between relative humidity and the abundance of Ae. aegypti and the presence of Ae. albopictus. These findings support laboratory and field observations showing climate-driven egg mortality, with greater desiccation resistance in Ae. aegypti than Ae. albopictus, and species-specific responses in occupancy of containers with drier conditions favoring Ae. aegypti [55][56][57]. Previous field studies have shown that dry periods are associated with disproportionately greater mortality of Ae. albopictus eggs than Ae. aegypti eggs in Florida [57]. Previous laboratory studies revealed desiccation stress on survival of adult Ae. aegypti and Ae. albopictus with mortality increasing non-linearly with decreasing relative humidity [58][59][60]. The complex relationship between adult survival, relative humidity, and the observed higher relative humidity in Florida could drive the negative association (S3 Fig). In addition, higher relative humidity was usually associated with greater precipitation, which was found to be positively correlated with the abundance of Ae. aegypti, but not the probability of occurrence of the two species in the sensitivity analysis (S2 Table). The effect of precipitation on the abundance of these two Aedes species was considered to be mediated by induced egg hatching in containers upon flooding and promotion of vegetation after raining [23,29]. The larger effect of precipitation on the abundance of Ae. aegypti than of Ae. albopictus could be due to this species' preference for ovipositing in artificial containers, which are prone to have more obvious influence from precipitation compared to vegetation.
The probability and efficacy of capturing Ae. aegypti and Ae. albopictus by a BG-sentinel trap was found to be greater compared to light traps (Table 2), which is consistent with previous findings [61,62]. We performed sensitivity analysis by fitting the model to data collected by BG sentinel traps only or light traps only (S9 Table). We found the robustness of our main results are seemingly unaffected by the spatial distribution of BG sentinel traps (S8 Fig). In addition, we were not able to assess the role of attractants due to limited data available, which are believed to increase the capture efficacy of mosquitoes [63].
Results from the assessment of goodness of fit and cross validation suggested that our model can provide highly accurate predictions on the presence and abundance of Ae. aegypti and Ae. albopictus, especially when the model incorporates the previous abundance of heterospecific and conspecific Aedes species at a trap. Analysis of long-term mosquito surveillance data is challenged by the excessive zero counts, which may be real absence, or absence due to trap failure or adverse environmental conditions. The ZINB regression can model the two scenarios of absence simultaneously. A larger rate of inaccurate predictions was observed during months when trap rates of both mosquito species were higher, which is due to the more dispersed variance of a higher trap rate and the exponential growth of mosquito populations. The relatively higher inconsistent proportion between observed and predicted occurrence in places with higher trap rates (Fig 2) is caused by the higher chance for false absences in places where Aedes can be found. In addition, spatial autocorrelation was found for the model of Ae. aegypti, which was mainly due to the high autocorrelation between observations in Miami-Dade. The estimates and predictions are however not affected by the spatial autocorrelation, as suggested by the model fit to the longitudinal training dataset removing data from Miami-Dade (S9 Fig  and S5 Table).
In order to account for the autocorrelations of the repeated measurements from the same trap and the potential unmeasured confounders that correlated with locations, we incorporated random effects in the model. However, we found a large proportion of spatial variations was explained by the random effects at county level. Despite the fact that the results from the "no abundance model" incorporating fixed effects are only consistent with previous predictions on the two Aedes vectors occurrence and survival probability in Florida [8,64,65], the maps are less heterogeneously predictive compared to models incorporating random effects. This is because both 5km × 5km climate and human population density demonstrate relatively less spatial variation in Florida, while the empirical data suggested great variation in the abundance captured across counties (Figs 1 and S1), which could be partially explained by climate data at a finer spatial resolution (as suggest by results using Daymet data). The county-level random effects could be largely mediated by the pre-existing niches of the two Aedes vectors, which are the result of unmeasured abiotic factors (S3 Table). Such unmeasured abiotic factors could be the systematic differences in mosquito surveillance across counties or other factors that are critical for mosquito survival, such as the micro-scale climate conditions, the density of mosquitoes favoring micro-environment (e.g. water containers and abandoned tires), the density of hosts other than human, the stochasticity of the establishment of habitat and the varying efforts of vector spraying across counties. Ideally, the model performance could be increased after considering these factors, but it is practically challenging to obtain such information.
Our no abundance model, which used random effect but not the prior abundance information, demonstrated good potential to provide real-time predictions on the occurrence and abundance of Aedes, especially in places with long-term mosquito surveillance. Although the performance on predictions (especially for occurrence/presence) for the two Aedes vectors could be affected by lack of information on local variations when extrapolating to outside counties, our model can still provide satisfactory predictions. In addition, the performance of temporal validations could be affected by the out-of-sample predictions as many counties started mosquito surveillance after 2017 (S1 and S3 Figs). The dynamic changes of niches of the mosquitoes may further hinder the distribution forecasting of the two vectors. Incorporating the species invasion process could help to improve modeling at places where rich spatial information is available.
Many efforts have been made to map the distribution of Ae. aegypti and Ae. albopictus at broad regional scales, which were highly dependent on vegetation and meteorological factors [7][8][9]64,65]. Places like Florida are theoretically suitable for Aedes survival and previous predictions using coarser scale climate and host density may have reduced utility for local mosquito control due to the relatively small spatial variations of such predictors at a fine scale. Results from this study suggested greater spatial variations that cannot be explained by climate factors and host population density alone, which calls for more detailed localized data to further aid predictions at fine scale. In addition, our study observes suppression of adult population of Ae. aegypti by Ae. albopictus, highlighting the importance of including species interactions in future mapping work as underscored by recent studies, especially when considering predictions at high spatial resolution [42]. Otherwise, the distribution of Ae. aegypti would likely be overestimated since the two Aedes vectors shared many common abiotic conditions. In addition, integrating standardized longitudinal mosquito surveillance could provide valuable information on absence and abundance, therefore reducing the sampling bias and disproportional weighting caused by presence only data [66]. By creating and making such a dataset available, it also enables the temporal predictions on presence and abundance, rather than a single prediction on occurrence.
There are several limitations to our study. First, our data has relatively more trap episodes during April to November, when the trap rate for these two vectors was often high. The estimated impact of low temperature on the presence and abundance of these two Aedes vectors may therefore be affected. Second, more than half of the records included in the main analysis are from Miami-Dade, St. Johns, Polk and Pinellas counties (S10 Table). We have modelled the random effects across both sites and counties to account for the potential spatial variations of surveillance, which may improve the generalization capability of our conclusions. We were not able to characterize specific details of trap locations or other aspects of the built environment, such as the micro-scale climate and environment. These details could explain more spatial variations in the distribution of the mosquitoes, as suggested by our sensitivity analysis that used a finer spatial resolution and had smaller site-and county-level variations.
Our models demonstrate potential for predicting the occurrence of Ae. aegypti and Ae. albopictus, to better inform targeted mosquito control efforts. Model predictions produced with and without the benefit of recent surveillance data were of high accuracy suggesting that real-time forecasts could be produced with just climate data alone. Our results, however, call for the need for additional local data to explain a large spatial variation in mosquito occurrence and abundance.
Supporting information S1 Text. Description of Mixed-effects zero-inflated negative binomial regression. (DOCX) S1  Table. Odds ratio (OR) and incidence rate ratio (IRR) estimate from mixed-effects zero-inflated negative binomial analysis of covariates of Aedes aegypti using mean temperature. (DOCX) S9 Table. Odds ratio (OR) and incidence rate ratio (IRR) estimate from mixed-effects zero-inflated negative binomial analysis of covariates of Aedes aegypti and Aedes albopictus collected from BG traps. (DOCX) S10 Table. Surveillance data by county.