Cutaneous Leishmaniasis and Sand Fly Fluctuations Are Associated with El Niño in Panamá

Background Cutaneous Leishmaniasis (CL) is a neglected tropical vector-borne disease. Sand fly vectors (SF) and Leishmania spp parasites are sensitive to changes in weather conditions, rendering disease transmission susceptible to changes in local and global scale climatic patterns. Nevertheless, it is unclear how SF abundance is impacted by El Niño Southern Oscillation (ENSO) and how these changes might relate to changes in CL transmission. Methodology and Findings We studied association patterns between monthly time series, from January 2000 to December 2010, of: CL cases, rainfall and temperature from Panamá, and an ENSO index. We employed autoregressive models and cross wavelet coherence, to quantify the seasonal and interannual impact of local climate and ENSO on CL dynamics. We employed Poisson Rate Generalized Linear Mixed Models to study SF abundance patterns across ENSO phases, seasons and eco-epidemiological settings, employing records from 640 night-trap sampling collections spanning 2000–2011. We found that ENSO, rainfall and temperature were associated with CL cycles at interannual scales, while seasonal patterns were mainly associated with rainfall and temperature. Sand fly (SF) vector abundance, on average, decreased during the hot and cold ENSO phases, when compared with the normal ENSO phase, yet variability in vector abundance was largest during the cold ENSO phase. Our results showed a three month lagged association between SF vector abundance and CL cases. Conclusion Association patterns of CL with ENSO and local climatic factors in Panamá indicate that interannual CL cycles might be driven by ENSO, while the CL seasonality was mainly associated with temperature and rainfall variability. CL cases and SF abundance were associated in a fashion suggesting that sudden extraordinary changes in vector abundance might increase the potential for CL epidemic outbreaks, given that CL epidemics occur during the cold ENSO phase, a time when SF abundance shows its highest fluctuations.


Introduction
Cutaneous leishmaniasis (CL) is a major neglected tropical disease [1] with a complex ecology [2], whose transmission, in the New World, requires the co-existence of vectors, reservoirs and humans [3,4].In Panama ´, detailed studies on the eco-epidemiology of the disease [2,5] described parasitological aspects of reservoirs [6] and vectors [7] and the environmental context of parasite-reservoir-vector interactions [8].These studies set several landmarks for understanding New World CL epidemiology, including the demonstration of two toed sloths, Choloepus hoffmanni [9,10] and other mammals [11,12] as reservoirs of Leishmania spp parasites.Insights on sand fly vector ecology included: catholic bloodfeeding patterns in dominant vector species [13,14], very limited dispersal of flying adults [15], high sensitivity of larval biology to environmental factors [16,17,18,19,20], differential adult resting behavior with species segregating along tree height [21], large diversity of co-occurring sand fly species across CL transmission foci [22,23,24] and heterogeneities in species composition across landscape gradients [23].Currently, the most common parasite causing CL in Panama ís Leishmania panamensis [25,26] and the resurgence and exacerbation of disease transmission has led to renewed efforts aimed at improving vector control [27] as a measure to reduce transmission at emerging transmission hotspots [25].However, larger questions about what is driving the resurgence of the disease, and how to best predict epidemic outbreaks remain unanswered [28,29].
The longitudinal nature of eco-epidemiological studies on CL in Panama ´ [2] revealed interannual patterns of variability in reservoir infection and abundance [10] and the sensitivity of sand fly density to weather fluctuations [30].Nevertheless, no study has addressed if large scale meteorological phenomena, such as El Nin ˜o Southern Oscillation (ENSO), are associated with interannual CL epidemic cycles, as observed in neighboring Costa Rica [28,29,31,32] nor the impacts of ENSO on Sand Fly populations.These issues are of special interest, given the increased reports of direct and indirect impacts of climatic variability patterns associated with global warming on the disease transmission, both in the New World [32,33,34,35,36,37,38] and the Old World [39,40].Here, we employ an 11 years long (2000-2010) monthly time series recording CL cases in the Repu ´blica de Panama ´to investigate seasonal and interannual cycles on this disease.During this time at the coarse spatial scale forest cover has marginally increased, yet locally some areas have seen deforestation [41], an activity usually linked with CL outbreaks [42].We specifically ask which climatic factors are associated with seasonal and interannual disease cycles, considering local temperature and rainfall records and sea surface temperature 4 (SST4), an index associated with ENSO activity in the region.We also ask whether dominant sand fly (SF) vector species undergo marked abundance changes associated with ENSO, that are ultimately reflected in CL transmission.We employed seasonal autoregressive models and cross wavelet coherence analysis to depict the association of CL with climatic factors, and found that while seasonal CL patterns were mainly associated with temperature and rainfall, interannual cycles of the disease were associated with SST4.Moreover, SST4 was also associated with interannual cycles in temperature and rainfall.SF vectors showed marked abundance changes associated with ENSO, where abundance in general decreased during the hot and cold phases of ENSO.Our results highlight both the general association of ENSO and weather patterns with CL dynamics in Central America, and how changes in SF vector abundance associated with ENSO might play a role in the emergence of CL epidemics.

Data
Monthly CL cases were compiled by the Epidemiology Department of Panama ´'s Ministry of Health for the period January 2000 -December 2010.Briefly data consisted of cases clinically diagnosed [43], and often confirmed by the microscopic examination of skin lesson scrappings/biopsies, Montenegro skin tests (MST) [44] or Indirect Immuno-Fluorescent Agglutination Tests (IFAT) [26].Data (Figure 1A) were collected from all the health facilities administered by Panama ´'s Ministry of Health and all data came from passive case detection.Reports were then compiled at the health area level (the operational geographical units of Panama ´'s Ministry of Health which are slightly different from Panama ´'s provinces and autonomous indigenous comarcas).Slightly over 80% of the cases came from West Atlantic Panama ´, the area facing the Caribbean Sea, west of the Panama ´Canal up to the border with Costa Rica [25].Representative samples from all over Panama ´ [25] indicate that over 95% (at least 90% for each reporting health area) of the CL cases in the time series are due to Leishmania panamensis.Cases due to Leishmania mexicana, Le amazonensis and Le colombiensis continue to be rare and sporadic, as observed in earlier epidemiological studies in Panama ´ [5].Moreover, all CL cases observed in migrants that likely acquired the infection in Panama ´have been typified as Le panamensis [45,46].Temperature data were obtained from the US National Oceanic and Atmospheric Administration, NOAA (ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2/) for Tocumen (Station 787920), Albrook (Station 783842 and 788060), Bocas del Toro (Station 787935), David (Station 787930) and Santiago (Station 787950).These daily time series were averaged per month, considering all values and the monthly maxima and minima (Figure 1B).Rainfall data, an average considering all weather stations, were obtained from ETESA, Panama ´'s electrical company (Figure 1C).Monthly SST 4, often referred as El Nin ˜o 4 (Figure 1D), was obtained from (http://www.cpc.ncep.noaa.gov/data/indices/ersst3b.nino.mth.81-10.ascii).The NOAA data for SST4 were collected from the area delimited by 5uNorth-5uSouth and 160uEast-150uWest of the Pacific Ocean.We also classified each month from the time series into the ENSO phases following the Oceanic Nin ˜o Index (ONI) from NOAA (http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ensoyears.shtml).Briefly, the ONI is estimated by detrending the monthly Sea Surface Temperature 3.4 (collected in the area defined by 5uN-5uS, 120u-170uW) over periods of 30 years, with residues with a value over 0.5 corresponding to the warm (or hot) ENSO phase, residues below 20.5 corresponding to the cold ENSO phase and residues in the interval [20.5,0.5]being normal conditions [47].
SF abundance data came from six studies performed either at Universidad de Panama ´or Instituto Conmemorativo Gorgas de Estudios de la Salud (ICGES).All these studies were performed within the Repu ´blica de Panama ´(Figure S1A).All the studies used a common sand fly sampling method, which consisted in the use of unbaited light traps placed at 1.5-2.0m height above ground, with traps operating from 6 pm to 6 am (a sampling effort referred as trap-night), and placed in three well defined eco-epidemiological environments: (i) domiciliary for samples from inside houses; (ii) peridomiciliary for samples collected in a radius of 100 m from a house; and (iii) forested areas for samples collected in areas with primary/secondary vegetation and outside a 100 m radius from a house.This standardized sampling, in principle, renders the comparison of the different datasets plausible.In most of the studies sand fly sampling was done with the purpose of describing the fauna at endemic locations [48,49,50], as part of the evaluation of SF control trials [27] and we also report unpublished data from entomological surveillance of endemic CL locations in 2007, 2009, 2010 and 2011.We focused our analysis in the abundance of Lutzomyia trapidoi, Lu gomezi and Lu panamensis (Figure S1B), the dominant vector species for CL in the Repu ´blica de Panama [2,7].This was done given the proven vectorial role of these species [2,7], and the lack of abundant records for other sand fly species.Figure S1C shows the eco-epidemiological environments sampled in each location and Figure S1D the year when samples

Author Summary
In this study we analyze data on sand fly (SF) abundance and cutaneous leishmaniasis (CL) cases from Panama ´.We asked whether weather patterns and climatic variability could have an impact on vector abundance that is ultimately reflected in CL transmission.We found that large epidemics of CL occur during the cold phase of El Nin ˜o Southern Oscillation (ENSO), with regular cycles coinciding with the oscillation patterns of ENSO.By contrast, and counterintuitively, we found that during the hot and cold phase of ENSO the average number of SF was reduced when compared with the normal phase of ENSO.Nevertheless, the cold ENSO phase also shows the largest variability in vector abundance, which is a likely indicative of sudden and extraordinary increases in SF abundance.We, therefore, propose that marked SF abundance changes, triggered by anomalous weather patterns associated with ENSO, likely play a major role in shaping CL interannual transmission cycles in Panama ´.
were collected at each location.A total of 12580 SF were collected over 640 trap-nights.In all the studies, the spermathecae or genitalia were inspected for species identification using the key by Young and Duncan [51].

Statistical Analysis
Leishmaniasis seasonality, interannual cycles and nonstationary associations with weather records and ENSO.Seasonality in the CL cases time series was assessed with monthly boxplots [52] and seasonal time series plot depicting the ENSO phase.The correlation structure in the CL time series was assessed by inspecting both the autocorrelation function, ACF, i.e., the time series correlation with itself through time, and the partial autocorrelation function, PACF, i.e., the correlation between consecutive time lags [53].Information from the ACF and PACF was used to choose the time lags necessary to fit a null model, i.e., without climatic covariates, of the CL time series.The null model coefficients were used to pre-whiten the climatic covariates SST4, Rain, TMAX and TMIN.Pre-whitening is a process that removes the common autoregressive structure of an ancillary time series (often referred as filtering) thus easing the study of association patterns with a focal time series [53].Null model residuals together with pre-whitened residuals from the climatic time series were then used to estimate cross-correlation functions, CCF, of CL with each one of the climatic covariates.This information was employed to build a full model that was simplified to avoid having more parameters than what is necessary to understand the dynamics of the CL time series, i.e., to avoid over-parameterization [54].We used backward elimination of climatic covariates for model simplification, i.e., taking the least significant covariate by rounds [53].We used the Akaike Information Criterion, AIC, to choose the best model in each round of backward elimination.AIC is a model selection criterion that picks the best model, in a group, once a trade-off function between the number of parameters and likelihood is minimized [53].Finally, in all cases, assumptions about model error, normality and independence, were verified using standard procedures for time series analysis [53].To better understand the process of model association we studied the correlation of the variables considered in the full model by estimating the correlation between all variables pairs and plotting the correlations in color and size categories according to their and magnitude, building a corrgram [55].Finally, time-frequency association of CL with climatic covariates was studied using continuous wavelet transforms [56] to estimate a cross-wavelet coherence.This analysis depicts the nonstationary association between time series [56,57], i.e., their non-constant association through time, specially by depicting the coherence, i.e., association of cycles between two time series over time [28].
Sand fly vector abundance and ENSO.We fitted a Poisson Rate Generalized Linear Mixed Model (PRGLMM) [58] to abundance records from each of the dominant SF vector species (Lu gomezi, Lu trapidoi and Lu panamensis).We chose PRGLMMs given the counting nature of our data [58] and to consider the lack of independence that emerges from observations collected across different studies [59].The model was a rate model given that for a couple of studies [30,49] we had records that came from several traps and, therefore, for parameter estimation we needed to include an offset variable to account for the heterogeneous number of traps generating the counts [60].In the models we considered the following fixed factors: (a) ENSO phase, to quantify patterns of interannual variability in SF abundance; (b) sampling month, to account for seasonal changes in SF abundance and (c) the sampling eco-epidemiological environment, to account for variability that can be attributed.We considered as random factors the variability that can be attributed to: (d) each study and the nested variability that can be attributed to each location within each study in order to consider the spatial variability that emerges from observations across different locations; (e) the temporal variability that could have emerged from samples collected in different years.For the analyses, we only considered SF abundance from traps used in the control treatment of the SF vector control trial study [27], i.e., records from houses no subjected to any pesticide, while we used data from all the light traps used in the other studies [47,48,49] and from unpublished vector surveillance records.

Sand fly vector abundance and cutaneous leishmaniasis
cases.To study the relationship between CL incidence and SF vector abundance, we employed the SF trap records employed to study the relationship between SF abundance and ENSO and estimated average numbers of SF/trap-night/month.We were able to obtain 35 Lu gomezi monthly abundance estimates (MAES), 36 Lu trapidoi MAES and 37 Lu panamensis MAES, in the period January 2000 -December 2010.We employed this data together with data from the CL time series to study their cross-correlation patterns [53].From our previous mathematical modelling efforts we expected a positive relationship between the number of cases and vector abundance [4,61,62,63,64], with the abundance of vectors leading the number of cases with a time lag [29].

Results
Figure 2 shows CL case seasonality.Figure 2A shows how CL cases peak at the beginning of the rainy season in Panama ´, in April and May [41].From January 2000 to December 2010 a total of 26140 CL cases were recorded in Panama ´. Figure 2B shows how ENSO phases have been uniformly distributed across the year, and how peaks in CL cases tend to occur during the cold ENSO phase, a pattern also shown in Figure S2 when pooling all months.Given the monthly nature of our data and the focus on the total number of reported CL cases, with observations spanning 132 months, we employed a battery of tools for time series analysis in the time domain, including temporal correlation functions and seasonal autoregressive models, and time-frequency domain, i.e., wavelets [65].Figure S3 show the different correlation functions employed to fit seasonal autoregressive models to the CL time series.The PACF (Figure S3A) suggested a seasonal autoregressive structure in the CL time series, where up to the first three lags and the seasonal lag (12 months) were associated with CL number at any time.The autoregressive pattern was also suggested by the ACF (Figure S4).Indeed a 3 rd order seasonal autoregressive model was selected as the best null model (Table S1) and employed for pre-whitening climatic covariates and for subsequent estimation of Cross Correlation Functions (CCFs) between CL and climatic covariates (Figure S3B, S3C, S3D, S3E, S3F).CL was negatively correlated with SST4 at lags 4 and 5 and positively correlated with temperature at lag 12 (Figure S3B).CL was negatively correlated with rainfall (Rain) at lag 15 (Figure S3C).CL was positively correlated with average (Temp) and maximum temperature at lag 13 (Figure S3D & S3E) autonomous from minimum temperature (Figure S3F).Thus, based on the significant correlations observed in Figure S3 the following null model was fitted: Where m is the mean value of the time series, w's indicate autoregressive parameters, b, c and a are parameters for climatic covariates and e indicates an error normally distributed and with variance s 2 e .After several rounds of backward elimination (Table S1) the following model was selected as best: Whose parameter estimates are presented in Table 1.The impacts of Temp(t-13) and SST4(t-12) on CL were positive, in contrast with SST4(t-5) which was negative, and assumptions about the error were met, ensuring a sound inference.The process of model selection (Table S1) left out variables that were strongly correlated (Figure 3) with the ones present in equation (2).For example, SST4(t-5), SST4(t-4) and Rain (t-15) were positively associated between them (Figure 3) and negatively with CL(t).Thus, SST4(t-5) was able to capture a common impact of ENSO on both rainfall and CL.Similarly, Temp(t-13) and CL(t-12) were positively associated between them and with CL(t) (Figure 3), rendering the inclusion of the former, in addition to SST4(t-5), enough to account for seasonality in CL(t).
The cross wavelet coherence analysis (Figure 4) confirms the outcome of autoregressive models, showing that CL was associated with: SST4 (Figure 4A), rainfall (Figure 4B) and temperature (Figure 4C) during the study period, all the three variables associated with the inter-annual variability in CL, specifically cycles with period between 2-4 years, the latter two also associated with the seasonal cycles (period of 1 year).Similarly, rainfall (Figure 4D) was associated with SST4 at inter-annual scales; and both rainfall (Figure 4D) and temperature (Figure 4E) were seasonally associated with SST4.In synthesis, ENSO, measured through SST4, shows an imprint on CL transmission robustly revealed by both the autoregressive model and the cross wavelet coherence analysis, where SST4 also impacts Rainfall and Temperature, which are, as well, associated with CL.
Data on SF vectors showed a common pattern where Lu gomezi (Figure 5A), Lu trapidoi (Figure 5B) and Lu panamensis (Figure 5C) reduced their average abundance during the hot and cold ENSO phases.Nevertheless, these three dominant SF vector species showed an increased variability in abundance during the cold ENSO phase, i.e., the boxes in the boxplots were longer, and large outliers were more frequent during the cold ENSO phase.After controlling for differences related to heterogeneity in eco-epidemiological sampling environment, and for seasonality associated with different sampling months, the GLMPMs for Lu trapidoi and Lu panamensis showed that abundance reductions were significant during both the hot and cold ENSO phases (P, 0.05), but not for Lu gomezi, (Table 2).The abundance across ecoepidemiological environments showed a similar pattern for Lu gomezi (Figure 5D) and Lu trapidoi (Figure 5E), where abundance was slightly larger in domiciliary than in peridomiciliary or forest environments, a pattern statistically significant (Table 2).By contrast, Lu panamensis (Figure 5F) was, respectively, about 9 and 2.5 times more abundant in forests and peridomiciles than inside the houses (P,0.05,Table 2).Table 2 shows that, in general, the unexplained variance in the SF abundance GLMPMs associated with spatial heterogeneity (Location Var) was about one order of magnitude higher than the unexplained temporal variance (Year Var).Similarly, model selection with AIC and BIC showed that it was not necessary to nest the spatial random effects within each study (Table S2).
Monthly CL case records were positively associated with the abundance of Lu gomezi (Figure 5SA) and Lu trapidoi (Figure 5SB) with the CCF peaking at three months lag, i.e., after a peak in SF vector abundance there was a peak in CL cases three months later.For Lu gomezi (Figure 6A) the 3 month-lagged pattern of association with CL cases was mainly linear, but for Lu trapidoi (Figure 6B) the number of cases seemed to flatten out at high SF abundances.By contrast, monthly CL cases and Lu panamensis abundance were negatively and significantly associated with a one and a two month lag (Figure 5SC).As expected from the CCF, the association between CL cases and the 3 month-lagged abundance of Lu panamensis had no clear pattern (Figure 6C).When the abundance of Lu gomezi and Lu trapidoi were added together the maximum positive association occurred at three months of lag, but the positive association was also significant at 2 and 4 months of lag (Figure 5SD).The 3 month-lagged pattern of association between CL cases and the combined abundance of Lu gomezi and Lu trapidoi (Figure 6D) resembled the one observed for Lu trapidoi alone (Figure 6B) with the number of cases flattening out at high SF abundance.

Discussion
A general criticism for studies addressing the impact of climate change on vector-borne disease transmission is that little to no attention has been given to what changes, if any, occur in entomological risk patterns [66], e.g., what happens to the vectors during ENSO?.Here, we tried to address that knowledge gap for CL by going beyond the description of the association between ENSO and weather patterns and CL epidemics in Panama ´, we inquired whether SF vectors change their abundance during ENSO.Our data showed that interannual cycles of CL transmission, as inferred from a CL case time series from Panama ´, were associated with ENSO, a pattern observed in neighboring Costa Rica [29], and also observed for malaria in the Repu ´blica de Panama ´ [67], highlighting the impacts of ENSO on vector-borne diseases in Central America.Large CL Epidemics were observed during the cold ENSO phase or shortly after it, where the delay might reflect the delay between transmission and clinical symptoms in American CL [29,43,68], a possibility further reinforced by the 3 month delayed association between vector abundance and CL incidence.
Seasonal (intrannual time scale) changes in CL transmission were associated with temperature, a weather component with low variability, i.e., low amplitude fluctuations, in the Panama ´isthmus.This pattern may make sense in light of Schmalhausen's law, the notion that biological systems are more sensitive to small changes in low variability factors when stressed by other environmental components [65].SF population dynamics may become more sensitive to changes in temperature given their need to cope with more marked changes in other weather factors, e.g., rainfall which has a more marked seasonal imprint than temperature in Panama [41].Indeed, the pattern of higher sensitivity to changes in temperature in places with marked seasonality in rainfall has been observed for other disease vectors [69].Here, we want to also note that large CL epidemics occurred during or shortly after the cold ENSO phase, a time when, on average, SF vector abundance is the smallest across ENSO phases.Nevertheless, as observed in the raw data, the cold ENSO phase is a time when SF vectors are also prone to show extremely large abundance records per trap-night, which might reflect insect population outbreaks [70,71], i.e., sudden extraordinary increases in vector abundance [72].Thus, the occurrence of large CL epidemics during or shortly after the cold ENSO phase might indicate a role for SF vector outbreaks on CL epidemics.Indeed, a detailed study in Venezuelan village showed that CL cases in an endemic village were associated with vector abundance [63,64].Nevertheless, the abundance of SF vectors in those studies didn't show potential ''outbreaks'' [63,64] in SF abundance, thus not allowing to assess whether CL case incidence flattens out with large vector abundance.This information is necessary to properly understand the role of climate on the entomological risk for CL transmission [29,62].This goal will require new longitudinal studies on vector abundance in Panama [30] where Leishmania spp infection in the vectors is also tracked [40], in order to better understand the relationship between vector abundance and vector infection, since constant or nearly constant infections rates in vectors have different implications to understand   the role of vectors on transmission patterns.For example, if infection rate decreases with vector abundance, such a densitydependent pattern might partially explain the flattening relationship between vector abundance and cases, such as observed with bloodfeeding success by SF vectors, which decreases with density [73].Nevertheless, the flattening can also emerge by, or in synergy with, the regulation in the recruitment of susceptible hosts [29] and/or the zoonotic reservoirs, some of which might also experience population outbreaks with ENSO [74,75] which can ultimately be linked to ENSO mediated changes in the resources sustaining wildlife reservoir populations in the neotropics [76,77].Indeed, the most studied wildlife spp reservoir in Panama ´, the two toed sloth, has shown relatively large interannual fluctuations in Leishmania (Viannia) sp.infection [10].
Our results also support the major role of Lu gomezi and Lu trapidoi as dominant SF vectors of CL in Panama ´ [2,43], given their ubiquity across domestic, peridomestic and forest environments.This ubiquity has implications for the role of these species in both the transmission to humans and as bridge of pathogens across vertebrate Leishmania spp hosts and eco-epidemiological environments [78].By contrast, Lu panamensis was mainly present in forest environments, which suggests that it might not be heavily involved in domiciliary/peridomiciliary CL transmission, a possibility also put forward by recent studies on spatial patterns of human infections [43] and dogs (unpublished data), which were mainly associated, respectively, with Lu gomezi and Lu trapidoi abundance, but not Lu panamensis.Nevertheless, our inferences are limited given our focus on acknowledged dominant vectors species, which is an approach that potentially biases the identification of other vectors present at CL transmission foci, a problem that can only be solved when studying the whole SF community [79].
Finally, our study has some limitations related to the nature of the data and its countrywide geographical scale.Although our CL and SF data have a consistent quality, there is ample room to improve CL and SF surveillance.For CL surveillance, an urgent need is to standardize diagnostics across the country using sensitive and specific methods [1,26].Even if it is impossible to standardize diagnostics at the health post level, an effort should be made to estimate the error in diagnostics, as done for malaria, where all clinically diagnosed cases are confirmed at the ICGES, and quality controls on the specificity of diagnosis are equally performed [67].The non-spatial nature of our analysis precludes the identification of transmission hotspots requiring attention [80] or zones were biases in case report might be occurring [81].Similarly, the role that patterns of socio-economic inequity might have in the impacts of climate change and weather variability on CL transmission [32] cannot be estimated.Nevertheless, the relative low variability in rainfall patterns across human inhabited zones in Panama ´ [41] suggests that a country-wise analysis is a sound method to make inferences about the relationship between CL transmission, ENSO and weather patterns.The CL cases time series showed no trend  which allowed us to ignore a denominator for the cases, nevertheless we cannot assert whether the lack of trends is due to stationary population patterns in the population at risk, or if they reflect other unknown changes in the populations at risk and/ or transmission.This point could be further clarified by the establishment of health demographic surveillance systems that could both improve the understanding of disease transmission patterns and the demography of populations living in CL endemic areas.Similarly, SF monitoring can be improved and more systematically done at endemic areas.This is an issue of major importance, since, given the delay between transmission and clinical CL, an early prediction of CL epidemics will be more robust if based on the monitoring of SF abundance and Leishmania sp.infection [29].Although each SF abundance estimate came, on average, from ten trap-nights, locations were variable and in some instances SF estimates came from as few as three tree-trap nights and one location, yet these scarce records are abundant in the context of entomological surveillance for CL, and for most neglected tropical diseases.In that sense, an effort could be made to establish sentinel posts in highly endemic counties, thus rendering feasible a highly standardized estimation of SF abundance across endemic areas, where ideally vector infection is also tracked and this information used for prediction and proactive vector control [31].

Conclusion
Our data clearly supports that changes in SF abundance and CL cases reported at health facilities in Panama ´are associated with ENSO.Interannual variability in CL cases is associated with ENSO, where large epidemics follow the cold ENSO phase, while seasonal patterns are associated with temperature and rainfall variability.CL cases were positively associated with 3-month lagged Lu gomezi and Lu trapidoi abundance estimates from light traps.SF vector abundance, on average, decreased during the hot and cold ENSO phases, when compared with the normal ENSO phase, yet variability in SF was largest during the cold ENSO phase suggesting that SF population outbreaks might play a role in CL epidemics, a subject deserving further research.

Figure 1 .
Figure 1.Monthly time series data (2000-2010).(A) Cutaneous Leishmaniasis cases in the Republic of Panama ´(B) Rainfall (C) Temperature.The solid line indicates the averages and dashed lines the extremes.(D) Sea Surface Temperature 4 (El Nin ˜o 4 Index).All time series start in January 2000 and end in December 2010.In the plots colors indicate the ENSO phase, for details refer to the inset legend in panel A. doi:10.1371/journal.pntd.0003210.g001

Figure 2 .
Figure 2. Cutaneous Leishmaniasis cases (CL) seasonality.(A) Boxplots of monthly incidence.Boxes contain data within the 25 th to 75 th quantiles.Lines inside the boxes show the median of the distribution for each month.(B) Seasonal (year-long) time series.Colors indicate the ENSO phase, see inset legend for details.doi:10.1371/journal.pntd.0003210.g002

Figure 3 .
Figure 3. Correlation between selected lags of the Cutaneous Leishmaniasis cases from Republic of Panama ´time series, and climatic covariates.Time lags are indicated inside parenthesis ''()'' and SST4, Temp and Rain are, respectively, abbreviations for Sea Surface Temperature 4 (El Nin ˜o 4 Index), Temperature and Rainfall.Circle size indicates the magnitude of the association, while color indicates the sign of the correlation, a scale is presented in the right margin of the figure.doi:10.1371/journal.pntd.0003210.g003

Figure 4 .
Figure 4. Cross-wavelet coherence analysis.Coherence between (A) Cutaneous Leishmaniasis cases in the Republic of Panama ´, Leish, and Sea Surface Temperature 4, a.k.a., El Nin ˜o 4 index, SST4 (B) Leish and Rainfall, Rain (C) Leish and Average Temperature, Temp (D) Rain and SST4 (E) Temp and SST4.A cross wavelet coherence scale is presented at the bottom of the figure, which goes from zero (blue) to one (red).Red regions in the plots indicate frequencies and times for which the two series share power (i.e., variability).The cone of influence (within which results are not influenced by the edges of the data) and the significant coherent time-frequency regions (p,0.05) are indicated by solid lines.doi:10.1371/journal.pntd.0003210.g004

Figure 5 .
Figure 5. Sand Fly vector species abundance during the different ENSO phases and by eco-epidemiological environment.Abundance by ENSO phase: (A) Lutzomyia gomezi (B) Lu trapidoi (C) Lu panamensis.Abundance by eco-epidemiological environment: (D) Lu gomezi (E) Lu trapidoi (F) Lu panamensis.Panels A, B and C show data only for April, October and November where the number of trap-nights was above 30.In all panels the y-axis is in a logarithmic scale.doi:10.1371/journal.pntd.0003210.g005

Figure
Figure S2 Boxplots of monthly Cutaneous Leishmaniasis cases as function of El Nin ˜o Southern Oscillation (ENSO) phase.Boxes contain data within the 25th to 75th quantiles.Lines inside the boxes show the median of the distribution for each month.(PDF) Figure S3 Cutaneous leishmaniasis cases and climatic covariates correlation functions.(A) Cutaneous Leishmaniasis cases, Leish, from Republic of Panama ´partial autocorrelation function, PACF.Cross-Correlation Functions, CCFs, between Leish and (B) Sea Surface Temperature 4, i.e., El Nin ˜o 4 Index (C) Rainfall (D) Maximum Temperature (E) Minimum Temperature and (F) Average Temperature.Blue dashed lines indicate the 95% confidence limits for correlations that can be expected by random.(PDF) Figure S4 Cutaneous Leishmaniasis cases from Republic of Panama ´autocorrelation function (ACF).The ACF is based on monthly data from January 2000 to December 2010.Blue dashed lines indicate the 95% confidence limits for correlations that can be expected by random.(PDF) Figure S5 Cutaneous Leishmaniasis cases and Sand Fly vector abundance Cross Correlation Functions (CCFs).CCFs between the number of monthly CL cases and sand fly vector abundance/trap night/month: (A) Lutzomyia gomezi (B) Lutzomyia trapidoi (C) Lutzomyia panamensis (D) Lutzomyia gomezi and Lutzomyia trapidoi.Blue dashed lines indicate the 95% confidence limits for correlations that can be expected by random.(PDF)

Table 1 .
Parameter estimates for the best model.

Table 2 .
Generalized linear mixed Poisson rate model parameter estimates.

Table S2 Sand
Fly vector species abundance model selection.AIC and BIC stand, respectively, for Akaike and Bayesian Information Criteria.Best model selection is guided by their minimization.Best models are bolded.(PDF)