Aedes (Stegomyia) aegypti (L.) (Diptera, Culicidae) is the main vector of dengue, chikungunya and zika viruses (Liang et al., 2015). According to the World Health Organization (WHO) these arboviruses (viruses transmitted by arthropods) infect more than 300 million people in the world every year, causing more than 20,000 deaths (WHO, 2016). During 2009, a dengue outbreak spread widely throughout Argentina causing five deaths and 27,752 non-lethal infections in 13 provinces. Since then, cases are continuously reported eventually leading to the strongest dengue outbreak in Argentina’s history during 2016, with 76,962 notified cases (data provided by the National Ministry of Health of Argentina). In addition, Ae. aegypti is also involved in the explosive epidemics of chikungunya (alphavirus) (Higgs et al., 2015) and zika (flavivirus) (Gatherer and Kohl, 2016), which reinforces its role as a vector of diseases with increasing importance in the Americas and the entire world. Considering this scenario, the prevention of infections transmitted by Ae. aegypti can only be achieved by reducing or eliminating bites by infected vector mosquitoes.
Control of mosquitoes can be directed against their adult or immature aquatic stages (larvae and pupae), with a number of methods available for each approach. The promotion of peridomestic clean-up to reduce mosquito breeding sites have been routine practices, mostly as one of many elements within multiple interventions including indoor and outdoor spraying, house screening, insecticide-treated materials and biological control (Bowman et al., 2016). Although it is not possible to quantify the specific contribution of such programmes to the reduction of vector populations or their impact on virus transmission, strong evidence from Cuba (Toledo et al., 2011) indicates that community working groups, including the promotion of environmental management, reduce the vector index significantly, even more than other Ae. aegypti control programmes. Independently of the control method applied, care in designing strategies is crucial, and a clear identification of vector hotspots in populated areas is a valuable tool for the prevention of disease transmission. According to Ostfeld et al. (2005), disease risk is more closely correlated with abundance of vectors than with their presence. Also, within the complex context of arthropod-borne diseases, studies must not only be focused on infected hosts, but also on the vector mosquitoes, to finally move into integrative approaches considering the interactions between them. Geographic information systems (GIS) and remote sensing can be useful in characterizing the environmental conditions that favour the development of the biological disease vector cycle (Reisen, 2010). In other words, these technologies offer the opportunity of mapping the habitat suitability of a particular area (Cleckner et al., 2011; Espinosa et al., 2016a).
Since 2008, Mundo Sano Foundation (www.mundosano.org) implemented an extensive programme for monitoring Ae. aegypti domestic breeding sites in Clorinda, a border town in the northern province of Formosa, Argentina. In a previous study, the changes in abundance of Ae. aegypti breeding sites infested with immature stages were analyzed and different temporal behaviours within each zone in the town were recognized. It was a remarkable finding considering the small size of Clorinda and the lack of systematic control activities during previous years. Focused on that issue, the main goal of the present study was to examine, in depth, the temporal evolution of the number of houses infested with immature stages of Ae. aegypti in each individual neighbourhood, and to explore if the different behavioural characteristics found could be attributed to areas that share the same environmental conditions. To that end, the information provided by remotely sensed variables was used to group all the sampled neighbourhoods into a reduced set of environmental clusters within Clorinda evaluating their utility as a proxy for the temporal evolution of dengue epidemics.
Materials and Methods
Clorinda (25°17’S, 57°43’W) is located in the easternmost limit of the province of Formosa in north-eastern Argentina four km from the Paraguayan border (Figure 1). The town lies on the right bank of the Pilcomayo River, 10 km before its confluence with the Paraguay River. This region presents a subtropical climate without a dry season. The annual median temperature is about 23°C and the mean annual precipitation 1,300 mm due to a rainy season stretching from October to May (Gurtler et al., 2009). According to the latest national census in Argentina, conducted in 2010, the town has 16,600 homes and its population is 23,290.
From October 2011 to November 2013, the presence and abundance of Ae. aegypti breeding sites at the household level were registered throughout seven entomological surveys in Clorinda. Each time, 20% of the homes were randomly selected for monitoring, including a total of 10,981 houses in the study. Verbal consent was obtained from the dwellers in order to perform the entomological sampling. The larvae and pupae collected during the survey were transported to the entomologic laboratory of Mundo Sano Foundation for taxonomic determination using a specific morphological key (Rossi and Almiron, 2004). Houses found to have at least one container with one or more Ae. aegypti larvae or pupae in it were considered positive sites. Spatially the counts were grouped by neighbourhood and temporally by period (2011-2012, 2012-2013) as well as by season (spring, summer, autumn).
The open-source Quantum GIS (QGIS), version 2.10.1 PISA (http://www.qgis.org/) software was used to build the point vector layers (.shp files) locating sampled points, recording sites of mosquito presence, collecting abundance data for each survey and grouping them by year. Sites with presence of Ae. aegypti immature stages were indicated as positive cases and those without as negative. Based on these layers, the QGIS heatmap tool was used to estimate the density of positive cases per year. To be able to count on a visual source for the evaluation of the spatial variability of densities, the difference between layers of different years was calculated as a raster image by season. Finally, the number of positive houses was calculated by neighbourhood and year, and this value was used in the spatio-temporal modelling.
Mapping of environmental variables
An image from the satellite SPOT 5 (http://www. geoimage.com.au/satellite/spot-5), obtained on April 9 2013, was used to generate the land cover classification and macro-environmental products described by Espinosa et al. (2016b). Briefly, unsupervised classification (k-means) was performed to classify the area under study obtaining seven land cover classes: bare soil, superficial water, wetlands, low vegetation (grass), high vegetation (bushes and trees), urban buildings and pasture/crops. Through the use of Google Earth, a set of 40 control points (ground truth) was created for each environmental cover to corroborate the accuracy of the classification. To this end, a confusion matrix was created, in which the classification results were compared to the additional ground truth information in order to identify the nature of any classification error and its quantity. QGIS and the image processing software ENvironment for Visualizing Images (ENVI), version 4.8 by Exelis Visual Information Solutions (https://exelisvis.com/) were used to create vectors and classification images, respectively, to assess the accuracy level of the classifications; then two different types of macro-environmental variables were generated for each class. The distance from each pixel to the class of interest was calculated using the ENVI buffer zone image tool resulting in a raster image in which every pixel had a value defined as the distance from that pixel to the nearest pixel of the selected class; then the percentage of each land cover class was calculated by applying the methodology described by Espinosa et al. (2016b). The window size for the generation of the percentage map was 31×31 pixels, attributing the average value of the window pixels to the central pixel. Shannon’s diversity index, a quantitative gauge of the different categories in a dataset, was calculated as a measure of landscape heterogeneity, using LecoS python plug-in for automated landscape analysis in QGIS (Jung, 2013). In addition, the normalized burn ratio thermal (NBRT) layer was obtained as proxy for the annual land surface temperature. This index includes near infrared, short wave infrared (SWIR) and thermal infrared Landsat bands (http://landsat.gsfc.nasa.gov/). These composites are made from Level L1T orthorectified scenes, using the computed top-of-atmosphere reflectance (Holden et al., 2005) and can be obtained from several repositories of Earth data, for example that in Google Earth Engine. A vector layer of the sampled neighbourhoods in Clorinda was used to extract the mean value of each macro-environmental variable, using the zonal statistics tool available in QGIS.
A subset of variables with low correlation (r<0.7) and/or biological importance in the Ae. aegypti ecology was chosen for subsequent analyses. All statistical analyses were performed in R (www.r-project.org) after standardizing the mean values of the chosen environmental variables in order to account for unit differences. The matrix of the values assigned to the different neighbourhoods was submitted to k-means clustering (Hartigan and Wong, 1979). The package cluster (Maechler et al., 2017) was used to graph the arrangement of neighbourhoods obtained from the analysis.
Generalized linear model
A generalized linear model (GLM) was fit to assess the contribution of each of the environmental clusters to explain any change in abundance of Ae. aegypti breeding sites. Absolute counts were corrected by the total number of sampled houses from one year to the next. As frequently occurs with counts, our data displayed overdispersion. Since the imposition of a Poisson distribution would be inappropriate, a negative binomial error distribution was chosen for handling this issue. The different clusters and the two years of sampling were set as factors and used as explanatory variables of the number of positive houses (the response variable). The model was run using the MASS library (Venables et al., 1999) and validated throughout the graphical exploration of the distribution of residuals (Zuur et al., 2010). Multiple comparisons with Bonferroni correction were performed between interaction terms using the multicomp (Hothorn et al., 2008) and factorplot (Armstrong, 2016) packages.
The variation in the observed densities of positive houses in the neighbourhoods sampled over six seasons during two consecutive years (2012 and 2013) is represented in Figure 2, which shows that different neighbourhoods had different temporal abundance patterns. For example, some near neighbourhoods, such as 7-8, 14-15-16 or 21-22, show differences between comparable seasons. The heatmaps shown in Figure 3 represent an interpolated surface showing the changes of density of positive houses from Clorinda between the two years with sampled data, which were also separated by seasons.
A total of 17 environmental layers were obtained after performing the unsupervised classification of the SPOT 5 scene containing Clorinda. The confusion matrix used with validation control points showed an overall classification accuracy of 77.6% and a Kappa coefficient of 0.74. Some classes (water, pasture/crops, low and high vegetation, wetlands) reached values above 80% accuracy, while others (bare soil and urban buildings) presented 50% and 76% accuracy, respectively. From the complete set of variables, the following eight showed low correlation (r<0.7) and were chosen for clustering the sampled neighbourhoods by k-means partitioning: distance to surface water, distance to wetlands, distance to areas with buildings, distance to tall vegetation, short-wave infrared, percentage of wetlands, percentage of tall vegetation and Shannon diversity index.
The k value that incorporates a high proportion of the variability among groups was set to 4 for the subsequent environmental class assignation of neighbourhoods, since higher values would lead to less parsimonious arrangements or overlapping among classes. Then, the matrix of values of the environmental variables per neighbourhood was classified by k-means partitioning, which resulted in the arrangement of the 32 sampled neighbourhoods into four environmental clusters (Figure 4).
Cluster k-I was composed by middleclass residential areas and a commercial zone mainly occupied by storehouses and car repair shops; the town’s cemetery is also located within this group (neighbourhood 11). Cluster k-II was represented by neighbourhoods with intense commercial activity in precarious houses, and floodable areas near the rivers, while cluster k-III was mostly composed by densely populated social housings (neighbourhoods 6, 7, 24, and 31) and a settlement of the native American Toba people (neighbourhood 26). Finally, cluster k-IV harboured peripheral neighbourhoods, such as farms and crops, all located in areas that must often be evacuated because of their proximity to the Pilcomayo River. The corresponding mean vectors of each cluster are shown as centroid values in Table 1.
The GLM results showed that the global interaction between years and environmental clusters was statistically significant (LR Chisq: 23.69, Df: 3, P<0.001), which means that the number of positive houses is jointly affected by the interaction between the environmental clusters and the year of sampling. When we performed multiple comparisons among the different interacting levels (each year with each environmental cluster) we noted that only the neighbourhoods included in cluster k-I were significantly different between the years (Figure 5). More precisely, the number of positive houses in k-I was 9.48 times higher (P<0.005, SE=0.37) in 2013 than in 2012. We did not observe any other significant increases among the pair-wise comparisons.
The proportion of variability explained by the model was satisfactory (D2= 0.75). The non-random distribution in the frequencies of infested houses among clusters was taken into consideration, thus the overdispersion (ф=1.19) was modelled with a negative binomial residual distribution. Also, the diagnostic plots of residuals from the model did not show any pattern that suggested violation of the normal error distribution and the homogeneity of error variance.
Discussion and Conclusions
It is widely accepted that the ecology of vector systems is complex and subject to rapid changes related to alterations in the environment, such as land use, urbanization, transport, and climate (Weaver, 2013; Estrada-Pena et al., 2014). It is therefore important to develop a good understanding of vector systems in their context, which involves adopting an ecological approach, which requires cooperation with specialists from diverse fields such as microbiology, entomology, climatology, ecology, urban planning, social science, political science, public health, etc. In the case of mosquito vectors, one of the challenges facing health authorities in urban areas is surveying potential breeding areas for habitats and introduction of control measures. In the present study, we assessed the utility of remotely sensed data to explore the use of environmental clustering of neighbourhoods in a small town as a proxy for increases in the number of houses with breeding sites infested with immature forms of Ae. aegypti.
Although Clorinda is a relatively small town, the differences in the number of houses that harbour containers infested with Ae. aegypti larvae (positive houses) between years provide evidence for the variation patterns structured throughout the town. In a previous study Espinosa et al. (2016b) analyzed the absolute values of positive houses by year, comparing five arbitrary geographical areas of Clorinda. For a more detailed overview of the abundance of infested breeding sites, we analyzed the differences among every sampled neighbourhood, taking into account the values obtained per season. The variation in the observed densities of positive houses not only showed that different areas are characterized by divergent behaviour (Figure 2), but also that the density increase of positive houses involves the whole town during springtime (Figure 3). These changes appeared to be more intense and localized in specific areas during summer and autumn, with maximum values tending to increase over time. According to data provided by the Argentinean National Weather Service, the maximum temperatures during summer and autumn were higher in 2013 than in 2012, with little variability in the average monthly temperatures, which is of interest as it is known that high temperatures have an effect on survival and demographic productivity (Brady et al., 2013). Although no relevant differences were observed with respect to total rainfall between the years (1,376 mm in 2012 and 1,561 in 2013), precipitations during the first year were accumulated in the summer-autumn period, while in 2013 they were concentrated during winter and spring. The fact that the rains started earlier in 2013 than in 2012 could explain the early increases in the number of positive houses during springtime; however, it does not account for the spatial structuring shown in the buffer distance maps. Taken together, these observations demonstrate that each particular region has its own spatio-temporal dynamics. The shifting hotspots of Ae. aegypti infestation significantly challenge the identification and targeting of key premises as demonstrated previously in Iquitos, Peru, where the identification of geographical areas with historically high levels of viral transmission may be more effective than targeting mosquito hotspots (LaCon et al., 2014). Based on this spatio-temporal dynamics, we oriented the subsequent analyses to investigate the possible association between the temporal evolution of the abundance of positive houses and a set of macro-environmental layers obtained from remotely sensed images. Interestingly, the clusters obtained by k-means partitioning (Figure 4) were consistent with field descriptions provided by the local headquarter of Mundo Sano Foundation, according to which these areas have their different characteristics. As stated above, one environmental cluster (k-I) showed a significant increase in the number of positive houses between the two years under study (Figure 5). The neighbourhoods belonging to this cluster were characterized by short distances to areas with buildings, far to tall vegetation and surface water and with a low percentage of this type of environment (Table 1). It is well accepted that, with the exception of the sylvatic ancestral form, Ae. aegypti mosquitoes are extremely anthropophilic, meaning that they find the favourable conditions for the development of their life cycle in human settlements, especially in domestic arrangements (Powell and Tabachnick, 2013). In this sense, it is expected that environmental variables related to areas with buildings are good proxies for Ae. aegypti abundance and distribution, as demonstrated in previous studies based on remote sensing (Fuller et al., 2010; Landau and van Leeuwen, 2012; Espinosa et al., 2016a). Vegetation affects the distribution patterns of Ae. aegypti depending on the urban features of the area under study, independently of the geographical region surrounding the study area. For example, Hayden et al. (2010) used oviposition traps to evaluate the importance of microclimate and human factors in predicting Ae. aegypti distribution in an arid environment. They found that the species presence was positively associated with highly vegetated areas within the neighbourhoods. In our study, cluster k-I was not associated with treecovered areas. Although Clorinda is located on the subtropical, humid Chaco region where the physical conditions are optimal for the development of mosquito life cycles, the typically thick vegetation in this phytogeographical region offers the ecological niche for sylvatic species, but not necessarily for domestic mosquitoes that find their breeding requirements within and around urban structures. Similarly, natural water bodies, such as the rivers and wetlands around and within Clorinda are not suitable as breeding sites for a domestic species that prefers artificial water containers for oviposition. Water storage in tanks due to irregular supply is a recurrent practice in Clorinda, but even more evident in neighbourhoods 4, 30 and 5, included in cluster k-I. In previously reported surveys, neighbourhood 4 (Primero de Mayo) showed high infestation levels relative to the rest of Clorinda (Garelli et al., 2009), and these values were associated with the common presence of ground-level water storage tanks (Garelli et al., 2011). As mentioned above, the town’s cemetery is also located within this cluster, offering high availability of containers for Ae. aegypti breeding. In addition, the presence of storehouses and car repair shops increase the environmental heterogeneity and the availability of breeding sites. These results could indicate that the synthesis achieved with macro-environmental predictors is a good proxy for events occurring at the micro-environmental level, and that a simple and robust environmental clustering could allow for easy and rapid identification of the most problematic areas within the town when operative decisions must be taken.
Although environmental variables can be useful to explain, at least partially, the temporal behaviour of mosquito populations, future studies should include the community participation in demographic and socioeconomic surveys during the monitoring of mosquito breeding sites (Obenauer et al., 2018). When studying anthropophilic vectors it is recommendable to include socioeconomic variables for modelling the spatiotemporal behaviour of local populations. According to the WHO, various factors other than environmental ones have been determined to influence a community’s vulnerability to dengue epidemics (World Health Organization, 2009). In this sense, variables like human population density, settlement characteristics, conditions of land tenure, housing styles, education, and socioeconomic status are important for planning and for assessing dengue risk. For example, in a previous study (Espinosa et al., 2016a), a deficit in the supply of potable water was associated with high abundance of mosquito breeding sites in particular areas of Clorinda. This type of information is recorded only every ten years in Argentina, by means of a national census, and the geographic scale in which it is taken is unsuitable for modelling mosquito species distribution at the micro-habitat level.
In summary, when sanitary decision-makers deal with problems related to vector-borne diseases, normally the main goal of geospatial analysts is to construct a risk map based on geographic in situ information and environmental variables obtained by different means, including Earth-observing satellites. A typically implicit assumption is that a given map, even if obtained for a specific time, could still be valid much longer, as if relative response values were temporally invariant. However, we observed that the pattern of abundance of positive houses is spatially and temporally variable. These differences could be attributed to: i) differential responses determined by the composition of macro-environmental variables, such as those analyzed in the present study; ii) the effect of processes only observable at the microhabitat scale (Phillips et al., 2006); iii) issues related to entomological surveys and the difficulty to obtain a precise picture of the entomologic situation across an entire town or city for a given time (Vezzani et al., 2005). In this sense, future studies should take into account observations of microhabitat characteristics that may affect the suitability of containers as breeding sites for Ae. aegypti, sustained longitudinal entomological surveys, including the use of oviposition traps, and the incorporation of sociodemographic explanatory variables.