Population Dynamics of Anopheles albimanus (Diptera: Culicidae) at Ipetí-Guna, a Village in a Region Targeted for Malaria Elimination in Panamá

Anopheles albimanus Wiedemann is a major malaria vector in Mesoamerica and the Caribbean whose population dynamics, in response to changing environments, has been relatively poorly studied. Here, we present monthly adult and larvae data collected from May 2016 to December 2017 in Ipetí-Guna, a village within an area targeted for malaria elimination in the República de Panamá. During the study period we collected a total of 1678 Anopheles spp. mosquitoes (1602 adults and 76 larvae). Over 95% of the collected Anopheles spp. mosquitoes were An. albimanus. Using time series analysis techniques, we found that population dynamics of larvae and adults were not significantly correlated with each other at any time lag, though correlations were highest at one month lag between larvae and adults and four months lag between adults and larvae. Larvae population dynamics had cycles of three months and were sensitive to changes in temperature with 5 months lag, while adult abundance was correlated with itself (1 month lag) and with the normalized difference vegetation index (NDVI) with three months lag. A key observation from our study is the absence of both larvae and adults of An. albimanus between January and April from environments associated with Guna population’s daily activities, which suggests this time window could be the best time to implement elimination campaigns aimed at clearing Plasmodium spp. parasites from Guna populations using, for example, mass drug administration.


Introduction
Panamá is the nation with the highest economic growth in Mesoamerica [1], with a long history of accomplishments [2,3] and landmarks in the study [4][5][6] and control of mosquito vectors of diseases [7] that occurred in Panamá's territory during the development of the Canal Zone as a colonial possession, Figure 1. República de Panamá map, highlighting East Panamá province and the location of Ipetí-Guna near the southeastern shore of Lake Bayano. This map was made using as a base a public domain map from the US National Park Service [41]

Mosquito Sampling
Mosquitoes, adults and larvae, were monthly collected from May 2016 to December 2017. Anopheles spp. adults mosquitoes were collected by human landing catch [22], a highly sensitive mosquito sampling method for anopheline mosquitoes [42], at the peridomicile of a housing cluster in Ipetí-Guna ( Figure 2). This housing cluster was chosen given the abundance of nearby positive larval habitats. It is worth highlighting that distance from the sampling location to the huts in the housing cluster ranged between 5 m and 10 m. The adult mosquito collections were made over three consecutive nights between 18:00 h and 21:00 h by two trained collectors that switched positions to minimize bias in the abundance estimates [41]. To minimize the infection risk, catches followed World Health Organization (WHO) biosecurity guidelines [43]. Anopheles spp. larvae were collected by dipping larval habitats near the households, but also on the Ipetí-Guna River banks and in the farmlands of the community (Figure 2), following the procedure recommended by WHO, of ten dips per square m at a sampling location [43]. Larvae and adults mosquitoes were placed in Petri dishes, which were coded by collection site and date, and then mosquitoes were morphologically identified using a dissection scope and taxonomic keys for the anophelines of Mesoamerica [44,45] and the reference collection at the Instituto Conmemorativo Gorgas de Estudios de la Salud (ICGES).

Mosquito Sampling
Mosquitoes, adults and larvae, were monthly collected from May 2016 to December 2017. Anopheles spp. adults mosquitoes were collected by human landing catch [22], a highly sensitive mosquito sampling method for anopheline mosquitoes [42], at the peridomicile of a housing cluster in Ipetí-Guna ( Figure 2). This housing cluster was chosen given the abundance of nearby positive larval habitats. It is worth highlighting that distance from the sampling location to the huts in the housing cluster ranged between 5 m and 10 m. The adult mosquito collections were made over three consecutive nights between 18:00 h and 21:00 h by two trained collectors that switched positions to minimize bias in the abundance estimates [41]. To minimize the infection risk, catches followed World Health Organization (WHO) biosecurity guidelines [43]. Anopheles spp. larvae were collected by dipping larval habitats near the households, but also on the Ipetí-Guna River banks and in the farmlands of the community (Figure 2), following the procedure recommended by WHO, of ten dips per square m at a sampling location [43]. Larvae and adults mosquitoes were placed in Petri dishes, which were coded by collection site and date, and then mosquitoes were morphologically identified using a dissection scope and taxonomic keys for the anophelines of Mesoamerica [44,45] and the reference collection at the Instituto Conmemorativo Gorgas de Estudios de la Salud (ICGES).

Weather Data
Temperature and rainfall data were obtained from gridded databases. The data were extracted from the Royal Netherlands Meteorological Institute website, using the KNMI explorer [47]. The specific gridded databases that we employed were: (i) NOAA CPC Morphing Technique ("CMORPH") database to obtain monthly rainfall data, which has a resolution of 0.25° [48]; (ii) NOAA Global Historical Climatology Network version 2 and the Climate Anomaly Monitoring System (GHCN_CAMS 2m model) for temperature, which has a spatial resolution of 0.5° [49]. To create a seasonal climate profile, with monthly boxplots for the 12 months of the year, monthly data were downloaded from 1998 to 2017 for the study area.

Vegetation Data
We employed images for the NDVI from the biweekly, 250 m resolution, vegetation product (M*D13A3), courtesy of the NASA Land Processes Distributed Active Archive Center, USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota [50]. We used the R package MODIStsp [51] to download and clip NDVI images from the study period. Further GIS procedures for the downloaded images were made using the package raster also in the statistical

Weather Data
Temperature and rainfall data were obtained from gridded databases. The data were extracted from the Royal Netherlands Meteorological Institute website, using the KNMI explorer [47]. The specific gridded databases that we employed were: (i) NOAA CPC Morphing Technique ("CMORPH") database to obtain monthly rainfall data, which has a resolution of 0.25 • [48]; (ii) NOAA Global Historical Climatology Network version 2 and the Climate Anomaly Monitoring System (GHCN_CAMS 2m model) for temperature, which has a spatial resolution of 0.5 • [49]. To create a seasonal climate profile, with monthly boxplots for the 12 months of the year, monthly data were downloaded from 1998 to 2017 for the study area.

Vegetation Data
We employed images for the NDVI from the biweekly, 250 m resolution, vegetation product (M*D13A3), courtesy of the NASA Land Processes Distributed Active Archive Center, USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota [50]. We used the R package MODIStsp [51] to download and clip NDVI images from the study period. Further GIS Insects 2018, 9, 164 5 of 17 procedures for the downloaded images were made using the package raster also in the statistical software R [52]. Briefly, all images from the study period were stacked into a geotiff, from which the monthly values were extracted, generating a time series in the process [20].

Time Series Analysis
To study the abundance patterns of Anopheles albimanus adults and larvae we followed a standard procedure for time series analysis, where we first inspected the autocorrelation patterns of both time series. We started by inspecting the autocorrelation function (ACF) and the partial autocorrelation function (PACF). Briefly, the ACF depicts the correlation of the time series with itself through time, while the PACF does so considering the correlation between consecutive time lags [53]. With information from both of these functions it was possible to identify time lags associated with population size, of both larvae and adults, at any given time. This information was then used to fit a null model that was used to pre-whiten time series from all the environmental covariates described before. Pre-whitening is a procedure that removes any potential similar autocorrelation from the environmental covariates time series, thus ensuring cross-correlation functions (CCF) show associations that do not emerge from time series having similar autocorrelation structure. To estimate CCFs, we used residuals from the null models fitted to the adults and larvae data and the pre-whitened time series from each environmental covariate. Based on the resulting CCFs we chose lags to fit time series models including the covariates that were associated with mosquito abundance by stage (adults and larvae). The time series models we fitted can be described by the following equation: where y t is the mosquito abundance at time t, x indicates the time lag, µ is the average mosquito abundance during the study period, ϕ x the autoregressive coefficients, cov indicates a covariate (weather and vegetation variables described above), α x the coefficients for the covariates and ε t is the error, which was assumed to be normal, identical and independent [53]. To select the environmental covariates, we used backward elimination guided by the minimization of the Akaike Information Criterion (AIC). AIC is a parameter that trades-off model goodness-of-fit and parameter number, and whose minimization can be used to select among models with a similar number of parameters [53]. For the best model, error assumptions were verified using standard procedures for time series analysis [53].

Results
During the study period we collected a total of 1678 mosquitoes, 1602 were adults and 76 larvae. Over 95% of the collected Anopheles spp. mosquitoes were An. albimanus, 1585 adults and 71 larvae ( Table 1). The monthly abundance (±S.D.) of An. albimanus was 79.25 ± 96.40 for adults and 3.55 ± 4.70 for larvae. Other anopheline mosquito species collected at the study site included: Anopheles punctimacula Dyar & Knab, and Anopheles punctipennis (Say), both of these species collected as adults and larvae, and Anopheles apicimacula Dyar and Knab, collected only as adult (Table 1). Table 1. Monthly mean abundance of Anopheles spp mosquitoes, by stage, collected between May 2016 and December 2017 in Ipetí-Guna, Comarca Madungandí, República de Panamá. The column "Total" indicates the total number of individuals collected for each species during the study period, and the column "Persistence" the proportion of the sampling sessions when a given mosquito species was found among the samples.

Adults
Larvae The mosquito An. albimanus was not only the most abundant anopheline mosquito species, but also the most persistent, being present at the study site 75% and 60% of the times mosquitoes were sampled, respectively, as adults and larvae. An. albimanus was absent in adult mosquito samples from April to May of 2016 and from February to April of 2017 ( Figure 3A), while larvae were absent in May and August 2016, and from January to May, and July of 2017 ( Figure 3B). Meanwhile, all other Anopheles spp., had very low persistence (Table 1) and were caught mainly during 2016 as both adults ( Figure 3C) and larvae ( Figure 3D). The mosquito An. albimanus was not only the most abundant anopheline mosquito species, but also the most persistent, being present at the study site 75% and 60% of the times mosquitoes were sampled, respectively, as adults and larvae. An. albimanus was absent in adult mosquito samples from April to May of 2016 and from February to April of 2017 ( Figure 3A), while larvae were absent in May and August 2016, and from January to May, and July of 2017 ( Figure 3B). Meanwhile, all other Anopheles spp., had very low persistence (Table 1) and were caught mainly during 2016 as both adults ( Figure 3C) and larvae ( Figure 3D). The persistence and abundance patterns of An. albimanus seemed to track changes in rainfall ( Figure 4A) and temperature ( Figure 4B) at the study site, with both adult and larvae being absent during the dry season that spans from January to April ( Figure 4C) and during the hottest season, with the highest temperatures, from February to May ( Figure 4D). It is worth highlighting that both rainfall and temperature had record high values during the study period, where rainfall in July 2017, at the peak of the rainy season ( Figure 4C), was well above the median, and temperatures were hotter  The persistence and abundance patterns of An. albimanus seemed to track changes in rainfall ( Figure 4A) and temperature ( Figure 4B) at the study site, with both adult and larvae being absent during the dry season that spans from January to April ( Figure 4C) and during the hottest season, with the highest temperatures, from February to May ( Figure 4D). It is worth highlighting that both rainfall and temperature had record high values during the study period, where rainfall in July 2017, at the peak of the rainy season ( Figure 4C), was well above the median, and temperatures were hotter than usual in January, April and December ( Figure 4D) when compared with the distribution of records from 1998 to 2017.
Insects 2018, 9, x FOR PEER REVIEW 7 of 16 first order autoregressive process for the adult abundance, i.e., the time series was significantly correlated (p < 0.05) up to one month lag ( Figure 6C), while, by contrast, larvae abundance was significantly associated (p < 0.05) in a pattern that suggested cycles with a three month lag ( Figure  6D). Regarding the environmental covariates, adults were significantly correlated (p < 0.05) with NDVI with a three month lag ( Figure 6E), while larvae were associated (p < 0.05) with temperature with a five month lag ( Figure 6F).     (Figure 2), had the highest NDVI during the entire study period, followed by the river site and the farm site, which reached low values during the dry season. Autocorrelation functions from the An. albimanus abundance time series suggest that adults and larvae had different patterns of temporal abundance association, where adults had a decaying autocorrelation function ( Figure 6A), but larvae had a cyclic autocorrelation function ( Figure 6B). Both of these were patterns confirmed by the partial autocorrelation functions, which suggested a first order autoregressive process for the adult abundance, i.e., the time series was significantly correlated (p < 0.05) up to one month lag ( Figure 6C), while, by contrast, larvae abundance was significantly associated (p < 0.05) in a pattern that suggested cycles with a three month lag ( Figure 6D). Regarding the environmental covariates, adults were significantly correlated (p < 0.05) with NDVI Insects 2018, 9,164 8 of 17 with a three month lag ( Figure 6E), while larvae were associated (p < 0.05) with temperature with a five month lag ( Figure 6F).   Meanwhile, the association between adults and larvae was not statistically significant at any time lag, yet it had a peak at four month lag ( Figure 7A). The association between larvae and adults was highest at one month lag ( Figure 7B), but not statistically significant, yet higher than the four month lag association between adults and larvae. Meanwhile, the association between adults and larvae was not statistically significant at any time lag, yet it had a peak at four month lag ( Figure 7A). The association between larvae and adults was highest at one month lag ( Figure 7B), but not statistically significant, yet higher than the four month lag association between adults and larvae.
(E) adults and rainfall, temperature and NDVI, based on the data for House (F) larvae and rainfall, temperature and NDVI, based on the data for house pools, farm pools and the river. The inset legend of panel (E) shows the color coding for the different environmental variables, while the inset legend of panel (F) shows the line coding for the NDVI measured at the different sampling locations. In each plot, blue dashed lines indicate the confidence intervals within which correlations are expected to arise by random (p > 0.05).
Meanwhile, the association between adults and larvae was not statistically significant at any time lag, yet it had a peak at four month lag ( Figure 7A). The association between larvae and adults was highest at one month lag ( Figure 7B), but not statistically significant, yet higher than the four month lag association between adults and larvae. The significant lags found in the different cross correlation functions were used to fit time series models described in the following lines. Table 2 shows the results of the model selection process for the time series model explaining An. albimanus adult abundance. As shown in Table 2, the The significant lags found in the different cross correlation functions were used to fit time series models described in the following lines. Table 2 shows the results of the model selection process for the time series model explaining An. albimanus adult abundance. As shown in Table 2, the autoregressive model including NDVI with a three month lag outperformed an autoregressive model.

Variables [Lag in Months] AIC
AR [1] 234.51 AR [1], NDVI [3] 227. 60 Parameter estimates for the best An. albimanus adult abundance model are presented in Table 3, showing that An. albimanus adult abundance was positively associated with NDVI (p < 0.05), indicating that mosquito abundance increases following vegetation growth at the study site. Table 3. Parameter estimates for the best time series model explaining Anopheles albimanus adult abundance at Ipetí-Guna, República de Panamá. In parameters AR stands for Autoregressive.

)
3.85 ± 1.15 * Error Variance (σ 2 ) 3320 * Statistically significant (p < 0.05). Table 4 shows the results of model selection for variables associated with An. albimanus larvae abundance. As shown in Table 4, the best model included temperature with five month lag as a covariate, adult abundance with one month lag and a seasonal autoregressive component with three month period. Although the associations between An. albimanus stages, i.e., adult and larvae, were not statistically significant, we considered the one month lag adult abundance in the An. albimanus larvae model because the correlation was slightly larger than that observed between adults and larvae. Table 4. Model selection for variables explaining Anopheles albimanus larvae abundance at Ipetí-Guna, República de Panamá. AIC indicates the Akaike Information Criterion for the model, and Variables the covariates considered in each model. In the variables, SAR indicates Seasonal Autoregressive. The best model, with the lowest AIC, is in bold font.

Discussion
Although malaria has been significantly decreasing in Mesoamerica [54], this vector-borne disease remains a public health problem in the region, particularly in Panamá, the only country in Mesoamerica that did not reach the millennium development goal of a 75% reduction of malaria transmission [55]. In Panamá, it has become increasingly clear that malaria mainly affects the Gunas [9,22]. The vulnerability of this group to malaria is such that transmission of highly pathogenic malaria parasites, often resistant to antimalarials [13,14], has been clonal, indicative of epidemic transmission as revealed by low parasitic genetic diversity, and geographically widespread [56]. Malaria transmission patterns in Panamá, although without a clear seasonality, have interannual transmission cycles linked to ENSO [18,20], a global climatic phenomenon associated with changes in malaria transmission worldwide [57][58][59][60][61]. This situation calls for increased efforts to define and implement focalized activities to control and eliminate malaria from the Gunas, where their overarching situation of social exclusion, poverty, and different cultural beliefs [9,22] on top of their poor housing quality [12,22], exclude a fast-paced change favoring improved living conditions. For example, housing improvement [62], a change that has been decisive for malaria disappearance from diverse settings around the world [8,[63][64][65][66][67][68][69][70][71] is out of reach for Guna communities in Panamá. In this context, insights from studying malaria mosquito vector ecology could be useful to guide malaria transmission reduction [72,73].
Our longitudinal study revealed a strong seasonality in An. albimanus; its abundance greatly decreases during the dry season, January to April, a pattern observed for anophelines elsewhere in the New World [74][75][76][77][78]. This reduced mosquito abundance pattern makes the dry season ideal for interventions aimed at eliminating parasites from the human population, for example, via mass drug administration (MDA) [79][80][81]. MDA can be adapted to reduce malaria transmission among indigenous populations with unique cultural characteristics [82], like the Gunas. Indeed, models of adaptive management in changing environments have already proved successful for natural resource extraction and food sovereignty among the Gunas [83]. Timing parasite clearance interventions from Guna populations during the dry season is expected to be optimal, mainly because the likelihood of any transmission, as parasites are cleared from the human population is greatly reduced when mosquitoes are scarce [76]. Moreover, MDA has been successfully applied for malaria elimination from indigenous populations elsewhere [80,84].
Our results confirm the importance that vegetation growth changes have on regulating the abundance of host seeking adult mosquitoes [42,74], given the time series model for adults indicated that peaks in mosquito abundance follow NDVI changes, with a lag of three months. The three month lag could reflect at least a couple of phenomena, delayed density dependence and resonant generational cycles. In delayed density dependence, population growth in the larvae, which are sensitive to vegetation growth changes [36][37][38], is seen with a delay in the adults [85]. Meanwhile, in resonant generational cycles [86][87][88], environmental impacts on population size have a build-up time to be observed in the adult mosquito population and can be enhanced by a positive autocorrelation structure in the environment, something that could also explain the 5 month lag for the impact of temperature in aquatic mosquito, i.e., larvae, population growth, as has been observed in other mosquito vectors of pathogens [89]. These relationships with environmental factors are also important to design interventions aimed at reducing the contact between the Gunas and malaria vector populations, since, for example, in Sub-Saharan Africa it is known that adherence to use insecticide treated nets is influenced by perceived mosquito nuisance [90] and the concomitant perceived malaria transmission risk [91]. In that sense, the best timing for interventions aimed at reducing contact between humans and vectors, or suppressing vector populations, will be during the wet season, for example, three months after peaks in NDVI, which are correlated with peaks in host-seeking mosquitoes, as shown in this study. Among the tools used to reduce the contact between Gunas and mosquitoes, it is unlikely that insecticide treated nets, a tool with an extraordinary success in reducing malaria transmission in Sub-Saharan Africa [92,93], will be effective among the Gunas, as they sleep in hammocks. However, the use of insecticide-treated hammocks [94,95] might be an option that can be adopted by the Gunas as it fits their cultural practices [82].
We think the entomological risk of malaria transmission posed by the three other collected species needs to be further evaluated because all the Anopheles spp. we collected have been found infected with malaria parasites [6,[96][97][98], and might surge in case strategies aimed at suppressing An. albimanus vector populations are implemented in Guna communities. This phenomenon has been commonly observed when knocking down the density of dominant vector species elsewhere [27,28,99].
Finally, some limitations of our study include the limited sampling which was constrained by the lack of economic resources to sample mosquitoes across the sparse Guna villages in Lake Bayano basin, that are extremely difficult to access during the rainy season. Our study was also limited by our ability to find larval habitats at the study site which, although productive, seemed to be very scarce compared with descriptions for African [100,101] and South American [75,78] settings. Although it has been suggested that human landing catch is the best method to sample adult mosquitoes in the Neotropics [42], we think our study could have benefited by using other sampling methods for the adult populations, such as Mosquito Magnet ® traps [102]. Also, our environmental data could have been more descriptive of the local environmental conditions of our sampling sites if our budget had allowed us to use data loggers measuring local weather variables [103] instead of resorting to spatially coarsely grained records from gridded databases. We also think that a more frequent sampling, for example, biweekly instead of monthly, could have been helpful to straightforwardly link patterns of larvae and adult abundance, which at the monthly time scale, were not significantly correlated with each other.

Conclusions
Our study describes basic features of An. albimanus population dynamics, a dominant malaria vector in Mesoamerica [6,21], highlighting An. albimanus abundance association with environmental variables at a site where this disease has a large burden on indigenous Guna populations [18]. Based on our observations, we can suggest the use of measures aimed at clearing parasites from the Gunas, for example, mass drug administration, during the dry season from January to April, when host-seeking malaria mosquito abundance is low and the likelihood of transmission is reduced. In contrast, the deployment of interventions aimed at reducing human and vector contact might be more successful during the wet season, when mosquito nuisance might promote the adherence to use protective measures against mosquito bites or lead to perceptible changes in nuisance, when using tools aimed at reducing vector populations, producing a positive feedback on their use and acceptance [104,105]. Fine tuning for the deployment time of interventions affecting malaria vectors can benefit from remotely sensed data of vegetation growth, for example, using NDVI, and also with the use of temperature records, as inferred from the lagged association between these variables and An. albimanus abundance, respectively, as adults and larvae. Given our previous findings about the association of ENSO and malaria transmission at Comarca Madungandi [18] and Comarca Guna Yala [20], where ENSO, respectively, has a positive and negative impact on malaria transmission, we think studies on malaria vector population ecology are a priority in Guna Yala, since local weather conditions might have different impacts on Anopheles spp. abundance that explain the contrasting impacts of ENSO on malaria transmission. We think that areas like Cartí in Comarca Guna Yala, a town with malaria transmission [20,22] and expected to receive an influx of Gunas from Caribbean islands that are disappearing as result of climate change [106], and for which no malaria related environmental health assessment has been made, are the priority for developing studies designed to understand malaria vector ecology, in order to propose interventions that ensure malaria transmission reduction and elimination from the Gunas in the República de Panamá.