Abstract

The purpose of this study is to create distribution models of seventeen Lutzomyia species in Venezuela. Presence records were obtained from field collections over 30 years by several research teams. We used maximum entropy method for model construction based on 30 arc-second resolution environmental layers: 19 bioclimatic variables, elevation, and land cover. Three species were distributed throughout north-central Venezuelan, two restricted to northern Venezuelan coast, and three throughout the west; five were restricted mainly to the Andean and finally two species within sparse pattern. The most important variables that contributed were related to precipitation. The environmental niche model of sandflies might be a useful tool to contribute to the understanding of the ecoepidemiological complexity of the transmission dynamics of the leishmaniases.

1. Introduction

The phlebotomine sandflies are important vectors of bartonellosis [1], flaviviruses, orbiviruses, phleboviruses, and vesiculoviruses [2, 3] as well as leishmaniases [4].

Cutaneous leishmaniasis in Venezuela is endemic and focal but is distributed almost all over the country where the population is about 25,000,000. The number of cases reported by the Ministry of Health during the period 1988–2007 was 52,402 with an average incidence of 10.23 in 100,000 inhabitants; 98.67% were diagnosed as localized cutaneous leishmaniasis, 1.1% as mucocutaneous leishmaniasis, and 0.23% as diffuse leishmaniasis [5, 6]. In relation to visceral leishmaniasis (VL), for the period 1990–2009, the number of registered cases was 597 and the average incidence was 0.12 per 100,000 inhabitants (unpublished files of the Department of Informatics, Instituto de Biomedicina, Caracas).

Feliciangeli [7] reported about 100 species of phlebotomine sandflies known in Venezuela, 30 of which are anthropophilic and eight are recognized as suspected or proven vectors of the leishmaniases. Following Young and Duncan’s classification [8], we list here the proven vectors of cutaneous leishmaniasis (CL): Lutzomyia ovallesi (Ortíz, 1952), L. gomezi (Nitzulescu, 1931), L. panamensis (Shannon, 1926), L. youngi Feliciangeli & Murillo, 1985, L. spinicrassa Morales, Osorno-Mesa, Osorno, and Muñoz de Hoyos, 1969, and L. migonei (França, 1920) and the proven vectors of visceral leishmaniasis (VL): L. longipalpis s.l. (Lutz & Neiva, 1912), L. pseudolongipalpis Arrivillaga & Feliciangeli, 2001, and L. evansi (Nuñez-Tovar, 1924).

Earlier studies have demonstrated that historical field and laboratory data regarding the ecology and distribution of vector-borne diseases can provide the basis for development of spatial-temporal environmental risk models by the use of standardized analysis of Geographical Information Systems, remote sensing data from earth-observing satellites, and ecological niche modeling [913]. Ecological niche modeling (ENM) uses record data in conjunction with environmental data to develop models of habitat range for a given organism [14, 15]. Some authors stated that ENM today plays a central role in many areas of ecology, conservation, and evolutionary biology, both because they can fill gaps in knowledge and allow a better estimate of multiple components of species diversity [1618]. This technique has been used in modeling distribution of diseases, such as dengue, malaria, and West Nile and yellow fever vectors, such as Aedes Meigen 1818 [19], Anopheles Meigen 1818 [2022], Culex Linnaeus 1758 [23], and Haemagogus Williston 1896 [24], respectively. In addition, ENM has been used to examine the potential distribution of Lutzomyia França 1924 in America [2530] and Phlebotomus Rondani 1840 in western Asia [31]. Recently, Foley et al. [32] proposed a new map service that allows free public online access to global sandfly collection records and habitat suitability models.

In this study, we use ENM to develop distribution models for seventeen species of phlebotomine sandflies (sixteen from Lutzomyia and one of Brumptomyia Franca & Parrot, 1921) in Venezuela. Using these models, we attempt to identify environmental factors which influence the distribution of these species, in order to contribute to leading to a more sophisticated risk analysis for leishmaniasis in Venezuela.

2. Material and Methods

2.1. Species Records

As part of a sandfly biogeographical study, presence data for the species were taken from Venezuelan collections performed by the Centro Nacional de Referencia de Flebótomos y otros Vectores (CNRFV), BIOMED, Universidad de Carabobo, Maracay, Venezuela, during 30 years, and bibliographic revision in scientific literature dating from 1970 through the present compiled by one of the authors (María Dora Feliciangeli). This study was restricted only to species with ten or more presence records and included six species of medical importance. Small sample sizes pose challenges to any statistical analysis and result in decreased predictive potential when compared to models developed with more occurrences. Hernandez et al. [33] evaluate several ENM algorithms and found that model accuracy increased with larger sample sizes for all modeling methods across the taxa tested. Nonetheless, useful models were produced with as few as 5–10 positive observations. Because of that, we included the sandflies species within small records sizes, in order to evaluate the accuracy of obtained models in such taxa. We compiled a country database (from 21 administrative states, excluding Monagas, Vargas, and Delta Amacuro) with 681 specimen records (98% represents CNRFV site collections) of 17 species: Brumptomyia beaupertuyi (Ortiz, 1954) (), Lutzomyia antunesi (Coutinho, 1939) (24), L. atroclavata (Knab, 1913) (80), L. cayennensis cayennensis (Floch & Abonnenc, 1941) (81), L. dubitans (Sherlock, 1962) (38), L. evansi (43), L. gomezi (92), L. lichyi (Floch & Abonnenc, 1950) (49), L. longipalpis s.l. (67), L. micropyga (Mangabeira, 1942) (19), L. migonei (23), L. nuneztovari (Ortiz, 1954) (12), L. olmeca bicolor Fairchild & Theodor, 1971 (12), L. ovallesi (47), L. panamensis (33), L. pilosa (Damasceno & Causey, 1944) (14), and L. punctigeniculata (Floch & Abonnenc, 1944) (29). Lutzomyia pseudolongipalpis was not included in this study because its known distribution is only restricted to three close localities in Lara state. The database coordinates were converted to the decimal degrees format, using GPS coordinates, 1 : 100,000 regional maps, and geographical gazetteers, and for database standardization protocols we followed [34, 35].

2.2. Environmental Layers

We used 19 bioclimatic variables [36] derived from climatic data from the 1950–2000 period: annual mean temperature (BIO1), mean diurnal range (BIO2), isothermality (BIO3), temperature seasonality (BIO4), maximum temperature of warmest month (BIO5), minimum temperature of coldest month (BIO6), temperature annual range (BIO7), mean temperature of wettest quarter (BIO8), mean temperature of driest quarter (BIO9), mean temperature of warmest quarter (BIO10), mean temperature of coldest quarter (BIO11), annual precipitation (BIO12), precipitation of wettest month (BIO13), precipitation of driest month (BIO14), precipitation seasonality (BIO15), precipitation of wettest quarter (BIO16), precipitation of driest quarter (BIO17), precipitation of warmest quarter (BIO18), and precipitation of coldest quarter (BIO19) and one topographic variable [37]: elevation (ELEV) and one ecological [38]: fifteen land cover classes (LAND). All variables were reduced to a grid resolution of 30 arc-seconds or 0.008333° (approximately 1 Km2) for the analysis.

2.3. Model Building and Evaluation

We modeled the species distribution using maximum entropy method [17, 39] with MaxEnt 2.3 [17]. The analysis was performed using the default parameters: maximum iterations to 500 and using convergence threshold in [40, 41]. Seventy-five percent of the data points for each species were randomly selected as training points and used in model building. The remaining 25% of the records were test points and used in model validation. Duplicate presence records were removed by the program prior to model development. The model output was set to logistic: it is an attempt to get as close as we can to estimate the relative habitat suitability that the species is present, given the environment. The program returns an estimated presence probability for a given location between cell values of 0 (no probability of species presence) and 1 (species is certain to be present); we based our predictions maps on a probability of 0.70 (values above this threshold were considered presence and they represented higher suitability). The model was evaluated using the area under the curve (AUC) of the receiver operating characteristic (ROC); this technique computes the total area under the curve created by plotting sensitivity against the fractional predicted area for the species. Also, a one-tailed binomial test was calculated with the null hypothesis being that the model does not predict the test points better than random [17, 31, 39, 42].

3. Results

Table 1 represents the accuracy of the niche models in the sandflies species. Of these 17 niche models all possessed an AUC greater than 0.80, with 15 of them being statistically significant () as indicated by the binomial test. Niche models were produced for Lutzomyia antunesi, L. atroclavata, L. c. cayennensis, L. dubitans, L. evansi, L. gomezi, L. lichyi, L. longipalpis s.l., L. micropyga, L. migonei, L. nuneztovari, L. olmeca bicolor, L. ovallesi, L. panamensis, and L. punctigeniculata. Figures 13 present maps of the predicted geographical distributions of sandflies species. Darker regions are those of high relative probability (>0.70) of occurrence while lighter areas are those in which the relative probability of occurrence is low (<0.70). The species present a variety of distributional patterns across Venezuela: Lutzomyia evansi and L. panamensis were distributed throughout north-central Venezuela; L. cayennensis and L. olmeca bicolor were restricted to northern Venezuelan coast; L. nuneztovari, L. longipalpis s.l., L. punctigeniculata, and L. gomezi were distributed throughout the west, L. micropyga, L. migonei, L. lichyi, and L. ovallesi were restricted mainly to the Andean region; L. atroclavata is distributed throughout northern Venezuela; finally L. dubitans and L. antunesi were present in a sparse pattern in the country.

Table 2 informs about the environmental variable contribution based on the Jackknife test; first, the AUC of each niche model was calculated using each of the environmental variables individually. Those variables that resulted in the highest AUC were interpreted as those which possessed the most information regarding the niche of a species. Second, the AUC of each niche model was calculated after omitting each of the environmental variables one at a time. In our study, the variables that increment the AUC value when included separately were precipitation of wettest quarter (BIO16), precipitation of wettest month (BIO13), annual precipitation (BIO12), and precipitation of coldest quarter (BIO19). The remaining variables showed smallest importance in the AUC contribution: precipitation of driest month (BIO14), elevation (ELEV), temperature annual range (BIO7), land cover (LAND), and so forth. However, some variables decrease the AUC when they are omitted from the models: isothermality (BIO3), BIO12, mean diurnal range (BIO2), LAND, and ELEV. For additional MaxEnt output, see the Supplementary Material available online at http://dx.doi.org/10.1155/2015/108306, within the omission and predicted area and Jackknife regularized trained gain figures for each species.

4. Discussion

The predicted distributions of Venezuela sandflies produce useful models within species that differ from 12 records (L. olmeca bicolor and L. nuneztovari) to 92 records (L. gomezi); however, in L. pilosa (14 records) and B. beaupertuyi (18) the models performance was poor. Some studies have demonstrated deterioration in predictive performance as sample sizes are decreased [26, 42]. In others, accuracy models were produced with few records maybe due to ecologically specialized species with smallest geographic extent of occurrence and very low tolerance [33]. Our records showed L. olmeca bicolor to be restricted to lowland deciduous forests (100–250 masl), L. nuneztovari to montane evergreen forests (500–1900 masl), L. gomezi to various land cover types and wide elevation range, L. pilosa to montane evergreen forest, agriculture types and wide elevation range, and finally B. beaupertuyi to various land cover types and wide elevation range. Our results are similar to ten Lutzomyia species potential distribution maps generated by the global geospatially referenced clearinghouse for sandfly disease vector species collection records (http://www.sandflymap.org/; access: 1/19/2015). However, the differences may be due to the environmental variables used to build the model and the number of records (presences) in each species. In Venezuela, there are only 40 records for seventeen species in the SandflyMap database, in contrast to the 681 records in our study.

In strict epidemiological sense, the Venezuelan Andes (Trujillo, Merida, Tachira states) and the Lara state showed to be the most important area with CL transmission [4346]. Feliciangeli et al. [47] stated that the distribution of Lutzomyia evansi (VL vector) and L. ovallesi almost overlaps the Andean lowland forests as well in the Northern Coast Cordillera; L. ovallesi is restricted to 0–800 masl and L. nuneztovari to areas close or to above 1,000 masl. In a review of the subgenus Micropygomyia of Lutzomyia, Feliciangeli [48] showed that the geographical distribution of L. micropyga to the Andean states may represent the natural continuity of the Colombian distribution; L. c. cayennensis occupies wide range of life zones (dry to moist tropical forests), whereas L. atroclavata is normally found in the lowlands (0–400 masl).

Rotureau [49] stated that Lutzomyia vectors are generally abundant year round in the Amazon. Nevertheless, the period encompassing the end of the dry season and the beginning of the rainy season is more propitious for their development. Peterson and Shaw [25] proposed a niche model for cutaneous leishmaniasis vectors from southern Brazil: they found that L. migonei had the broadest ecological distribution, and Lutzomyia whitmani (Antunes & Coutinho 1939) tended to be the most tropical in distribution (high precipitation, high temperature, wet climates, and high vapour pressure). More recently, a niche model for Leishmaniasis in north America suggests that the most important environmental variables were mean temperatures of the wettest and warmest quarters in Lutzomyia anthophora (Addis, 1945) and annual mean temperature and the minimum temperature of the coldest month in L. diabolica (Hall, 1936) [27]. Colacicco-Mayhugh et al. [31] determined the niche model in Phlebotomus alexandri Sinton 1928 and P. papatasi (Scopoli, 1786), using land cover classes derived from Normalized Difference Vegetation Index (NVDI) and 19 bioclimatic variables. The variables with high probabilities presence were urban, field/woody savanna, and woody savannah coverage, while the remaining variables have very modest training gains. They consider that this may be partly due to sampling bias, as collections of phlebotomine sand flies tend to be associated with research related to human leishmaniasis and anthropogenic biomes. Almeida et al. [50] observed in Mato Grosso state (Brazil) similar predicted areas between niche models for L. longipalpis and L. cruzi (Mangabeira 1938) when NVDI land cover classes and climatic variables were included together or in a separate analysis. Finally, they reported high abundance of L. longipalpis in the rainy season, and in L. cruzi followed moderate to regular precipitation and minimum temperature elevation. Quintana et al. [29] showed that precipitation seasonality and precipitation in the warmest quarter were the variables that most contributed to the ENM of Lutzomyia neivai and Lutzomyia migonei in Argentina. Our results suggest that precipitation of driest quarter and precipitation of driest month were the most important variables in L. migonei model prediction. In Colombia, [30] noted that distributions of L. longipalpis and L. evansi are strongly correlated with precipitation; their breeding success depends on the availability of rich humid soils, and precipitation during dry months, when the river flow can decrease, can become a limiting factor. These findings are similar to ours; annual precipitation is the most important variable for these species that contributes to the model construction.

The VL vectors L. evansi and L. longipalpis s.l. showed to be sympatric in northern Venezuela [49, 51]. Then L. panamensis and L. evansi were also reported to be focused in Caribbean coast of eastern Venezuela where cutaneous and visceral leishmaniasis coexist [52]; Lutzomyia gomezi and L. ovallesi in the north-central Venezuela are the main anthropophilic species implicated as vectors of this disease [49, 53]. Feliciangeli [7] reported the geographic distribution of eight leishmaniasis vectors; of these, five species showed the widest distribution within fourteen to twenty administrative states, and the remaining species were restricted to one to four states. Our results are congruent with this author, showing wide predicted distributions for L. longipalpis s.l., L. evansi, L. ovallesi, and L. gomezi; however, the potential distribution in L. panamensis was limited to fewer states.

Recently, in Mexico, [27] estimate the probable distribution of all sandfly species with potential medical importance for leishmaniasis and related it to the transmission areas.

In conclusion, although we are aware about the limitations of using secondary data collected with different methods across many years, our results show that the ENM might be a useful tool to contribute to the understanding of the ecoepidemiological complexity of the transmission dynamics of the leishmaniases. To reach this goal, further studies are needed which must include epidemiological variables, such as the distribution of CL/VL cases and reservoirs.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

The authors thank the Venezuelan Ministry of Science and Technology, for supporting by means of the Projects MISION CIENCIA: 2008000911-2 and 2008000911-4.

Supplementary Materials

The following figures shows the omission rate and the jackknife test of variable importance for each sandfly species. The omission rate vs the predicted area as a function of the cumulative threshold. The omission rate is calculated both on the training presence records (75%), and on the test records (25%) for each sandfly species. The omission rate is the fraction of the test localities that fall into pixels not predicted as suitable for the sandfly species, and the proportional predicted area, which is the fraction of all the pixels that are predicted as suitable for the species. In order to determine which variables contribute most to the model development, was calculated the jackknife tests of variable importance. The jackknife procedure produces three different types of models: (1) models created with one variable at a time excluded and all other variables included, (2) models created with only one variable included, and (3) a model created with all variables. Variables that are most important to model development are those that decrease the training gain when removed from the model and show gain when the model is developed with only one variable.

  1. Supplementary Material