One piece of the puzzle: Modeling vector presence and environment reveals seasonality, distribution, and prevalence of sandflies and Leishmania in an expansion area

The recent geographic spread of Leishmania infantum along the borders of Argentina, Brazil and Paraguay has been highlighted. In our previous study, Lutzomyia longipalpis was found in 55 of 123 patches surveyed, and in some patches, sandflies were found at higher densities, forming hotspots. Based on the One Health approach, we investigated the seasonality of the vector, the presence of parasite DNA, and the environmental factors that contribute to vector and parasite dispersal in these previously described hotspots in Foz do Iguaçu, Brazil. Entomological surveys were conducted monthly for one year. Fourteen hotspots peridomicile and six intradomicile were sampled. PCR was used to assess the prevalence of Leishmania DNA in sandflies. Zero-inflated negative binomial regression was used to determine the association of micro- and mesoscale environmental variables with the occurrence and abundance of the three most abundant sandfly species sampled. A total of 3543 species were captured, with Lutzomyia longipalpis being the predominant species (71.78%) of the 13 species found. Evandromyia edwardsi, Expapillata firmatoi, Micropygomyia ferreirana and Pintomyia christenseni were reported for the first time in the region. NDVI, distance to water, precipitation, west-to-east wind, wind speed, maximum and minimum relative humidity, and sex were significant variables associated with vector presence/abundance in the environment. Vector presence/abundance in the peridomicile was associated with precipitation, altitude, maximum temperature, minimum and maximum relative humidity, west-to-east wind, wind speed, and sex. Leishmania DNA was detected in an average of 21% of Lu. longipalpis throughout the year. Vector abundance is concentrated in urban and peri-urban areas, with some specimens present in different parts of the city and some sites with high vector abundance. This distribution suggests that the risk of actual contact between humans and parasite vectors in urban areas during the epidemic period is associated with patches of peri-urban vegetation and then extends into urban areas.


Introduction
Leishmaniasis is a neglected disease that spreads slowly and adapts to new environments and regions of the world. During the 20th century, infection usually occurred in forests or rural areas and was transmitted by the bite of infected female sandflies. However, environmental changes, including agricultural encroachment, deforestation, road and dam construction, irrigation, mining, urban expansion, and, more recently, pipeline construction, have altered the behavior and distribution of vectors and parasites. Endemic transmission has changed, and many infectious disease outbreaks have occurred in new areas [1][2][3][4]. In the last four decades, rapid urbanization of VL has been observed, combined with unplanned urbanization, expanding into adjacent rural areas where zoonotic cycles of leishmaniasis occur, promoting its rapid spread and human infection [5][6][7][8][9]. A major environmental change in Brazil linked to visceral leishmaniasis (VL) epidemic was the construction of gas pipelines in the 1990s [10], which significantly affected the geographic spread of leishmaniasis [2,11,12]. Since the state of Paraná (previously free of leishmaniasis) borders the states of Mato Grosso and São Paulo (with leishmaniasis presence), it was expected that VL would be entering through one of these borders. But in 2011, health authorities received a report [13] indicating the presence of Lu. longipalpis on the Argentine side bordering the state. In 2012, this species was detected in the western region of Paraná State [14] leading to the proposal of a project for the tri-border area (Argentina, Brazil and Paraguay) to evaluate vectors, reservoirs and the possible presence of human cases [13]. On the Brazilian side, the first step was a cross-sectional study to investigate the presence of vectors, Leishmania infantum infected dogs and the spatial distribution of the pathogens. We showed that different species of sandflies are present in many areas of the city of Foz do Iguaçu, where Lu. longipalpis was present in 45% of patches surveyed with Leishmania DNA at high prevalence [15]. In Foz do Iguaçu, 24% of the dogs were infected with L. infantum in 54% of the 123 patches examined [15]. Based on this alarming scenario, it is essential to understand how the abundance of vectors and their infection with L. infantum fluctuates throughout the year and what environmental factors influence it. Moreover, adopting an approach based on One Health helps us to identify transmission risk scenarios that can support specific prevention and control strategies. This study is part of the International Development Research Centre (IDRC) research project #107577-2, which aims to study leishmaniasis epidemics in the triborder region.

Study area, collection and identification of sandflies
This entomological study was conducted in Foz do Iguaçu, (25 • 32 ′ 52"S, 54 • 35 ′ 17"W), Brazil. The sandflies were collected from November 2014 to October 2015. Hotspots present in four different units, from 14 in the peridomicile area and six in the intradomicile area, with the highest sandfly abundance previously reported [15] were selected for collection (Fig. 1) [16,17] were set up 1.5 m above the ground for three consecutive nights in each month, from 05:30 p.m. to 07:30 a.m. After collection and screening, all sandflies were separated by sex [15] and identified morphologically according to the Galati identification keys [18]. After identification the female sandflies, including engorged and nulliparous, were placed in 2-mL tubes containing 70% ethanol and maintained at − 20 • C [15].

Detection of Leishmania DNA in wild-caught female sandflies
For molecular assessment of Leishmania presence, we worked with a pool of a maximum of three individuals of the same species. DNA was extracted [15] and amplification of ribosomal internal transcribed spacer 1 (ITS1) was performed using LITSR and L5.8S primers [19]. As an internal control, primers were used to amplify the IVS6 region of the cacophony gene [20]. Positive control (MHOM/FR/78/LEM75 strain) and negative control (ultrapure water) were added to each PCR reaction. Amplification products were separated on 1.5% agarose gel and visualized after staining with ethidium bromide. Positive electrophoretic PCR products were purified and sent to Macrogen (Korea) for sequencing. The electropherograms were checked using BioEdit and compared with sequences deposited in GenBank [21] using Blastn. The prevalence of L. infantum DNA in the population was determined as [22] with a significance level of p < 0.05. Statistical analysis was performed with GraphPad software.

Environmental abiotic variables
The climate was classified corresponding to a temperate oceanic or subtropical highland climate (Cfb) [23]. To assess the environmental variables influencing the presence and abundance of sandfly species, fourteen micro-and mesoscale environmental abiotic variables selected for this study are shown in Fig. 2. The maximum (max) and minimum (min) temperatures (T) and the relative humidity (RH) were recorded during the sampling period with digital thermo-hygrometers (TFA, Germany) in each domicile unit. Rainfall and wind velocity data were obtained from a weather station in the city of Foz do Iguaçu. The Normalized Difference Vegetation Index (NDVI) was used to highlight the presence of vegetation and plant biomass in each studied area [24]. Three LandSat 8 satellite images from February, August, and October 2014 were used to generate the NDVI. The NDVI values ranged from − 1.0 to +1.0, values closer to − 1.0 representing the absence of vegetation. Maps were made based on bands 6 (B6 -near infrared band) and 5 (B5 -near red band) of the LandSat 8 satellite and processed using the open software QGIS, version 2.18 [25]. The normalized difference water index (NDWI) allowed us to estimate the imperviousness of the areas by assessing degree of anthropization. The NDWI was generated based on the mid-infrared (6 TM and 5 OLI) and near-infrared bands (5 -OLI and 4 -TM) [26]. To test which environmental variables predicted the presence and abundance of sandfly species found at each hotspot, we applied zero inflated negative binomial (ZINB) regression using the pscl package [27] in the R environment [28]. The ZINB model considers the excessive number of zeros and the overdispersion found in the data by combining a logistic model to test for the presence and absence of individuals, and a count component that assumes that the predicted abundance values are drawn from a negative binomial distribution [29,30]. As there was a high degree of multicollinearity between the variables (Fig. 2), we used bidirectional stepwise regression with backward elimination to select the predictors to be included in the model [31]. The non-significant variables (p > 0.05) for each model were excluded from the analysis. Then, we created a correlation matrix between all pairs of environmental factors and, among the significant factors, the variables with the highest p-value of each collinear pair (p > 0.7 or p < − 0.7) were also excluded. Subsequently, we regressed all remaining candidate variables against the presence and abundance of specimens using the ZINB model. To avoid over-parameterization, variables that did not increase the predictability of the model were subsequently removed. Therefore, we tested the elimination of variables with the highest p-value that were not significant (p > 0.05) in each component of the model. This approach led to three different possibilities: 1) removal of the variable from the logistic component of the model, 2) removal from the count component of the model, and 3) simultaneous removal of variables from both components of the model. The models were compared using the Akaike Information Criterion (AIC), and the model with the lowest AIC value was selected. The elimination process was repeated until the elimination of candidate variables did not decrease the AIC value for the tested models. Analyses were performed separately for three species (Lu. longipalpis, Nyssomyia whitmani, and Ny. neivai), and in the intradomicile and peridomicile traps.

Sandfly fauna and seasonality
During the 12 months period, the total sampling effort amounted to 10,122 h, and 3554 sandflies were collected and classified into 13 different species (Table 1). The most abundant species was Lu. longipalpis (71%), followed by Ny. whitmani (19.62%) and Ny. neivai (1.92%). The present study reports, for the first time, Evandromyia edwardsi, Expapillata firmatoi, Micropygomyia ferreirana, and Pintomyia christenseni in Foz do Iguaçu. Unidentifiable specimens were named Phlebotominae sp. (Table 1). The average male to female ratio for all species together was 3.47:1.00, but for Lu. longipalpis it was 7.31:1.00. For the three most abundant species (Lu. longipalpis, Ny. whitmani, and Ny. neivai), peridomicile and intradomicile sites were analyzed for male and female behavior (Fig. 3). Lu. longipalpis was present over 10 months (excluding June and July) and occurred in all 20 intra-and peridomicile hotspots. The highest abundance in intra-and peridomicile occurred in March and April. Both sexes were more abundant between January and April. Lu. longipalpis enters the houses from February to June. In the entomological survey previously mentioned, the abundance of this species was 55.7% on the Brazilian side [15], and 74.2 and 47.9% in Argentina and Paraguay, respectively [32,33]. In different parts of Brazil, the abundance of Lu. longipalpis, vary from 25 to 97%, depending on the region [34][35][36][37], and males generally outnumber females [35,36,38,39].
In Foz do Iguaçu, Lu. longipalpis occurred mainly in urbanized areas. The same was observed on the Argentinean side, but the period of highest abundance of Lu. longipalpis was in early fall [9], just after the rainy months. In Latin America, several studies have reported a higher prevalence of Lu. longipalpis associated with the rainy season, especially in northeastern Brazil [34,[40][41][42][43][44]. In an outbreak of VL in the city of Belo Horizonte, Minas Gerais (central Brazil), this vector showed the highest abundance from October to March, increasing progressively until February, coinciding with the rainy months. Then, the population began to decline in April, reaching its lowest level from June to August, the driest months. In the Brazilian state of Mato Grosso do Sul, a survey conducted for two consecutive years showed the occurrence of several peaks of Lu. longipalpis: the first in February and the second in April, with a higher frequency (72%) in the rainy season compared to the dry season (28%) [45]. The climate in Paraná State is subtropical and the seasons are not as distinct as in other regions of the world or northern Brazil, which justifies the presence of Lu. longipalpis in all seasons. Ny. whitmani accounted for 19.62% of the sandfly fauna, whereas Ny. neivai accounted for only 1.92%, with the females of these species predominating in colder months. In the cross-sectional study [15], the abundance values were 28.3 and 11.6 for Ny. whitmani and Ny. neivai respectively. In Argentina and Paraguay, the abundance of Ny. whitmani was 25% (Puerto Iguazú-Argentina) and 38.8% (Alto Paraná Department-Paraguay), respectively [46,47]. Nyssomyia whitmani has an increase in the male abundance in August, and January in the peridomicile. However, in the intradomicile, the number of individuals sampled increased in April, May and August. Males and females entered the houses and were prevalent in the peridomicile between April and October, decreasing significantly in July. Nyssomyia neivai was more abundant in the peridomicile than intradomicile. The highest abundance in the peridomicile environment was observed in unit C, followed by units A and D. Among the 20 hotspots, five of them (329, 264, 448, 470, and 616) had the highest abundances. In the intradomicile, hotspots 329 and 470 had the highest numbers of captured females, especially in the autumn and winter. In our study, Ny. whitmani females were found in the intradomicile mainly in fall and winter. In the State of Rio de Janeiro, Ny. whitmani was abundant during the coolest months (June, July, and August), although both occurred throughout the year [48]. In Argentina, Ny. whitmani is mainly present in less urbanized areas, with peaks of abundance in early spring and summer [46].

Environmental variables
From the 14 environmental variables studied, those that were significant in at least one of the ZINB models (count or logistic) were summarized ( Table 2 and Fig. 3). In the intradomicile, no variable was associated with the presence of these species, and the following variables were significantly correlated with sandfly abundance: distance from water (+, positively), NDVI (− , negatively), west-to-east wind (+), wind speed (− ), and sex (+) for Lu. longipalpis; maximum (+) and minimum moist (− ) for Ny. whimani; and rain (− ) for Ny. neivai.
For the peridomicile, altitude (− ), maximum temperature (− ), and minimum (+) and maximum (− ) moist predict the presence of Lu.  longipalpis, whereas the altitude (+), west-to-east wind (+) and sex (+) were associated with its abundance. Rain (− ) and wind speed (− ) was correlated with the abundance of Ny. whitmani, and the wind speed (− ) with its presence in peridomicile. Finally, wind speed (+) predicts the peridomicile abundance of Ny. neivai. The important goal was to know which variables were significant for the presence of Lu. longipalpis in the intra-and peridomicile. The significant variables for intradomicile were distance from water, NDVI, WE wind and wind speed like, as proposed by Salomón [7]. In the peridomicile habitat, altitude, WE wind, maximum temperature, and minimum and maximum relative humidity were significant. Other studies have supported that the distribution of Lu. longipalpis is associated with environmental factors such as temperature, humidity, and wind [48,49]. In the present study, Ny. whitmani intradomicile abundance was significantly associated with minimum and maximum relative humidity, while only wind speed was related to both abundance and presence in the peridomicile.

Leishmania infection rate in females
Among the 712 females of the main species, 361 pools with three females were tested for DNA presence of Leishmania spp. In some pools only two females are present, if less than three females were found in a trap or in a single night. Of these, 136 pools belonged to Lu. longipalpis, 126 pools to Ny. whitmani, 77 pools to Phlebotominae sp. (non-identified), 17 pools to Ny. neivai, two pools to Pintomyia pessoai and Nyssomyia intermedia, and one pool to Migonemyia migonei. The mean positivity for Leishmania DNA was 12.7% throughout the year and most of the females were identified as Lu. longipalpis (56 positive pools), followed by Ny. whitmani (15 positives pools) and Ny. neivai (3 positives pools). Among the unidentified females, Leishmania spp. DNA was detected in 10 pools. The highest number of females with Leishmania DNA was found in the trap-329 in both peri-and intradomicile. Table 3 shows Leishmania DNA prevalence by species, by area and by season. Fifty-six of the 85 pools with positive PCRs were sequenced, allowing species identification of the parasite in 34 of them. The remaining 22 pools showed a mixed of Leishmania spp. (L. infantum, and L. braziliensis sequences in the electropherogram, which prevented species identification through Sanger sequencing. Leishmania infantum accounted for 94% of the Leishmania DNA presence among the different phlebotomine species. Only one pool of Ny. neivai presented Leishmania braziliensis DNA. For Lu. longipalpis, the highest prevalence of positive females was in autumn (44.9%). The percentage prevalence of positive females in winter (30.8%) attracted our attention. However, it should be noted that in winter represented only 11% of the total number examined. In Argentina a prevalence of 3.9% was reported [32]. In Alto Paraná Department-Paraguay, this value was 23.4% for Lu. longipalpis [33]. In the Americas, studies conducted so far have shown that infection rates can vary depending on the country/region, season or vector species. In Brazil, the prevalence of Leishmania DNA ranged from 0.2 to 36.5%. Considering the different states that constitute the Brazilian territory, a prevalence of 0.2% of infected females was observed in Bahia, 2.6% in Mato Grosso do Sul, 2.7 to 3.9% in Minas Gerais, 3.7% in Maranhão, 36.5% in Ceará and 4.8 to 7.2% in São Paulo [50]. The situation is not different in the Old World, where Italy has a prevalence of 2.9% [51]. In Morocco, PCR with ITS as the target has shown a prevalence of 7.3% [52]. Also, in Madrid, the infection rate was 58.5% of Phlebotomus perniciosus samples for L. infantum using kDNA-PCR methods [53], while in Tunisia, the prevalence of L. infantum within P. perniciosus was 0.2% using nested ITS-PCR [54]. It would certainly be interesting to use the same methodology to compare infection rates. Therefore, it is important to know the infection rates and the periods in which they occur.
Additionally, Nyssomyia whitmani and Ny. neivai have been implicated in the transmission of L. braziliensis under sympatric conditions [15]. Thus, we show that the pathogenic complex (Leishmania/hosts) can consist of several vectors and genotypes of the parasite and a main reservoir: the dog. Another aspect that has attracted attention is that Table 2 Zero these two elements (vector and reservoir) are closely related to humans. In this context, the theoretical (genesis of the focus) and practical (combat strategy) interests are fundamental to propose control measures. In this study, the period with the lowest prevalence of DNA in sandflies was spring. However, females with Leishmania DNA were present in 10 of the 12 months studied. Chemical control can be a way to control vector populations when the first population appears (August, September), especially in the hotspot areas identified in this study. In winter, chemical control can be used to reduce the population in human residences and animal shelters. Vulnerable populations need support and education to understand the life cycle of Leishmania and sandflies and to avoid transmission and aggravation of leishmaniasis.
The "macrofocus" of L. infantum can be maintained and rapidly spread below the 54S parallel in areas where conditions are now favorable and climate change has already been observed through temperature increases of about 2 • C [55]. The scenario observed for the triborder area indicates that vector control measures must involve all three countries. In addition, the tri-border area shows a scenario that favors the flow of people and goods between Brazil, Paraguay and Argentina, which to this day is characterized by an intense flow of commuters. Another aspect to consider is the presence of a hydroelectric power plant in the tri-border region. The environmental changes caused by the discharge of the Itaipu dam (related to increased water inflow into the dam, and consequent increase in lake surface area and evaporation) should be considered when proposing vector control measures during these periods. It is important to note that extreme weather events and accelerated climate change may impose limitations on the model used. Therefore, eco-epidemiological surveillance of these vectors in hotspot areas is an important strategy to reduce the transmission of leishmaniasis in all its forms. Because of the strong link between leishmaniasis and environmental factors, One Health approaches are key to leishmaniasis control. For example, educational practices cannot be based solely on providing information about the disease and risk factors, but require a more comprehensive intervention. The perspective of aligning surveillance and health education activities with health promotion is essential for the success of community-based interventions. Therefore, it is important to broaden the discussion beyond the scientific environment or health services and to involve all professional categories and population groups that can be active in community change (e.g., NGOs, professional associations, community leaders). It is also important to note that when pipelines or hydroelectric dams that alter the environment are planned (a feature of the region studied), the health sector of each country must be informed of the risks involved. The multidisciplinary nature of diseases such as leishmaniasis must be respected, and teams from all disciplines, not only parasitologists and entomologists, but also anthropologists, sociologists and urban planners, must be involved in the risk analysis.

Conclusions
The results presented in our work show that the new sandfly colonization is concentrated in urban and peri-urban hotspots, with some individuals present in different parts of Foz do Iguaçu, Brazil. However, few sites have high vector abundance. Our data and modeling suggest that the risk of actual contact between humans and parasite vectors in urban areas during epidemic periods is associated with peri-urban vegetation patches, with subsequent dispersal into urban areas; the period of highest transmission of L. infantum in the study area was January to May. If measures are taken to reduce access to sources of infection, there will be fewer sandfly infections and consequently fewer human infections. The implementation of vector control measures on the brazilian side of the border from September to December will reduce the sandfly population and the risk of transmission of protozoa to humans. Finally, policy makers are encouraged to use the applied methodology to implement strategies to control and reduce Leishmania vector populations. Finally, we emphasize that the fight against cutaneous and visceral leishmaniasis cannot be limited to vector control, but must be based on the One Health concept.

Table 3
Number of females assessed for the presence of Leishmania sp.. DNA, and its prevalence for each sandfly species and area surveyed in Foz do Iguaçu, Paraná. The prevalence of Leishmania sp. in the different sandfly species was calculated using the average number of individuals in the pool. The minimum, average (bold), and maximum infection rates are displayed. Species with negative results are not included.