Integrating environmental, entomological, animal, and human data to model the Leishmania infantum transmission risk in a newly endemic area in Northern Italy.

INTRODUCTION: Historically, leishmaniasis in Italy was constrained to areas with Mediterranean climate. In the last 20 years, sand fly vectors (Phlebotomus perniciosus), cases of canine leishmaniasis (CanL) and cases of human visceral leishmaniasis (VL) have been observed in Northern Italian regions, traditionally classified as cold areas unsuitable for sand fly survival. AIM: We aim to evaluate through a One-Health approach the risk of endemic transmission of Leishmania infantum in the Piedmont Region, Northern Italy.
METHODS: We collected environmental, entomological, animal, and human data. We applied a geostatistical binomial model to map the probability of P. perniciosus presence in the study area, using selected environmental parameters as predictors. We evaluated the spatial relationship between the probability of P. perniciosus presence and the geographical distribution of CanL and VL cases observed between 1999 and 2013.
RESULTS: Between 1999 and 2003, 142 sampling sets (17%) out of 839 resulted positive for P. perniciosus. Elevation, degree of slope, normalized difference vegetation index (NDVI) and summer temperatures were associated with positive sampling sets. During the study period, 164 (13.6%) of Piedmont municipalities reported at least one autochthonous case of CanL, while 89 VL cases were observed in 54 municipalities (4.5%). We observed an association between municipalities affected by autochthonous CanL cases and the estimated probability of P. perniciosus presence (Odds Ratio for 10% increase of probability: 2.66; 95% confidence intervals (CI): 2.16-3.37). We found that human VL incident cases were positively associated with the probability of the municipality of residence of being endemic for CanL (Incidence Rate Ratio for 10% increase of probability: 1.49; 95% CI 1.02-2.16).
CONCLUSIONS: Using a One-Health approach, we quantified the spatial association between the distribution of P. perniciosus, municipalities endemic for CanL and incident cases of human VL, suggesting that the disease has become endemic in the Piedmont region.

Introduction

Human visceral leishmaniasis (VL) is a vector-borne disease caused by a parasite of the Leishmania genus [1]. In the Mediterranean basin, Leishmania infantum represents the main putative agent of VL [2]. The parasites are transmitted to susceptible organisms, such as humans and dogs, by the bite of female phlebotomine sand flies [1]. Human VL is usually diagnosed clinically among subjects that present irregular fever, anaemia, leukopenia and hepatosplenomegaly. Whereas most immunocompetent individuals will not develop any symptoms after Leishmania infection, immunocompromised subjects develop a VL, which is typically fatal if untreated [3]. In the canine host, L. infantum causes the canine leishmaniasis (CanL), a severe disease that occurs in fewer than 50% of infected animals. Both symptomatic and asymptomatic infected dogs also act as reservoir hosts for the human disease [4].

Transmission of L. infantum is associated with the geographical distribution and abundance of competent vectors, mainly represented in Italy by Phlebotomus perniciosus [5,6]. The geographical distribution of phlebotomine sand flies is affected by several environmental parameters, such as altitude, land-cover, vegetation, air temperature and relative humidity [7,8]. For instance, air temperature is believed to influence both the vector distribution and the seasonal activity of the sand flies [9,10]. Consistently, cases of both canine and human leishmaniasis have been historically reported in central and southern peninsular Italian regions, characterized by hot summers and mild winter temperatures [11,12]. However, in the last 20 years, both sand fly vectors and cases of CanL and human VL have been observed in Northern Italian regions, traditionally classified as colder areas characterised by continental climate [13].

Assessing and measuring the endemic status of human VL in recently affected areas, such as Northern Italy, raises difficulties. Due to the long VL incubation period and the possibility of frequent travels to endemic areas (such as trips to the Mediterranean coastal areas of Northern Italy), the occurrence of human cases in non-endemic regions could be an unreliable indicator of new endemicity. On the contrary, assessing the geographical distribution of the vectors and canine cases does not directly allow evaluation of the transmission risk to human populations. For these reasons, a One-Health approach, integrating environmental, entomological, animal, and human data has been proposed as a valid option to study and control the transmission of Leishmaniasis [14,15]. With this study, we have evaluated, through a One-Health approach, the risk of transmission of L. infantum in the Piedmont region, an area of northern Italy with an increasing number of human VL cases [13,16].

Methods

Study area

Piedmont is a region in the north-west of Italy with a surface of 25,387 km2. The altitudinal range varies from 4632 m a.s.l. of the highest peak to 72 m a.s.l. of the Padan Plain. Piedmont surface topography is characterized by 43% mountains, 30% hills and 26% plains. According to the Köppen and Geiger classification, in plains and hilly areas the climate is continental temperate or subcontinental temperate with hot summers and generally cold winters [17]. However, the area presents different climatic features because of the presence of several microclimates determined by the presence of the Alps mountains, hills, the Padan Plains and Lake Maggiore. The air temperatures show a regular decrease with elevation, except for changes due to geographical setting [18].

The region is divided in 1206 municipalities and 8 provinces and has roughly 4.3 million inhabitants. The population density is on average 174 inhabitants per km2, but is highly heterogeneous, reflecting the orography of the area. Population density steadily increases from low densely inhabited mountainous areas, to hilly areas, and then to the densely inhabited Padan Plains, reaching its peak in the urbanized area of the Metropolitan City of Turin.

Data collection

Entomological data

We obtained information about the presence/absence of P. perniciosus from 839 sampling sets randomly identified from potentially suitable areas across Piedmont Region between 1999 and 2003. In total 208 out of 1206 municipalities (23%) of Piedmont region were covered by the sampling activities. Briefly, each sampling set consisted in 10 to 20 sticky traps (total area: from 0.4 to 0.8 m2) embedded with castor oil placed in suitable sand fly diurnal resting sites (animal shelters, houses, and scarp wall cracks) for one night, from sunset to sunrise [7]. Sampling activities were performed during the summer months (from the second half of June to the first half of August) when the highest sand fly densities are expected [9]. Collected insects were stored in 70% ethanol until further analysis. All specimens were cleared in chloral-lactophenol and mounted in Hoyer medium and morphologically identified according to taxonomical keys [19,20].

Environmental data

We selected four environmental predictors potentially associated with the presence of P. perniciosus from the literature [5,7,9,10,21], namely elevation, slope, normalized difference vegetation index (NDVI) and summer maximum temperatures. Data about elevation were obtained from the Digitalized Elevation Model (DEM) raster with a resolution of 1 × 1 km released by the Italian Institute of Environmental Protection and Research (ISPRA) [22]. Slope was computed from the DEM raster in GIS environment and was expressed in degrees to describe the angle of inclination of the terrain. The NDVI annual average values in the period 2001–2005 were computed using data from the Moderate Resolution Imaging Spectroradiometer (MODIS) Terra satellite with 1 × 1 km resolution and used as a proxy of greenness. Monthly average temperatures for the period 2001–2017 were obtained from MODIS Land Surface Temperature (LST) with an approximate resolution of 1 × 1 km. The monthly daytime mean temperatures for June, July and August were averaged and used as a proxy for the daytime summer temperatures.

Canine data

For each municipality of the Piedmont region, we obtained information about the presence or absence of detected autochthonous cases of CanL between 1998 and 2013. Data were retrieved from an online map produced within a study conducted by the National Institute of Health. In this study, municipalities were classified in relation to the canine leishmaniasis as “disease free” or “affected” [23]. Presence of endemic CanL was assessed by testing with the Indirect Fluorescent Antibody Test (IFAT) serum samples of dogs that were likely to be infected or sick and that had not travelled to endemic areas.

Human data

The Hospital Discharge Register (HDR) records were used to identify incident cases of VL among residents of the Piedmont region, using the International Classification of Diseases, 9th Revision, Clinical Modification (ICD 9 CM) codes. Hospital discharges for VL (ICD 9 CM = 085.0 in any diagnostic position) recorded between January 1, 1995 and December 31, 2013 (19 years) were selected. In this database, each hospitalized patient is identified by an anonymous ID code. Patients with VL can be readmitted to the hospital multiple times, since disease relapses are frequent. We used the ID code to identify, for each patient, only the first VL-related hospitalization and to exclude readmissions. Since VL can lasts many months, we introduced a lag time of 4 years to define incident cases. By introducing a 4-year lag, we identified incident cases for the period 1999–2013 (15 years) excluding subjects that had been already admitted for VL between 1995 and 1998. In addition, we restricted our study to Italian citizens. In the study area, indeed, VL is a rare disease and almost all identified VL cases with foreign nationality were natives of countries where VL is endemic. To perform spatial analysis, cases were attributed to their residential address at the time of hospitalization. Population data for the Piedmont region were obtained from the Italian National Statistical Institute (ISTAT) [24].

Statistical analysis

Geostatistical modelling of the presence of P. pernicious

An initial exploratory analysis was carried out to assess the relationship between the environmental predictors and the proportion of traps positive for P. perniciosus. Through a multivariable logistic regression, we modelled the association between selected environmental predictors and the proportion of traps positive for P. perniciosus. Based on the results of an exploratory analysis, the effect of NDVI was modelled with a linear function, while the effect of elevation and summer temperatures were modelled through quadratic splines. Slope was dichotomised into flat (slope < 2°) and steep areas (slope ≥ 2°). Scatter plots and equations included in the final model are reported in the supplementary materials. In the final multivariable analyses, geostatistical binomial regression with a logit-link function was used to model the proportion of positive sampling sets. In order to lower the computational burden of the geostatistical models, sampling sites located less than 250 m from each other were attributed to the same sampling site. More specifically, the proportion of positive sampling sets (N = 839), p(x), at site xi (N = 749) was modelled as:where: d(x) is the vector of the georeferenced explanatory variables; S(x) is a random effect that accounts for the spatial correlation between observations induced by unmeasured factors affecting the presence of P. perniciosus. In this case, S(x) was modelled as a zero-mean stationary and isotropic Gaussian process with spatial correlation function given by: Cor {S(x), S(x’)} = exp. (−||x − x′||/ϕ). Where ϕ is the scale parameters which regulates the rate of decay of the spatial correlation for increasing distance, and ||x − x′|| is the distance in space between the two spatial points x and x′. From the fitted geostatistical model, we plotted a predictive map of the probability of the presence of P. perniciosus over a 1 × 1 km grid. The gridded probability was aggregated at the municipality level to evaluate its relationship with the presence of VL and CanL cases. Given that the inhabitants are not distributed homogeneously within a municipality area, we computed the population-weighted average of the estimated gridded probability, using a population density raster with 1 × 1 km resolution [25]. Parameter estimation and spatial prediction were carried out using the PrevMap package in the R software environment [26]. Maps presented were created using Quantum GIS 3.4.

Statistical modelling of human VL cases, endemic CanL and P. perniciosus presence

For each municipality, we studied the association between the estimated probability of P. perniciosus presence and the endemic status for CanL. We fitted a logistic regression model with CanL (0: municipality free from CanL; 1: municipality with autochthonous cases of CanL) as the outcome and p(Phl) (predicted probability of P. perniciosus presence at the municipality level) as the explanatory variable:

From the fitted model we estimate p(CanL), namely the probability of each municipality of being endemic for CanL given the P. perniciosus presence.

To model human cases, we used a negative binomial regression. We modelled age- (0–25, 25–75, >75 years) and sex-standardized incidence rates of VL cases at the municipality level, computing the expected incident cases for each stratum. We fitted a negative binomial regression to incident cases at municipality level, using p(CanL) as the explanatory variables:where μ and E are respectively the observed and the expected number of human VL cases at municipality level, while p(CanL) is the estimated probability for each municipality of being endemic for CanL estimated from Model 1.

Results

Between June 1999 and September 2003, 839 sampling sets in the study area were tested (Fig. 1). Overall, 142 sets (17%) resulted positive for P. perniciosus. All environmental predictors resulted associated with the probability of positive traps sets. The relationships between the proportion of positive sets and elevation, slope, NDVI and summer temperatures are shown in Fig. 2. At the municipality level, the median value of the estimated probability of the presence of P. perniciosus was 11.3% (5th–95th centiles: 0.3%–30.5%). The geographical distribution of the estimated probability is shown in Fig. 3. The highest probability was observed in the central hilly area of the region and in the pre-Alps, while the lowest probability was predicted in the mountain areas and in the eastern Padan Plains. Between 1998 and 2013, 164 municipalities of the Piedmont region (13.6%) reported at least one autochthonous case of CanL (Fig. 4). Between 1999 and 2013, 89 human VL cases were observed in 54 municipalities (4.5%) of the study area (Fig. 5). At the municipality level, we observed a positive association between the presence of autochthonous CanL cases and the estimated probability of P. perniciosus presence (OR for 10% increase of p(Phl): 2.66; 95% CI: 2.16–3.37). When we fitted the negative binomial regression to VL cases, we found that the incident cases were positively associated with the probability of the municipality of residence of being endemic for CanL (IRR for 10% increase in p(CanL): 1.49; 95% CI 1.02–2.16) (Table 1).

Fig. 1

Study area and distribution of sampling sets.

Left Panel: the framed area includes the study area. Right panel: P. perniciosus sampling sites. White dots: negative sampling sets, red dots: positive sampling sets. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2

Proportion of positive sampling sets and continuous environmental predictors (altitude, NDVI, temperature).

Upper Panel: Fig. A: Logit transformation of proportion of positive sampling sets in relation to elevation, size of circle is proportional to the number of observations in the range shown in the x axis. Fig. B: Fitted value of Logit transformation of proportion of positive sampling sets and 95% CIs.

Central Panel: Fig. A: Logit transformation of proportion of positive sampling sets in relation to NDVI, size of circle is proportional to the number of observations in the range shown in the x axis. Fig. B: Fitted value of Logit transformation of proportion of sampling sets and 95% CIs.

Lower Panel: Fig. A: Logit transformation of proportion of positive sampling sets in relation to summer temperatures (average 2000–2017), size of circle is proportional to the number of observations in the range shown in the x axis. Fig. B: Fitted value of Logit transformation of proportion of sampling sets and 95% CIs.

Fig. 3

Geographical distribution of the estimated probability of the presence of P. perniciosus aggregated at the municipality level.

Fig. 4

Geographical distribution of the municipalities affected by autochthonous CanL cases.

Fig. 5

Geographical distribution of human VL incident cases by municipality.

All cases recorded in a given municipality have been geographically attributed to the centroid of the municipality of residence. Dot size is proportional to the log incidence rate at the municipality level.

Table 1

Association between the probability of the presence of P. perniciosus, the presence of canine leishmaniasis and human visceral leishmaniasis in Piedmont municipalities.

Model 1: logistic regression
Outcome	Exposure	OR	95% CI
CanLm	p(Phlm)	2.66	2.16–3.37

OR: Odd Ratios, CI: Confidence Intervals, IRR: Incidence Rate Ratio, CanL: Municipalities affected by autochthonous CanL cases; μ: Incident VL cases at the municipality level; p(Phl): estimated probability of P.perniciosus presence aggregated at the municipality level (expressed as 10% increase); p(CanL): probability of a municipality of being endemic for CanL (expressed as 10% increase).

Discussion

The Piedmont Region, north-west Italy, was not previously considered endemic for L. infantum, but recent studies have documented new foci of infection [13,27,28]. With the current study, we used a One-Health approach, integrating environmental, entomological, animal and human data, to systematically assess and map the risk of endemic transmission of L. infantum in that region, covering a 15 years timespan.

Results provided by the geostatistical prediction of P. perniciosus presence showed a heterogeneous distribution of the vector in the study area. The vector distribution was limited to specific areas, characterized by suitable environmental features. Altitude is one of the most important factors influencing the distribution of P. perniciosus [5]. In the present study, the probability of observing P. perniciosus resulted higher among hilly and premountainous areas (elevation range 250–750 m) and lower in high mountainous areas (elevation > 1000 m) and in flat areas of the Padan Plains, in line with previous findings [5]. In addition, we observed that the proportion of positive traps sets for P. perniciosus was higher in areas characterized by warmer summer temperatures and higher NDVI, which is a proxy of vegetation cover. Summer temperatures greatly influence the ability of the insect to fulfil the lifecycle during the activity season and the L. infantum development within the vector. The highest probability of the presence for P. perniciosus was found in areas with temperature ranges between 26 and 32 °C, confirming previous experimental and field data [5,29]. The positive association between P. perniciosus presence and vegetation cover, described by NDVI in the model, could be explained by the higher presence of organic material that ultimately represents a suitable site for winter larval survival and development. Results provided by the geostatistical model are also consistent with the geographical distribution of canine and human cases. We observed a positive association between municipalities affected by autochthonous CanL cases and the estimated probability of P. perniciosus presence. These results were further strengthened by the positive association observed between the incident VL cases and the geographical distribution of CanL. Taken together, these results show that the transmission of Leishmaniasis has become endemic in Piedmont in recent years.

The main limitation of the current study is that the entomological samples are based on monitoring activities performed between 1999 and 2003, while data about CanL cases and VL cases are referred to a timespan that runs until 2013. It cannot be excluded that the geographical distribution of P. perniciosus might have slightly changed after sampling, however it is likely that this change has not altered sensibly the geographical heterogeneity of P. perniciosus presence in the area. Saying it another way, areas that presented a higher presence of P. perniciosus during the sampling period are likely to have kept a higher probability also in the following period.

The main strength of our study lies in the integration of several sources of information, including environmental, entomological, animal and human data, using a One-Health perspective to evaluate the risk of endemic transmission of L. infantum in the Piedmont region. We found positive associations between the spatial distribution of the vector, infected reservoir hosts and human cases. These findings indicate that L. infantum strains have been circulating in some areas of the region during the study period. Our results also suggest that it is likely that a proportion of the 89 VL cases observed in Piedmont between 1999 and 2013 have contracted the infection in Piedmont, in contrast to the hypothesis that the majority of VL cases reported in Piedmont region are imported from neighbouring areas, such as coastal areas, historically considered endemic for L. infantum. The observation of the occurrence of human VL cases alone would not have been sufficient to identify areas at high risk of endemic L. infantum transmission, given the low occurrence of human cases in the study area. On the other hand, our study suggests that, since the majority of VL cases occur in areas where autochthonous cases of CanL are observed and infected dogs are easier to identify (up to 50% of infected dogs develop symptoms [30]), monitoring canine cases of leishmaniasis together with the sampling of sand fly may be a useful tool for public health authorities to identify high-risk areas.

The emerging transmission of L. infantum in new affected areas could be explained by the introduction of infected dogs from endemic areas, the expansion of the vector distribution, or the combination of the two [2]. It has been argued that the geographical distribution of Leishmaniasis and its vectors in Europe has been changing over the last decades [8,31,32] and, in particular, that the recent northward spread might be a consequence of the increase in temperatures values and extreme climate indices driven by climate change [33]. Climate change might affect leishmaniasis distribution by the effect of air temperatures on the survival and abundance of its vector species. For instance, recent studies have predicted that climate change could trigger the spread of sand fly vectors in areas outside the present geographical range and lead to an increased risk of contracting leishmaniasis for the human population in several areas of the world [6,34,35]. Our study cannot directly elucidate whether the endemization process observed in the study area is attributable to climate change. However, the expansion of the vector distribution induced by climate change remains one plausible mechanism to explain the northwards shift of L. infantum observed in Italy. Cold temperatures play a fundamental role in limiting the geographical distribution of the parasite and its vectors [6,34,35] and continental winters might have prevented the Piedmont region from the endemization of Phlebotomus species until the end of the last century. Piedmont, however, in the last decades experienced a non-negligible increasing trend in temperatures, with an increase of 2.3 °C in the average maximum temperatures and an increase of 1.5 °C in the average minimum temperatures in the last 60 years [36]. Observed increasing temperatures could have triggered the expansion of sand fly vectors starting from the end of the last century in Piedmont. To verify this hypothesis, repeated sand fly sampling in old and new sites in the coming years may clarify if and how the geographical distribution and abundance of this vector is increasing and may help to identify the main drivers of such a change.

In conclusion, the current study has shown the use of a One-health approach to identify the geographical distribution of endemic transmission of L. infantum in recently affected areas. The results support the hypothesis that an endemization process for L. infantum has occurred over the last two decades in the Piedmont region and can be helpful in addressing public health interventions to control disease transmission in emerging high-risk areas, where both human and canine cases occur.

Funding

None.

Declaration of Competing Interest

None.