Spatial distribution and cluster analysis of a leishmaniasis outbreak in the south-western Madrid region , Spain , September 2009 to April 2013

D Gomez-Barroso1,2,3, Z Herrador (zherrador@isciii.es)3,4,5, J V San Martín5,6, A Gherasim4, M Aguado6, A Romero-Maté6, L Molina6, P Aparicio4,5, A Benito4,5 1. Network Biomedical Research Centre in Epidemiology and Public Health (CIBERESP in Spanish), Madrid, Spain 2. Health Institute Carlos III (ISCIII), Madrid, Spain 3. These authors contributed equally to this article 4. National Centre for Tropical Medicine, ISCIII, Madrid, Spain 5. Network Biomedical Research on Tropical Diseases (RICET in Spanish), Madrid, Spain 6. Fuenlabrada University Hospital, Madrid, Spain

Since July 2009, there has been a community outbreak of leishmaniasis in south-west Madrid, Spain.The present study used the spatial distribution of cases to investigate the connection between the outbreak and a recently built peri-urban park.We included 157 cases of cutaneous (CL) and 90 cases of visceral (VL) leishmaniasis diagnosed at Fuenlabrada University Hospital between July 2009 and April 2013.CL and VL cases were geo-referenced and incidence rates by census tract were calculated.To identify high-risk areas, the spatial autocorrelation between individual cases was estimated.In a next step, areas where risk of disease was significantly increased were identified by cluster analysis.Higher incidence rates and the areas with highest intensity of CL and VL were located in the north-western part of the municipality.The most likely cluster of CL comprised three census tracks with relative risk (RR) = 11.5 (95% confidence interval (CI): 9.2-13.6).Two additional significant VL clusters were detected, the most likely one with RR = 9.2 (95% CI: 7.3-11.1).In addition, we found one significant VL cluster in the immigrant population (RR = 12.8; 95% CI: 9.3-16.1).The spatial pattern of leishmaniasis transmission revealed a relation between the outbreak and the suspected risk area.

Background
Leishmaniasis is a parasitic disease caused by more than 20 protozoan species of the genus Leishmania.It is transmitted through the bite of female sandflies of the genera Phlebotomus (Old World) and Lutzomyia (New World).An ecological system in which a Leishmania species is maintained indefinitely is usually formed by one principal reservoir host, e.g.dogs for Leishmania infantum in both the Old and New Worlds [1].Recent reports indicate that in Spain, other animals such as wild carnivores, rabbits and hares may play a role in the maintenance of the system, occasionally bringing the parasite from its enzootic focus into contact with humans [2,3].
In Spain, leishmaniasis is endemic and both the visceral (VL) and the cutaneous (CL) forms are caused by L. infantum.Phlebotomine sandflies belonging to the subgenus Larroussius serve as vectors, while the dog is the main reservoir [4].The real prevalence of leishmaniasis in Spain is unknown.Notification of leishmaniasis is not mandatory at national level (mandatory only in 12 of 17 regions) and there is no national leishmaniasis control programme [5].Moreover, Leishmanias sp.cases are underascertained and underreported in Spain [6].
In the autonomous region of Madrid, notification of leishmaniasis has been mandatory since 1997 and is the most common zoonotic disease affecting infants and in patients with human immunodeficiency virus (HIV) infection [7,8].Between 2000 and 2009, the annual incidence rate was around 0.5 per 100,000 population (between 12 and 25 leishmaniasis cases per year) [9].During the last quarter of 2010, the case number increased fivefold.Subsequent research confirmed that an outbreak of leishmaniasis had been ongoing since July 2009 in the south-western area of the region of Madrid, with most cases occurring in the city of Fuenlabrada [5,10].
Most of the L. infantum isolates obtained from patients in the outbreak area presented an uncommon strain: the combined genotype L-920 [2].The initial hypothesis to explain this outbreak postulated risk areas for bites by Phlebotomus sandflies infected with this new strain, creating a temporospatial clustering of infected individuals.Because leishmaniasis was found in only 7% of the dogs examined in the outbreak area, other reservoirs were suspected and investigated.Around 30% of the hares studied were positive for Leishmania parasites [3,5], suggestive of a sylvatic transmission cycle possibly linked to a recently built park called Bosquesur, adjacent to the urban area of Fuenlabrada [11].
Spatial pattern analysis has been found to be useful to better understand disease transmission of parasitic diseases when there is a strong correlation between the spatial distribution of the disease and its hosts [12,13].We aimed to assess the spatial distribution of CL and VL cases and the cluster occurrence within the city of Fuenlabrada (Madrid) as well as the distribution of the cases in relation with the park Bosquesur.

Study area
Fuenlabrada is a city and municipality located in the Madrid metropolitan area in central Spain.It is the fourth biggest town in the community of Madrid and is located in the south-west, ca 22 km from the capital.The climate is Mediterranean, with an estimated average annual rainfall of 450-500 mm.The yearly mean maximum and minimum temperatures are 20.2 °C and 7.6 °C, respectively.The mean altitude of the area is 664 m above sea level and the municipality is divided into 108 census track units.Figure 1 shows the municipality's land cover (based on data from the European Environment Agency [14]) and population density at census track level to describe the study area.

Cases
CL and VL data from September 2009 to April 2013 were supplied by the Internal Medicine and Dermatology Departments of Fuenlabrada University Hospital.Diagnosis was confirmed by direct observation of the parasite in a skin biopsy specimen or by a positive polymerase chain reaction (PCR) as well as by isolation in Novy-MacNeal-Nicolle (NNN) culture.Further details on laboratory methods are provided elsewhere [2].Data contained information on a number of clinical, epidemiological and demographic variables; those considered in this study were sex, age, address of residence, migrant status and date of symptom onset.Frequencies and percentages were used to summarise data.The differences were assessed by Student's t-test and chi-square test for continuous and categorical variables, respectively.

Population
The population data was obtained from the 2011 annual statistic published by the City of Fuenlabrada, stratified by sex and age, at census track level [15].Fuenlabrada had an area of 39.1 km 2 and an estimated population of 204,100 in 2011, of whom 75.3% were in the age group of 16-64 years and 15.4% were immigrants.

Spatial analysis
Geo-referencing of all VL and CL cases was carried out using Google Earth to obtain the geographical

Figure 1
Land use and population density, Fuenlabrada (Madrid), Spain, September 2009-April 2013 Maps were drafted using the free-licence source: Spanish National Institute of Geography.National Plan for Aerial Orthography (PNOA project).Available from: http://www.ign.es/ign/layoutIn/actividadesFotoTelePNOA.do[15].coordinates of patients' residences.In order to study the features of the territory, we obtained the map of Fuenlabrada by census track from the National Institute of Statistics (2011) [16] and downloaded images from the National Plan for Aerial Orthography (PNOA) project [17].

Spatial distribution
We calculated the incidence rates for CL and VL by census track adjusted by four age groups (0-14, 15-44, 45-64, ≥ 65 years) and by sex.We plotted the map with rates to understand the distribution of the disease at census track level.We calculated the local Moran's index in order to study the local indicator of spatial autocorrelation (LISA) between the rates.This index assesses local associations by comparing local averages to global average [18].Its significance is estimated by generating a reference distribution using 999 random permutations.The LISA significance map includes the following categories: 'high-high' indicates clustering of high value rates (positive spatial autocorrelation), 'low-high' indicates that low value rates are adjacent to high value rates (negative spatial autocorrelation), 'low-low' indicates clustering of low value rates (positive spatial autocorrelation), 'high-low' indicates that high values are adjacent to low value rates (negative spatial autocorrelation), and 'not significant indicates that there is no spatial autocorrelation.
In order to understand the risk distribution, we estimated the spatial smoothing distribution of CL and VL by means of the kernel density function.This function is based on the quadratic kernel function described by Silverman [19].This method is an interpolation and smoothing tool used to generalise the position of a point to an area.Kernel density estimation fits a curved surface over each case such that the surface is highest above the case and zero at a specified distance (the bandwidth) from the case.
To estimate the distance to the park Bosquesur, which was the main suspected risk area according to previous investigations [3,11,20], we drew a polygon around this park.Distances to other landscape elements were also assessed.We created buffers with different distances: 500, 1,000 and 2,000 m around the park Bosquesur to measure the distances of the cases´ addresses.Then, we calculated the average nearest neighbour index (ANNI) to calculate the distance between the location of each case and their nearest neighbour.This method was used as first approach to test whether or not the cases were clustering.If the average distance was below the average for a hypothetical random distribution, the distribution of the analysed characteristics was considered to be clustered.The index is expressed as the ratio of the observed distance divided by the expected distance.Thus, if the index is < 1, the pattern exhibits clustering, while index of > 1 indicates a trend towards dispersion.A significance level of 99% was chosen in the analysis.

Cluster analysis
Spatial clusters were analysed using the SaTScan spatial statistic estimator developed by Kulldorff [21].To assess CL and VL spatial clusters in the entire population, and separately in the immigrant population, we Maps were drafted using the free-licence source: Spanish National Institute of Geography.National Plan for Aerial Orthography (PNOA project).Available from: http://www.ign.es/ign/layoutIn/actividadesFotoTelePNOA.do[15].

A.Cutaneous leishmaniasis B. Visceral leishmaniasis
Bosquesur park LISA Not Significant High-High High-Low Low-High Low-Low used the scan statistic estimator to perform a purely spatial analysis, based on the assumption of a Poisson distribution.This method consists of creating a circular window which scans the entire study area.In our study we restricted the spatial window to a maximum radius of the average distance between cases (250 m).
The radius of the centroid varies continuously in size from zero to the specified upper limit, in our case 250 m.The circle with maximum likelihood and containing more cases than expected is denominated the most likely cluster.
An increase in observed cases above the number expected was assessed using Monte Carlo test simulations (999 replications) with a 95% confidence interval.
For the spatial analysis, we used Arcgis version 10.0, free software GeoDa and SaTScan.Data analysis was performed using SPSS version 18.0.

Spatial distribution
The incidence rates of CL and LV in the different census tracks varied between zero and 1,003 and between zero and 613 cases per 100,000 population, respectively.Figure 3A shows the distribution of CL rates.Higher rates of CL (between 500 and 1,100 cases/100,000 population) were observed in three census track in the north of the municipality, while the highest rate for VL was found in a different northern census track (Figure 3B).
The LISA for CL and VL rates are shown in Figure 4. Significant clusters of high values for CL and VL were detected: a CL hotspot (p < 0.005) in a census track located in the north-west of the municipality, and two significant VL hotspots (p < 0.005) in two census tracks with high values adjacent to census tracks with low value rates.
The average distance among the cases (bandwidth) was 250 m.Following spatial smoothing, areas with high intensity for CL were identified in the north of the city, close to the park Bosquesur (Figure 5A).The distribution of the VL cases is shown in Figure 5B; high intensity areas were located in the northern part of the municipality.
Table 1 summarises the cases distribution and distance to the park.We observed a decreasing trend in CL and VL occurrence in relation to the distance to the park.Up to 75% of CL cases and 70% of VL cases lived in places less than 1 km from this suspected risk area.
The ANNI was 0.505 (p < 0.001) for CL and 0.582 (p < 0.001) for VL throughout the study area.Both for CL and VL, the pattern exhibited clustering.The observed average distances were 98 and 270 m and the expected average distances were 193 and 463 m for CL and VL, respectively.

Cluster analysis
Four significant CL clusters were identified through spatial cluster analysis (p < 0.001, Figure 6A).The most likely cluster comprised three census tracks in which 24 cases were diagnosed during the study period, while the number of expected cases was 2.46 (relative risk (RR) = 11.50 and p < 0.005).Other secondary clusters were located in the north of the municipality.
The VL spatial cluster analysis detected two significant clusters (p < 0.001, Figure 6B; Table 2).The most likely cluster comprised one census track with eight cases and an expected number of cases of 0.95.The RR was 9.15 (p < 0.005).The other secondary cluster was located in the north-eastern part of the municipality.
When we analysed the data according to the migrant status, we found one significant cluster for VL in the immigrant population, with four observed cases, 0.35 expected cases and a RR of 12.75 (p < 0.005) (Figure 6C).This cluster included one census track located in the north of the municipality, different from the one detected for VL in overall population.No significant cluster of CL was identified through spatial cluster analysis in the immigrant population.

Discussion
Our study revealed the spatial characteristics of human CL and VL during an outbreak in Fuenlabrada (Madrid) using geographic information system (GIS) tools and spatial statistical analysis.Similar approaches have been used in Spain to investigate the spatial distribution of the sandfly vector and canine leishmaniasis [4,22].However, to our knowledge this is the first attempt to implement spatial techniques to assess the distribution and cluster occurrence of human CL and VL cases in an outbreak in Europe.
We have described the evolution of the VL and CL cases in Fuenlabrada separately, as analysing them together may have resulted in misinterpretation.VL was more common in men than in women, while no sex differences were found for CL.This seems to be concordant with findings from previous epidemiological studies indicating that VL occurs more frequently among adult men [20,23].Although it has been hypothesised that sex hormones play a role in the modulation of immunity against leishmaniasis [24], the explanation for this trend still remains uncertain.
In our research, the census tracks with the highest incidence of CL were different from the ones with highest incidence of VL.It is difficult to determine whether the different distribution of CL and VL cases determined by disease onset, vector distribution, disease transmission pattern and/or presence of the host.
Leishmaniasis is known to be a diverse and complex disease [1], therefore further studies may be needed to investigate this particular finding.
Although the spatial distribution of CL and VL leishmaniasis did not overlap perfectly in space, both were spatially clustered along the border of a park.Risk associated with other landscape elements within the municipality were analysed without significant results.Leishmaniasis is one of the main parasitic diseases of the world for which the transmission profile includes landscape elements and environment [25].In our research, all spatial methods deployed showed that the northern peripheral census tracks were the most heavily affected.Moreover, CL and VL incidence rates and the distance to the park Bosquesur were spatially correlated.Finally, the cluster analysis also showed that the most likely CL and VL clusters were close to this park.
Our results were in accordance with previous research that emphasises that this disease can be associated with urbanisation close to vegetation areas [23,26,27].Worldwide, leishmaniasis outbreaks have been related to human activities close to or within forested areas: road construction, building of dams, irrigation schemes and horticulture development, and establishment of new residential colonies lead to intrusion into the sylvatic cycle of the disease [27,28].The park Bosquesur is a man-made natural area, located between the towns of Alcorcón, Leganés, Fuenlabrada, Getafe and Pinto.During the outbreak period, leishmaniasis cases were reported from these five towns [10].Adjacent to this newly planted area, public construction work was carried out on one of the main roads (M-407) to Fuenlabrada in 2009-10.All these recent changes at peri-urban level may have altered the ecology of this area and the transmission dynamics of leishmaniasis.
The sum of observed cases in the four significant clusters of CL was 61 (39% of the total 157 CL cases included in the analysis), while the percentage of VL cases included in the significant clusters was even lower, 14% (n=13) of the 90 VL cases.This difference in case clustering can also be observed in the epidemiological curve (Figure 2), although we cannot draw conclusions in this regard as we did not carry out a temporospatial analysis.Other exposure pathways are also plausible.As shown by de Almeida et al., the pattern of leishmaniasis is not static and the disease may spread from one area of a municipality to another [29].
On the other hand, the use of a patient's home residence as a marker of the place of infection may not be accurate because contact between host and vector can have occurred outside the home.
The park has several recreational areas with footpaths and bicycle trails [30].Social networks, recreational activities and other interactions between human settlements and the peri-urban natural environment may play a role in the distribution of the disease in this context [26,28].Nevertheless, these establishments are attended mostly during daytime.First appearance of active P. perniciosus females occurs at sunset and the density peak is usually reached around 23:00-24:00 [5,11].Therefore, the highest probability of transmission may be not associated with the presence of individuals inside the park, but rather in homes or outdoor places close to, but outside the park.
According to Arce et al., patients were asked during the outbreak investigation about recreational habits, but neither particular areas nor classic environmental risk factors were identified [10].It should be noticed that only cases were assessed during the official investigation.We believe that the outbreak investigation would have benefited from a control group in order to evaluate exposures and risk factors.
A high percentage of VL cases (39%) were immigrants and we identified a cluster of VL in this particular group.This cluster did not overlap with the VL cluster detected in overall population.We also carried out the analysis in the non-immigrant population only (data not shown) and did not find any relevant difference compared to the analysis on entire population.According to official data from the local authorities, 15.4% of the Fuenlabrada population in 2011 were immigrants, of which 14.9% are from Sub-Saharan Africa [15].The percentage of immigrants in the study area was similar to that of other urban areas in Spain [31].Several factors could have been responsible for the high percentage of VL in immigrants: migration of non-immune people from areas where L. infantum is not endemic [1,32], differences in host immune responses [33], and differences in living habits and/or health-seeking behaviour.These factors should be explored in further investigations.Our study has several limitations.One relates to the long incubation period of the disease, rendering determining and interpreting the information related to transmission site difficult.Control measures may have affected the spatial distribution of cases.We are aware of disinfection activities undertaken since 2011; although no information on the targeted sites is publically available, we cannot discard their impact on the case distribution.
We chose Kulldorff's method for the spacial cluster analysis because it has several advantages [34]: it adjusts for population density and confounding variables (e.g.age and sex); there is no pre-selection bias since the clusters are selected without prior hypothesis on their location, size or time period; the statistical test takes into account multiple testing and delivers a single p value; if a cluster is detected, its location is specified.
In a molecular typing study, Chicharro et al. found that the outbreak was not caused by a single parasite strain, as four combined genotypes were found.At Fuenlabrada Hospital, the ITS-LOMBARDI type was isolated from all serotyped cases.This L. infantum ITS type has been present in this region since at least 1992 [2].A high density of P. perniciosus has also been observed in the park [10].In previous research carried out in central Spain, a correlation was found between vector density and cases living in areas between villages or at the limits of a village [22].Future entomological and molecular typing studies could benefit from our results, as a combined methodology could allow more precise conclusions regarding the transmission patterns.

Conclusion
Although our study design did not allow establishing causal associations, the methodology can be considered useful in generating hypotheses during an outbreak investigation.Future work should examine the role of vector density, seroprevalence of Leishmania in canine and other possible reservoirs, climate variability, socio-economic conditions, land use and changes made by humans to the habitat over a longer time span in the study area.This will allow building accurate risk maps and targeting prevention and treatment interventions in these high-risk areas in a timely manner.