The global distribution of the arbovirus vectors Aedes aegypti and Ae. albopictus

1The mosquitoes Aedes aegypti [= Stegomyia aegypti] and Aedes albopictus [= Stegomyia albopicta] are vectors of several globally important arboviruses, including dengue virus (DENV), yellow fever virus, and chikungunya virus (CHIKV). The public health impact of DENV and CHIKV has increased dramatically over the last 50 years, with both diseases spreading to new geographic locations and increasing in incidence within their range. The remaining burden of vaccine-preventable yellow fever is similarly likely to be dramatically underestimated. DENV, with a nearly ubiquitous distribution in the tropics and more recently introduced to Europe, is the most prevalent human arboviral infection causing 100 million apparent annual infections world-wide with almost half of the world's population at risk of infection. CHIKV recently received considerable public health attention due to the outbreaks in Réunion in 2005–2006 (225,000 infections), Italy in 2007 (205 infections), and France in 2010 and 2014 (2 and 11 locally transmitted cases, respectively) as well as its recent invasion into the Americas with over 1 million cases recorded to date. Increases in distribution and intensity of transmission are compounded by the lack of commercially available antivirals or vaccines for either disease, although new therapeutics and vaccines are in development. Similarly, while yellow fever infections have been on the decline due to extensive vector control and an effective vaccine developed more than 70 years ago, it still causes a significant disease burden in Africa and South America. Given the public health impact of these diseases and their rapid global spread, understanding the current and future distribution, and determining the geographic limits of transmission and transmission intensity, will enable more efficient planning for disease control. Because these diseases can only persist where their mosquito vectors, Ae. aegypti and Ae. albopictus are present, understanding the distributions of these two species underpins this strategy.

2The global expansion of these arboviruses was preceded by the global spread of their vectors. Ae. aegypti originated in Africa where its ancestral form was a zoophilic treehole mosquito named Ae. aegypti formosus. The domestic form Ae. aegypti is genetically distinct with discrete geographic niches. It was hypothesised that due to harsh conditions coupled with the onset of the slave trade, Ae. aegypti were introduced into the New World from Africa, from where it subsequently spread globally to tropical and sub-tropical regions of the world. Ae. albopictus, originally a zoophilic forest species from Asia, spread to islands in the Indian and Pacific Oceans. During the 1980s it rapidly expanded its range to Europe, the United States and Brazil. Today both Ae. aegypti and Ae. albopictus are present in most Asian cities and large parts of the Americas. Ae. aegypti feed almost exclusively on humans in daylight hours and typically rest indoors. In contrast Ae. albopictus is usually exophagic and bites humans and animals opportunistically but has also been shown to exhibit strongly anthropophilic behavior similar to Ae. aegypti in specific contexts.

3A number of previous studies have mapped the global or regional distributions of Ae. aegypti and Ae. albopictus by focusing on different aspects of their ecology. The majority examined the impacts of climatic conditions, often with an exclusive focus on temperature. and used results from laboratory studies to identify potential limits of establishment in Japan and Asia suggesting a minimum mean temperature in the coldest months of −2°C and −5°C respectively limits their distribution. extended that work by modeling the adult survival of both species under laboratory and field conditions, indicating that Ae. albopictus has higher survival rates than Ae. aegypti, though adults of the latter can tolerate a wider range of temperatures. Applying these results to global temperature data, produced maps indicating areas where the temperature is suitable for these vectors to persist. Whilst temperature is clearly a crucial factor constraining the distribution of the two species, these results alone are not sufficient to discriminate between areas where the species can and cannot persist. Other studies went further using statistical models, predicting the distributions of both species (though particularly Ae. albopictus) using a broader range of climatic variables including precipitation.

4Whilst these studies incorporated several generic climatic factors to predict the current and future distribution of the species, we were able to integrate a bespoke species-specific temperature suitability covariate and account for anthropogenic factors that are known to influence Ae. aegypti and Ae. albopictus distributions. Both species are container-inhabiting but differ in their behaviour and biology so that they occupy different niches. A few local studies showed, however, that local spread of Ae. albopictus and declining Ae. aegypti populations might be linked to inter-species competition and/or non-reciprocal cross-species inseminations. Socio-economic factors affecting the distribution of the Aedes mosquitoes other than the use of containers to store water, include the use of air-conditioning, housing quality, and the rate of urbanisation. In addition to exclusively focusing on meteorological factors in determining the spatial extent of the Aedes mosquitoes, many models used small sets of input occurrence data, which were biased towards particular countries with well-developed surveillance systems, such as, Brazil and Taiwan.

5In this context, we set out to model the global distribution of these two important vector species, compiling the most comprehensive occurrence dataset to date from published literature and national entomological surveys. To overcome previous modelling limitations, a probabilistic species distribution model using Boosted Regression Trees (BRT) was produced for each vector. Our models combine environmental and, for the first time, land-cover variables to predict the global distribution of both species at high spatial resolution. Importantly, the models quantify prediction uncertainty and aim at identifying key contributing factors and inter-species differences in their environmental niches.

6In total, data collection yielded 19,930 and 22,137 spatially unique occurrence records for Ae. aegypti and Ae. albopictus respectively, which were used to train the distribution models. This includes up-to date records from national entomological surveys from Brazil and Taiwan for both species. For Ae. aegypti, >60% of all occurrence records are from Asia and Oceania, 35% are from the Americas and only 575 unique occurrences are available for Africa and Europe (Table 1a). Similarly for Ae. albopictus, most of the occurrences are from Asia (75%), 23% are from the Americas and only 542 records are available from Europe and Africa (Table 1b). For each continent the top 10 countries in terms of occurrences recorded are shown for both species (Table 1). The geographic distribution of the occurrence records is the widest ever recorded with particularly high spatial and temporal resolution in Taiwan and Brazil for both species and in the United States for Ae. albopictus. All occurrence data have been made openly available through an online data repository to ensure consistency and reproducibility.

7Maps showing the predicted global distribution for Ae. aegypti and Ae. albopictus are presented in Figures 1, 2, respectively. The distributions of the two species differ markedly in a number of places. Ae. aegypti is predicted to occur primarily in the tropics and sub-tropics, with concentrations in northern Brazil and southeast Asia including all of India, but with relatively few areas of suitability in Europe (only Spain and Greece) and temperate North America. In Australia, however, Ae. aegypti shows a wider geographic distribution than Ae. albopictus, which is confined to the east coast, largely reflecting the known historic distribution of Ae. aegypti. By contrast, the distribution of Ae. albopictus extends into southern Europe (Figure 3A), northern China, southern Brazil, northern United States (3b), and Japan. Again, this reflects the current and historic distribution of Ae. albopictus and the ability of the species to tolerate lower temperatures.

8In Europe, the predicted potential distribution of Ae. albopictus contains most of the known occurrence points, but suitability is also predicted in Portugal and the west of Spain, and in much of south-eastern Europe and the Balkans, where the species has yet to be reported. Similarly, in China Ae. albopictus has yet to be reported from much of the area predicted to be environmentally suitable. By contrast, in the United States the species has been reported from almost all of the predicted suitable areas, with the exception of a small band of predicted suitability on the western slope of the Sierra Nevada. Due to the relatively sparse reporting from Africa it remains uncertain whether areas predicted to be highly suitable are already infested or have yet to be colonized by the species. Ae. albopictus for example has only been reported from some West African countries (Nigeria, Cameroon, Gabon, the Central African Republic, Congo, Côte d'Ivoire) and Madagascar, and South Africa (as well as some islands in the Indian Ocean). The distribution of Ae. aegypti in Africa seems to be much wider, with reports of species occurrence in over 30 countries.

9For both species, the most important predictor was temperature. Temperature suitability indices had high relative influence statistics for both species; this variable was selected in approximately half of regression tree decisions for Ae. aegypti (54.9%, CI = 53.7–56%) and Ae. albopictus (44.3%, CI = 42.7–45.6%). The full definition of a relative influence statistic is given in the ‘Materials and methods’ section under the heading Predictive performance and relative influence of covariates. Precipitation and vegetation indices made up the remainder of predictors. Urban land cover made very little contribution to either model (Table 2). Model evaluation statistics under cross-validation were high (AUC: 0.87 and 0.9 respectively) for both model ensembles, indicating high predictive performance of the model. Effect plots for each covariate are shown in Figure 1—figure supplement 2. Maps of uncertainty associated with these predictions are presented in Figure 1—figure supplement 3.

10By combining the most comprehensive dataset of occurrence records with an advanced modelling approach and a bespoke set of environmental and land-cover correlates, we have produced contemporary high-resolution probability of occurrence maps for Ae. aegypti and Ae. albopictus, two of the most important disease vectors globally. Dengue and chikungunya, pathogens transmitted by these vectors and rapidly expanding in their distributions, are increasingly prominent in public health agendas and pose significant health threats to humans. In common with previous work to map the global distributions of the dominant vectors of malaria, the maps will improve efforts to understand the spatial epidemiology of associated arboviruses, and to predict how these could change in the future. Specifically, these maps may be used to prioritize surveillance for these vector species and the diseases caused by the viruses they transmit in areas where disease and entomological reporting remains poor. For example, in parts of Asia and Africa where there is a mismatch between predicted environmental suitability and reported occurrences, these maps could be used to determine whether the vector has yet to fill its niche or if it is present but has not been reported due to limited entomological surveillance. They may also be used to identify areas where the species could persist but has yet to be reported, in order to proactively prevent vector establishment.

11The relative contributions of each of the environmental covariates to the global models concur with our theoretical and experimental understanding of each species' biology. Both species' distributions are highly dependent on the limiting factor temperature places on survival of the adult mosquitoes and on the gonotrophic cycle (Table 2). The inclusion of a bespoke temperature suitability index, both in defining the pseudo-absences and as a covariate, allowed us to capture both geographic and temporal variations in the species-specific effects of temperature in a single variable, leading to improved predictive skill of the models. As both Ae. aegypti and Ae. albopictus lay their eggs in small water-filled containers, it is encouraging that precipitation also has a strong influence on the model's predictions. The stronger influence of minimum precipitation for Ae. albopictus than for Ae. aegypti (16.1% vs 9.1%, Table 2) may reflect the former species' preference for non-domestic juvenile habitats, which are solely reliant on filling via precipitation. By contrast, Ae. aegypti primarily inhabits domestic water-holding containers that are maintained in low-precipitation environments by water storage activities. The greater importance of enhanced vegetation index (EVI) for Ae. albopictus than for Ae. aegypti (15.3% vs 12.1%, Table 2) also supports the hypothesis that Ae. albopictus tends to prefer non-domestic juvenile sites. This does not, however, rule out the possibility that the two species can overlap. Additional finer scale studies need to be conducted to investigate if competitive exclusion for hosts and/or habitat occurs between Ae. aegypti and Ae. albopictus. The effect of urbanicity was surprisingly low for both species (2% and 1.1% for Ae. albopictus and Ae. aegypti, respectively). As both species have been shown to inhabit a wide variety of urban and peri-urban settings with various degrees of intensity, it is likely that the simple urban/rural distinction of our urbanicity covariate did not sufficiently capture this variation and instead continuous covariates such as EVI allow to better distinguish the respective habitat types and were thus chosen more frequently by the model. Incorporating a larger set of covariates allowed us to investigate not only the effect of temperature on survival but for additional variance as shown in the relative influence plots (Figure 1—figure supplement 1). Future Aedes species distribution models could be improved by including a comprehensive global covariate that distinguishes human settlements using complex satellite imagery processing tools.

12Our maps are based on covariates where each 5 km × 5 km pixel represents yearly mean average values. We therefore produce maps that represent the long-term average distribution of both species. However, this does not allow us to directly infer seasonal patterns of distributions which might be of importance on the periphery of the species distributions. With a more temporally resolved dataset it may be possible to capture the effects of intra-annual seasonality on the species' distributions. Adding mechanistic determinants, such as survival, have previously been used to combine seasonal patterns with global distribution maps. To make best use of the comprehensive set of data collected, we construct models and maps at a global scale, allowing the model to share information across the whole spectrum of environmental regions. However, given the scale at which this study was performed, there is always the possibility that variation in microclimate or local adaptive strategies of both species may have a significant impact in some locations.

13Previous studies have discussed the risk of pathogen importation and autochthonous transmission of DENV and CHIKV in Europe and the Americas without comprehensively accounting for the distribution of the vectors. These freely available vector distributions maps (http://goo.gl/Zl2P7J) can now be used as covariates to refine these studies and to generate high-resolution maps of the risk of possible local DENV and CHIKV transmission in currently non-endemic settings. Such maps would be useful for prioritizing surveillance in areas where there is a risk of disease importation. This will be especially important in areas where sporadic cases of related viruses have been reported, such as Europe, the United States, Argentina, and China.

14Both Ae. aegypti and Ae. albopictus have a history of global expansion associated with trade and travel. Introductions of the species over long distances and between continents has been associated with international trade routes via shipping and overland spread driven by human movement and transport routes, both facilitated by the endophilic behavior of the two species. The global spread of the associated pathogens has undoubtedly been a consequence of increasing global connectedness. As these processes continue and the world becomes increasingly connected and urbanized, risk of importation and subsequent autochthonous transmission of DENV and CHIKV will continue to increase. The true distribution of both species is influenced by a variety of factors, not just the ones presented here. Nevertheless, this study represents an important baseline for further refinements. For instance, our maps can be used to indicate areas where the species are likely to become established if introduced. Accurately predicting the future distributions of these species will also require model-based estimates of the rate at which these species colonize new areas. Such predictions can be informed by human and trade mobility patterns between endemic and non-endemic regions as well as data on the past spread of the vectors. Improving our ability to predict rates of vector importation will therefore be crucial to inferring future risk.

15Previous studies have provided crucial information on genetic variation both within and between populations of these two vector species. As the volume of georeferenced information on the population genetics of Ae. aegypti and Ae. albopictus increases, the potential to incorporate this information into mapping analyses to understand the current and future distribution of disease risk also increases. Phylogeographic analyses offer a unique way to infer the recent patterns of vector spread and to identify the major routes of importation. This information is crucial to inform models that predict the risk of vector introductions.

16Phylogenetic information could also be used to inform future iterations of the species distribution models used here by enabling the model to characterize and map environmental suitability for different vector subspecies. This could be particularly useful in the case of Ae. albopictus where genetic variation is known to underlie the ability to undergo diapause and therefore to overwinter in colder locations. Mapping the distributions of distinct genetic subgroups could also improve our understanding of the complex interactions between mosquito vector populations and virus strains and how this relates to spatial variation in transmission intensity.

17The maps presented comprise a contemporary estimate of the current and potential future distribution of Ae. aegypti and Ae. albopictus. As more occurrence data become available, these maps can be refined to incorporate recent importation and establishment events and corresponding improvements in predictions. By disseminating both the occurrence data and the predictive maps on an open-access basis we hope to facilitate both the future development of these maps and their uptake by the global public health community.

18A BRT modelling approach was applied to derive probabilistic global environmental risk maps for Ae. aegypti and Ae. albopictus. BRT models are machine-learning model ensembles commonly used in species distribution modelling (SDM) and show strong predictive performance due to their ability to handle complex non-linear relationships between probability of species occurrence and multiple environmental correlates. Our model required the following sets of input data in order to make accurate predictions of the distribution of these two species: (i) a temperature suitability mask defining the fundamental limits of both species; (ii) a globally comprehensive dataset of geo-positioned occurrence points for both species; (iii) appropriate land-cover and environmental covariate datasets that help explain the current distribution of the species; and (iv) a set of species absence records that further refine the species range and reduce sampling bias. Details regarding the specific attributes of the model and data generation are outlined below and maps of each of the covariates are shown in Figure 1—figure supplement 2.

19While the niche of a species is determined by a host of environmental, ecological and socio-economic factors of unknown influence and interaction strength, it is possible to exclude parts of the niche if the direct effects of one factor on a step rate-limiting to population persistence are well known. One such example for mosquito population persistence is whether temperature permits adult females to survive long enough to complete their first gonotrophic cycle and thus oviposit. Both adult female longevity and length of first gonotrophic cycle are temperature dependent. Combining these two relationships with a dynamic population-level simulation,, evaluated the thermal limits to persistence of Ae. aegypti and Ae. albopictus populations on a global scale. The binary outputs of this model are used as a mask to sample pseudo-absence points in locations known to be unsuitable–thereby informing the statistical model using mechanistic model outputs. The temperature suitability index developed by Brady et al. is also used in a continuous variable form (i.e., the relative number of ovipositions of parous females permitted by temperature) as a covariate in the BRT model.

20The database used for this study contains information on the known global occurrences of the adults, pupae, larvae or eggs of Ae. aegypti and Ae. albopictus globally from 1960–2014. We included data from a variety of sources, including (1) published literature and (2) primary and unpublished occurrence data from national and international entomological surveys. To our knowledge this is the largest, most comprehensive global dataset for both Ae. aegypti and Ae. albopictus. Confirmed Aedes occurrences were entered in the database after a comprehensive literature search using methods described elsewhere (;; http://dx.doi.org/10.5061/dryad.47v3c). In short, this included extracting all available location (latitude and longitude) information from the relevant articles, primarily using Google Maps (http://www.google.com/maps) so that it matched the spatial resolution of our covariate datasets of approximately 5 km × 5 km. Primary and unpublished data sources were obtained from Brazil, Europe, Indonesia, Taiwan, and the United States. After consolidating all data into two large databases for each species, independently they underwent spatial and temporal standardization. An occurrence record was defined as a single occurrence at a given unique location within one calendar year. This was important to avoid over-representation in regions where multiple surveys per year were performed, such as Taiwan or Brazil. To ensure the accuracy of the data we overlaid the geolocated occurrence points with a raster that distinguished land from water. Any records that were positioned outside the land area were subsequently removed. In total we assembled 19,930 and 22,137 occurrence records for Ae. aegypti and Ae. albopictus respectively. The distribution of occurrence points are plotted in Figure 1—figure supplement 4.

21The distribution of both species considered in this study are known to be influenced by environmental factors such as temperature and demographic factors such as urbanisation. Global gridded maps of such variables are becoming ever more available and have been commonly applied in SDM and disease mapping. The rationale for the inclusion of each variable we used is described below.

22Survival of Ae. aegypti and Ae. albopictus is highly dependent on temperature and water availability. EVI measures vegetation canopy greenness and can be used as a proxy for soil surface-level moisture that are associated with the availability of mosquito larval development sites. Eggs and adults require moisture to survive, with low dry season moisture levels affecting adult mortality. Vegetation canopy cover reduces evaporation and wind speed in the sub-canopy, which protects mosquito development sites. We used range and mean values of MODIS EVI after processing through a gap-filling algorithm described elsewhere.

23The principal larval habitats of both species are man-made containers that are used for water storage or accumulate rain. Some local studies have shown that there is a relationship between precipitation and vector abundance. To account for the availability of water-filled containers a maximum and minimum annual precipitation layer was extracted from the WorldClim database and projected for the year 2015 (http://www.worldclim.org).

24Ae. aegypti adults are highly domesticated mosquitoes feeding almost exclusively on humans, larvae develop preferentially in artificial containers in close association with human habitation, often in urban settings. Ae. albopictus are more commonly found in rural and peri-urban settings, feeding readily on a variety of mammalian and avian species, although Ae. albopictus shows similar larval development behavior in artificial containers. To account for differences in urban, peri-urban and rural environments we built a categorical variable by supplementing the projected 2010 Global Rural Urban Mapping Project (GRUMP) urban and rural categories with land-cover classes using night-time light satellite imagery and population density, using the most up-to-date national censuses available to the smallest available administrative unit available. A gridded surface of 5 km × 5 km cells was generated with each pixel representing either urban, peri-urban, or rural areas.

25BRT models consistently outperform other species distribution models such as maximum entropy (Maxent), GARP, and BIOCLIM in their predictive performance. BRT combines the strengths of regression trees (i.e., the omission of irrelevant variables and the ability to model complex interactions) with machine learning techniques (i.e., the building of an ensemble of models that approximate the true response surface). To prevent overfitting, the model used a penalized forward stepwise search and cross-validation method to identify the optimal number of decision trees. Modelling was performed using the gbm, dismo, raster and seegSDM R packages using the R v 3.1.1 environment.

26Pseudo-absence (also referred to as background) records provide a sample of the set of conditions available to the species in the region rather than actual absences. These records are needed because true absences are generally unavailable in large composite datasets such as the one used in this study. To account for reporting bias in presence data, a common problem with presence-only SDM, which if not accounted for can lead to biases in the resulting predictions, we follow in sampling pseudo-absence points according to the same reporting bias likely to be present in occurrence records (namely spatial variation in reporting of mosquito occurrence). Firstly, we selected 10,000 occurrence records of Aedes species from the Global Biodiversity Information Facility (http://www.gbif.org), omitting all records of Ae. aegypti and Ae. albopictus. This dataset is intended to reflect biases in mosquito reporting in areas which are suitable for Aedes mosquitoes. Secondly, to reflect areas where habitats are biologically not suitable for Aedes occurrence we sampled an additional 10,000 pseudo-absence points at random locations, with sampling probability greatest in areas that the biologically-based temperature suitability index predicted to be unsuitable. Thirdly, sampling of occurrence points was also biased towards oversampled regions such as Brazil and Taiwan in which there were a large number of reported occurrence records due to the inclusion of results of large national entomological surveys (Table 1). Therefore, we weighted occurrence records from these locations so that the density of occurrence records per country matched the density of all other records globally by dividing the number of occurrence points by the size of the respective countries.

27An ensemble BRT was constructed using 120 sub-models to derive uncertainty distributions of the prediction map. Each of the 120 sub-models was fitted to a separate bootstrap resampling of the dataset and used to generate a probability map for each individual species on a 5 km × 5 km resolution. The mean of these 120 sub-models was used as the final Aedes risk maps. Pixel based uncertainty was estimated by calculating the 95% confidence interval from the 120 sub-models.

28The variables used as land-cover and environmental correlates used in this study are quantified based on their relative influence (0–100) on explaining the variance in the models calculated as the sum of the number of times a particular variable is selected for splitting the decision tree, weighted by the squared improvement to the overall model averaged over all trees. Note that in a BRT, non-informative predictors are largely ignored. Predictive performance of each sub-model was evaluated using the area under curve (AUC) statistic calculated as the mean AUC for each of the ten cross-validation folds evaluated against the other 90% of the data under the pairwise distance sampling procedure of. The overall predictive accuracy of the model was measured as the mean and standard deviation of these AUCs across all 120 sub-models.