Genetic diversity, phylogeography and molecular clock of the Lutzomyia longipalpis complex (Diptera: Psychodidae)

Background The Lutzomyia longipalpis complex has a wide but discontinuous distribution in Latin America, extending throughout the Neotropical realm between Mexico and northern Argentina and Uruguay. In the Americas, this sandfly is the main vector of Leishmania infantum, the parasite responsible for Visceral Leishmaniasis (VL). The Lu. longipalpis complex consists of at least four sibling species, however, there is no current consensus on the number of haplogroups, or on their divergence. Particularly in Argentina, there have been few genetic analyses of Lu. longipalpis, despite its southern expansion and recent colonization of urban environments. The aim of this study was to analyze the genetic diversity and structure of Lu. longipalpis from Argentina, and to integrate these data to re-evaluate the phylogeography of the Lu. longipalpis complex using mitochondrial markers at a Latin American scale. Methodology/Principal findings Genetic diversity was estimated from six sites in Argentina, using a fragment of the ND4 and the 3´ extreme of the cyt b genes. Greatest genetic diversity was found in Tartagal, Santo Tomé and San Ignacio. There was high genetic differentiation of Lu. longipalpis in Argentina using both markers: ND4 (FST = 0.452, p < 0.0001), cyt b (FST = 0.201, p < 0.0001). Genetic and spatial Geneland analyses reveal the existence of two primary genetic clusters in Argentina, cluster 1: Tartagal, Santo Tomé, and San Ignacio; cluster 2: Puerto Iguazú, Clorinda, and Corrientes city. Phylogeographic analyses using ND4 and cyt b gene sequences available in GenBank from diverse geographic sites suggest greater divergence than previously reported. At least eight haplogroups (three of these identified in Argentina), each separated by multiple mutational steps using the ND4, are differentiated across the Neotropical realm. The divergence of the Lu. longipalpis complex from its most recent common ancestor (MRCA) was estimated to have occurred 0.70 MYA (95% HPD interval = 0.48–0.99 MYA). Conclusions/Significance This study provides new evidence supporting two Lu. longipalpis genetic clusters and three of the total eight haplogroups circulating in Argentina. There was a high level of phylogeographic divergence among the eight haplogroups of the Lu. longipalpis complex across the Neotropical realm. These findings suggest the need to analyze vector competence, among other parameters intrinsic to a zoonosis, according to vector haplogroup, and to consider these in the design and surveillance of vector and transmission control strategies.


Introduction
Visceral leishmaniasis (VL) is a parasitic disease caused in the American continent, by Leishmania infantum (syn. Le. chagasi). VL has an estimated global incidence of 500,000 cases and 59,000 deaths per year [1]. Human cases of VL in Latin America are reported in 12 countries, while Argentina reports an expanding epidemiological scenario [2]. The first autochthonous VL cases in both humans and canines in Argentina were reported in 2006, in the city of Posadas, Misiones [3]. Currently, the country has reported 106 human cases of VL, with 7.7% mortality, between 2006 and 2016 [4]. Misiones is the province with most human VL cases; the age group with highest incidence is children 0-15 yrs [5]. In Argentina, the main vector of Le. infantum is Lutzomyia longipalpis (Lutz & Neiva) [3]. This phlebotomine was first reported in 1951 in Candelaria, and later in 2000, in Corpus, both in Misiones [6]; there were no VL cases reported at either time. The first autochthonous VL human cases with parasite presence confirmed were reported from Posadas, in 2006 [3]. At present, the vector has dispersed to southern Argentina while also expanding in the north. Currently, Lu. longipalpis is also reported in Corrientes, Entre Ríos, Chaco, Formosa, and Salta provinces [5,[7][8][9][10].
There are few studies on the genetic status of Lu. longipalpis in Argentina, although it is clear that the species is expanding southward and is colonizing urban environments [52]. Salomón et al. [33] using the "per" gene, found high genetic differentiation between Lu. longipalpis from Posadas (Misiones province, Argentina) and the north and southeast regions of Brazil, suggesting they might represent another sibling species. In contrast, males from Argentina secrete the same male pheromone as those from Paraguay [53] and other populations from Brazil [27,29,54]. It is therefore important to be able to identify and differentiate haplogroups, since they may differ in importance as diseases vectors. Haplogroups may differ in terms of vector competence, vector capacity, or other epidemiological aspects [55], or the clinical expression of disease, based on maxadilan salivary content [31]. Although at least four sibling species have been proposed for the Lu. longipalpis complex [16,18,19,56,57], there is insufficient evidence currently for a general consensus on the exact number of haplogroups, or their divergence. The aim of this study was to analyze the genetic diversity and structure of Lu. longipalpis from Argentina, and to integrate these data to re-evaluate the phylogeography of the Lu. longipalpis complex using mitochondrial markers at a Latin America scale.

Study area and taxonomic identification
Lutzomyia longipalpis was collected in six sites in Argentina by the Leishmaniasis Research Network (REDILA). Four of these localities have reported human VL cases (Puerto Iguazú, San Ignacio, Santo Tomé, and Corrientes city), one has reported only canine cases (Clorinda), and one has not reported any VL cases to date (Tartagal) (Fig 1, S1 Table). Clorinda and Corrientes city are located in the humid Chaco ecoregion, which has a warm subtropical climate. Puerto Iguazú and San Ignacio are located in the Paranense forest ecoregion with a humid subtropical climate, and Santo Tomé is found in grassland and forest ecoregion also in a humid subtropical climate. Tartagal is located in the Yunga forest ecoregion with a subhumid warm climate [58]. Sandfly collections were carried out using REDILA-BL light traps [59], and specimens were preserved in 90% ethanol and transported to the laboratory. The last three Sampling sites for Argentinian specimens for which there are sequences for both gene fragments (triangles), for those sequences analyzed from GenBank for ND4 [16,42] (circles), and those from GenBank for cyt b sequences [20,37,40,41] (square). This map was generated using QGIS v.2.6.1 [61]. segments of males were dissected for taxonomic identification, while the rest of the insects were used for molecular analyses. The dissected segments were clarified with lacto-phenol and mounted on slides, and identified using optical microscopy using the Galati taxonomic key [60].

DNA extraction, amplification and sequencing
Genomic DNA was extracted from individual sandflies using DNA Puriprep-S kit Highway (Inbio, Argentina), according to manufacturer's instructions. The ND4 fragment was amplified using primers reported by Soto et al. [16]. PCR reactions were carried out in 50 μl volumes, containing 25 μl of Gotaq green Master Mix (Promega, USA), 50 pmol/μl of each primer, 2.5 mM of MgCl 2 and 50 ng of template DNA. Amplification conditions were: 95˚C for 5 min, followed by 35 cycles of 95˚C (30 sec), 48˚C (45 sec), 72˚C (1 min) and a final elongation step at 72˚C for 10 min. The 3´end of the cyt b fragment was amplified using primers and amplification conditions reported by Hodgkinson et al. [40,41] and the PCR was performed according to Pech May et al. [43]. PCR products were separated using 1.5% agarose gel electrophoresis, stained with 0.5 μg/ml Sybr Green (Invitrogen, USA), and visualized under UV light. PCR products were purified by HW DNA Puriprep-GP kit (Inbio, Argentina), and sequenced using both forward and reverse primers in a CEQ 2000XL automated sequencer (Beckman Coulter, USA).

Genetic analyses of Argentinian populations
Forward and reverse sequences from all samples were used to generate consensus sequences and were manually aligned and edited using MEGA v.7 [62]. All haplotypes were deposited in GenBank (accession numbers for the ND4 fragment: MH166358-MH166392; accession numbers for the cyt b fragment: MH166339-MH166356). Intra-population and global genetic diversity was analyzed using the number of mutations (η), the number of segregating sites (S), the number of unique sites (Su), the mean number of pairwise differences (k), haplotype number (Nh), haplotype diversity (Hd), nucleotide diversity (π), and the nucleotide polymorphism index (θ), using DnaSP v.5.10 [63]. Neutrality tests based on Tajima's D and Fu's Fs [64,65] were calculated based on segregating sites, using the same software. In addition, the mismatch distribution test was also analyzed using ARLEQUIN v.3.5 [66]. The goodness-of-fit of the observed data to a simulated expansion model was tested using both Harpending's raggedness index (r) [67] and the sum of squares deviations (SSD), using 10,000 replicates. Molecular variance (AMOVA) was used to evaluate population genetic differentiation, using 10,000 random permutations in ARLEQUIN v.3.5. The p values of pairwise F ST were adjusted using the Holm-Bonferroni sequential correction "p´" [68,69]. Association between genetic and geographic distance was analyzed using a Mantel test [70], implemented in a trial version of XLTAT (https://www.xlstat.com). Linear geographic distances between sites were estimated using QGIS v.2.6.1 [61].
ND4 and cyt b fragments were concatenated and analyzed using GENELAND package in R to infer the number of Lu. longipalpis genetic clusters among Argentinian sites [71][72][73][74]. This program uses a spatial statistical model and Markov chain Monte Carlo sampling with GPS coordinates to estimate the number of populations or genetic clusters (K). Preliminarily, we estimated the number of K from 1 to 10, using 10 million MCMC iterations and 1,000 thinnings. Five independent runs with fixed "K" were run (to avoid ghost populations), assuming an uncorrelated allelic frequency and spatial model. For each run, the posterior probability (PP) of subpopulation membership was computed for each pixel of the spatial domain (100 x 100 pixels), using a burn-in of 1,000 iterations. Variation among and within populations of genetic clusters was analyzed using 10,000 random permutations for the AMOVA in ARLE-QUIN v.3.5.
Divergence times among haplogroups were estimated using BEAST v.1.8.4 [78]. We calibrated using the divergence rate reported by Esseghir et al. [36], and corrected as suggested by Ho et al. [79], assuming 0.105 substitution/site/million years. An uncorrelated lognormal relaxed-clock model was used to allow rate variation among branches and the coalescent exponential growth option. Convergence of the MCMC chains was checked using TRACER v.1.6, with 200 as the minimum effective sample size. The length of the runs were: a) 100 million generations for the ND4 fragment, and b) 500 million generations for the cyt b fragment; sampling was every 10,000 and 50,000 generations, respectively. The first 10% of the samples were discarded as burn-in and a maximum clade credibility tree, reflecting divergence times and their 95% highest posterior densities, was estimated using TreeAnotator of the BEAST package and visualized using FigTree v.1.4.
Nucleotide divergence (Da) for the haplogroups was estimated based on the number of net nucleotide substitutions, using the Jukes and Cantor (JC) correction [80,81] with DnaSP v.5.1.0. The pairwise F ST comparison between haplogroups was performed using 10,000 random permutations in ARLEQUIN v.3.5. The p values were adjusted using a Holm-Bonferroni sequential correction "p´" [68,69]. The single sequence from Venezuela was not included in these latter analyses.
The global genetic differentiation of Lu. longipalpis from Argentina using the ND4 fragment is high (F ST = 0.452, p < 0.0001), although most of the genetic differentiation was within populations (54.8%). Greatest genetic difference was between San Ignacio and Corrientes city, with an F ST of 0.86 (p < 0.0001) ( Table 2). There was no evidence for genetic isolation associated with geographic distance (r = 0.24, p = 0.37).

Lu. longipalpis genetic diversity in Argentina using the cyt b fragment
Seventy-six specimens of Lu. longipalpis were sequenced from all six sites in Argentina, using the 261 bp fragment of the cyt b gene (Table 1B). This fragment had 20 polymorphic nucleotides (7.66%), and the A-T composition was 76.3%. The cyt b fragment identified 18 haplotypes, ranging from 2 to 12 per site. Globally, haplotype diversity was high (Hd ± SD = 0.772 ± 0.036), similar to that for the ND4 fragment, as were low nucleotide diversity and nucleotide polymorphism index (π ± SD = 0.008 ± 0.0009 and θ ± SD = 0.015 ± 0.003, respectively). San Ignacio, Santo Tomé and Tartagal had highest cyt b diversity (Table 1B) Table 2). Again, there was no evidence for genetic isolation associated with geographic distance (r = 0.49, p = 0.06).

Phylogeographic analysis using the ND4 fragment
A total of 91 haplotypes of Lu. longipalpis reported from seven Latin American countries were included in ND4 phylogeographic analyses (Guatemala, Honduras, Costa Rica, Venezuela, Colombia, Brazil and Argentina) (S1 Table). The Median-joining network identifies at least  (Fig 3, for details of haplotypes see S1 Table).
The BI for ND4 phylogeographic analysis was constructed using the HKY + G model as the most appropriate for the data (-lnL = 2913.8833; BIC = 7035.9464) with gamma of 0.1510. The tree revealed eight haplogroups with maximum support for two major clades (PP = 1.0; Fig 4). One of the major clades is composed of Ar1 (Argentina), Ar-Bra (Argentina-Brazil), Ar2 . Net nucleotide substitutions (Da) between haplogroups ranged from 1.2% to 7.5%, while the pairwise F ST between haplogroups was high, with genetic differentiation ranging from 0.468 to 0.966 (Table 3).

Phylogeographic analysis using the cyt b fragment
Phylogeographic analyses using the cyt b fragment included 58 haplotypes of Lu. longipalpis from Venezuela, Brazil, and Argentina (S1 Table). The Median-joining network, had few mutational steps separating haplotypes, identifying five haplogroups (BI was also taken into account for haplogroup definition; Figs 5 and 6). The Bra haplogroup has a difference of two mutational steps with H46/Tartagal/Argentina (Ar1 haplogroup), while the Ar1 haplogroup  https://doi.org/10.1371/journal.pntd.0006614.g004 has haplotypes only from Tartagal, Santo Tomé, and San Ignacio. The Ar2 haplogroup has cyt b haplotypes from all six Argentinian sites, the most frequent haplotype being H39 (present in five sites), although this has only a two mutational step difference with H5 (Brazil). The Ar-Bra haplogroup has haplotypes from Tartagal, Santo Tomé and from Juazeiro (Brazil). The single haplotype from Altagracia de Orituco, Venezuela, has a 9 mutational step difference with other frequent haplotypes shared with four Argentinian populations (H40) (Fig 5). The BI from cyt b sequences was constructed using the HKY + G model as the most appropriate for the data (-lnL = 1073.7679; BIC = 2793.0201) with gamma 0.1910. Two major clades were identified, with PP ranging from 0.6 to 1.0: a) one with the Ar1, Ar2, Ar-Bra and Bra haplogroups, and b) the clade with the single haplotype from Venezuela. One haplotype from Juazeiro (H5/Brazil) grouped with Argentinian haplotypes into the Ar-Bra haplogroup (Fig 6). Divergence of the complete Lu. longipalpis complex from the MRCA was estimated at 0.45 MYA using the cyt b fragment (95% HPD interval = 0.22-0.83 MYA). The MRCA for the Argentinian and Brazilian haplogroups was estimated at 0.28 MYA (95% HPD interval = 0.15-0.45 MYA). Net nucleotide substitutions (Da) between haplogroups varied from 0.4%-4.7%, while the pairwise F ST between haplogroups was high, genetic differentiation ranging from 0.381 to 0.521 (Table 4).

Discussion
The two mitochondrial fragments from the ND4 and cyt b genes identified high overall haplotype diversity, but relatively low nucleotide diversity and nucleotide polymorphism for Lu. longipalpis in Argentinian populations. Similar indices for the ND4 gene were reported by Soto et al. [16] in Lu. longipalpis populations from Honduras, Central America (Honduras + Guatemala) and Colombia. Likewise, Coutinho-Abreu et al. [20] using the 3´of the cyt b Genetic diversity and phylogeography of the Lutzomyia longipalpis complex  Genetic diversity and phylogeography of the Lutzomyia longipalpis complex gene from nine Brazilian Lu. longipalpis populations reported similarly low nucleotide diversity, but higher haplotype diversity. Variability in Lu. longipalpis is similar to that within other species having wide distributions such as P. papatasi in the Old World using the ND4 gene [38]. A greater number of ND4 haplotypes were found (Nh = 35), than that for the 3´of cyt b gene (Nh = 18) in the present study. The populations with highest genetic diversity, using both markers, were Tartagal, Santo Tomé, and San Ignacio. Previous studies have demonstrated that in general, older established populations have higher genetic diversity, which could be related in part to a relatively constant population size [82]. Neutrality indices were not significant for the three populations, which coincides with little change in population size. However, Tartagal had contradictory results, since the cyt b neutrality test (Fs = -5.598) and its mismatch distribution were consistent with an excess of low-frequency haplotypes, both characteristic of relatively recent population expansion, or of selective sweep/hitch-hiking [83]. In contrast, neutrality tests for the ND4 in Tartagal were not significant and the mismatch distribution was multimodal, which is consistent with demographic equilibrium [84]. Both fragment mismatch distributions from Santo Tomé and San Ignacio indicate that few individuals share each haplotype, which along with high haplotype diversity and a large number of unique haplotypes, may indicate that their populations were historically established. In comparison, Corrientes city, Puerto Iguazú, and Clorinda had lower genetic diversity, which may be the result of a reduced effective population size and capacity for dispersal (causing loss of genetic diversity and increase of differentiation among the populations). These two features would favor genetic drift, as was proposed for Lu. cruciata [43], Lu. anduzei [44], and Lu. olmeca olmeca [45]. Additionally, populations with low haplotype and nucleotide diversity may have recently experienced prolonged or severe bottlenecks [85]. Results from both gene fragments are contradictory from Corrientes city, the neutrality test with cyt b suggesting an important population size reduction (D = 2.083, p < 0.05, and positive Fs), but the ND4 indicating population expansion (significant and negative Fs and unimodal mismatch distribution). A greater sampling effort will be required for Puerto Iguazú, Corrientes city, and Clorinda, from which too few haplotypes were identified for analyses and comparisons. These populations should be continuously monitored to analyze the impacts of environmental factors (abiotic, biotic, fragmentation) and healthcare or agricultural interventions (insecticides) over time. Landscape modification or fragmentation caused by anthropogenic actions (e.g. agricultural practices, deforestation) has been shown to be associated with the loss of genetic diversity in Lu. cruciata [43], Lu. gomezi [51], Lu. anduzei [44], and P. papatasi [39].
Argentinian Lu. longipalpis populations have high global genetic differentiation, with both markers, perhaps associated with the species´reduced dispersal capacity. Sandfly dispersal depends on several factors, including average flight distance, wind speed and variation, and distance to resources [86]. It also depends on landscape heterogeneity (abiotic and biotic conditions), which influence female host, resting, mating, and egg laying site selection [87]. Phlebotomines are generally poor fliers, with movement restricted to short, flight-assisted hopping. Lu. longipalpis flight range is from 1 m to 500 m, but no more than 1,000 m around their breeding sites [86]. This flight range limitation would lead to rapid local population differentiation, due to genetic drift or to population fragmentation [43,44,51]. In Argentina, two specific Lu. longipalpis genetic clusters were identified. Given shared haplotypes in Tartagal and Santo Tomé, cluster 1 may have occurred via a colonization event from Brazil (Ar-Bra haplogroup) or from the Bolivia-Brazil-Paraguay Gran Chaco eco-region, which is associated with sporadic cases of rural human VL [6]. Curiously, there are no sequences of Lu. longipalpis from these latter countries in GenBank. In the case of Santo Tomé, the most likely route of entry to Argentina is São Borja, Brazil, where an outbreak of VL was reported in 2008 [88]. A similar situation could have occurred with rapid dispersal of Lu. longipalpis-VL and urban colonization in Campo Grande (Mato Grosso do Sul-Brazil). Subsequent colonization may have occurred via Paraguay and the Argentinian border areas of Clorinda and Posadas. In the present study, we detected genetic similarities among the populations of Clorinda, Corrientes city and Puerto Iguazú (these three populations share several haplotypes and low genetic differentiation), which is in agreement with rapid invasive scenarios of urban Lu. longipalpis-VL [52,89].
Bayesian inference using the cyt b gene was not well supported, probably due to its low substitution rate [44]. However, too few haplogroups may have been analyzed (Argentinian, Brazilian and Venezuelan), and hence analysis of a larger sample set is recommended. The ND4 gene inference, in contrast, strongly supports monophyly for the Lu. longipalpis complex. Currently, the exact number of sibling species in the complex remains tentative, due to its broad distribution and absence of sufficient analyses. Soto et al. [16] identified four clades using the ND4 gene (Brazil, Central America, and laboratory colony population from Colombia and Venezuela), and similarly, Arrivillaga et al. [18,19] using the COI gene also identified four clades (Laran: Venezuela; cis-Andean: Venezuela, Colombia, northern Brazil; trans-Andean: Venezuela, Colombia, Central America; and Brazilian: Brazil). Analyses of the ND4 sequences from the present study, in addition to those reported by Soto et al. [16] and Sonoda [42], identified eight haplogroups: Ar1 (Argentina), Ar2 (Argentina), Ar-Bra (Argentina, Brazil), Bra (Brazil), CA (Guatemala, Honduras, Costa Rica), Col1 (Colombia), Col2 (Colombia) and Ven (Venezuela). The latter haplogroup is interpreted with caution since only one sequence was analyzed despite the MRCA support with 1.0 PP. These eight haplogroups are highly distant by a minimum of several mutational steps, and they have high pairwise F ST and nucleotide divergence, suggesting negligible gene exchange, similar to that suggested in a previous study [16]. This high genetic differentiation among haplogroups may potentially be due to vicariance, and/or climatic tolerance limits, in addition to low dispersal capacity. Geographic and/or climatic barriers have already been associated with diversification not just for Lu. longipalpis [16,20,42], but also for Lu. whitmani [48,49], Lu. cruciata [43], Lu. gomezi [51], Lu. umbratillis [90], and P. papatasi [91].
The present study provides evidence that there are at least three Lu. longipalpis haplogroups in Argentina. Salomón et al. [33] were the first to provide evidence that Argentinian populations may be sibling species to those reported from the northeast and southeast of Brazil [29,92,93]. The present study documents significantly high genetic differentiation between the Argentinian and Brazilian haplogroups, despite the Ar-Bra haplogroup, which groups haplotypes from Argentina (Tartagal and Santo Tomé) and Brazil (Jacobina, Lapinha caves and Juazeiro). This genetic difference contrasts with similar complex-level compounds produced by Lu. longipalpis from Posadas, Argentina [33], Asunción, Paraguay [53], and many populations from Brazil including Lapinha [54].
Divergence of the Lu. longipalpis complex from its MRCA occurred approximately 0.70 MYA, and resulted in two principal clades, one located east and south of the Amazon basin, giving rise to the principal South American (SA) haplogroups (Brazil and Argentina), and another north from northern South America (NSA), through the Mesoamerican corridor (CA) to Mexico. These two clades are similar to those reported by Arrivillaga et al. [18,19], although in contrast to those reported by Soto et al. [16]. Divergence times indicate that this latter cluster from Mesoamerica and NSA, was the first to diverge (0.58 MYA), significantly earlier than the cluster from east and south SA (0.36 MYA). Secondary divergence between the Col2 and the Ven haplogroups (0.25 MYA) and between Ar2 and Ar-Bra/Ar1 haplogroups (0.22 MYA), occurred on similar time scales. Diversification of all haplogroups occurred after the middle Pleistocene, probably during inter-glacial periods, when landscape fragmentation probably provoked diversification hotspots which conserved high diversity [13]. Extreme climate changes that occurred in the Pleistocene forced adaptation/selection of most biota, including phlebotomine sandflies, to resulting biotic and abiotic conditions [94]. Sandflies may have been subjected to local cycles of dry and cold periods, taking refuge and adapting to the more permanent humid resting fragments, which allowed these haplogroups and diversity to evolve [13,95]. Indeed, plant species' distribution shifts resulting from climate variation during the Pleistocene, have also been associated with changes in Lu. longipalpis complex diversification [18,95,96]. Arrivillaga et al., [18,19] suggest that divergence was probably a result of vicariance events that occurred throughout the late Pliocene and Pleistocene (e.g. Andean orogeny). However, recent analyses indicate a highly conserved geographic coverage of the ecological niche of Lu. longipalpis from the Last Glacial Maximum (LGM) in northern Mesoamerica (88.1%), the species complex´northern limit [13]. Hence, genetic diversity is highest in the most temporally conserved landscapes. This hypothesis can now be analyzed between and within the major Lu longipalpis haplogroup clades.
This is the first report on the genetic diversity of Lu. longipalpis from Argentinian populations, which has high genetic differentiation, two genetic clusters, and three haplogroups. Phylogeographic results provide evidence for a high level of divergence among the eight haplogroups identified for the Lu. longipalpis complex using ND4. Finally, the findings represent only the first stage of future studies required to include a more balanced sampling across Lu. longipalpis distributions and a greater number of samples not only within Argentina, but in all continental subregions. Multilocus genetic analyses will also be required in order to more completely understand evolutionary processes in this important vector species complex, and the impact of environmental change on vector transmission risk of VL. The understanding of the biological and evolutionary aspects of this species complex at the micro and macro-evolutionary levels are central to understanding the interplay between vector capacity, vectorial competence for different Leishmania parasites, urban colonization potential (micro-environment adaptation), demographic aspects of Leishmania transmission, and the clinical expression of disease [30,55]. and Marcio Antúnez for their help in field and laboratory work. We thank Mundo Sano for the samples from Clorinda.