Leishmaniasis is a parasitic disease caused by several different species of protozoan parasites (Holakouie-Naieni et al., 2017), and there is no definitive treatment for this disease (Nilforoushzadeh et al., 2007). Two different hosts are required to complete the life cycle of the Leishmania parasite, where humans and some other mammals, especially canines (that also act as reservoirs), are the definitive hosts with sand flies belonging to the genus Phlebotomus in Old World and Lutzomyia in New World acting as intermediate hosts (Bates, 2007; Akhoundi et al., 2016; Bari and Rahman, 2016).
The least fatal form of the disease, cutaneous leishmaniasis (CL), is caused by Leishmania major and is the most common clinical form of disease in the Middle East. However, the disease affects an additional 60+ countries in tropical and subtropical regions (Jacobson, 2011). The global incidence is 0.7-1.2 million new cases per year worldwide (Alvar et al., 2012) and around 1.7 billion people are estimated to live in areas where leishmaniasis is common (Pigott et al., 2014). CL is a serious public health problem in many rural areas of Iran (Akhavan et al., 2010). A total number of 589,913 cases of CL was reported between 1998 and 2013, with the annual incidence of 30.9 per 100,000 in the Iranian population (Holakouie-Naieni et al., 2017). Isfahan Province in Iran is an endemic region for CL and the main vector with respect to humans is the Phlebotomus papatasi sand fly with the gerbil Rhombomys opimus being the main nonhuman reservoir (Yaghoobi-Ershadi et al., 2005; Yaghoobi-Ershadi, 2012). Distribution and abundance of vectors and reservoirs of this disease is affected by different climatic, socioeconomic and cultural factors (Desjeux, 2001; Mollalo et al., 2014). Moreover, unplanned urbanization and environmental changes such as irrigation, dam construction and desertification increase the risk of infection (WHO, 2008; Nilforoushzadeh et al., 2014).
Recent advances in geographic information systems (GIS) and remote sensing (RS) have promoted the study of spatial epidemiology and environmental factors affecting the vector-borne diseases (Kassem et al., 2012). GIS can determine resource allocation and is a valuable approach in implementation of control measures (Barbosa et al., 2014). A number of epidemiological studies around the world have used GIS-based methods in risk mapping studies and identifying endemic area of diseases (Sudhakar et al., 2006; Bhunia et al., 2010; Salahi-Moghaddam et al., 2010; Mollallo et al., 2014, 2015; Rajabi et al., 2014; Hagenlocher and Castro, 2015). These techniques have been used for developing landscape predictors of sand fly abundance as an indicator of human vector contact. Bhunia et al. (2010) studied the effects of several topographic indexes on the prevalence of visceral leishmaniasis (VL) using GIS and RS. They investigated the relation between altitude, temperature, humidity, rainfall and the normalised difference vegetation index (NDVI) with the prevalence of the disease in the north-eastern Indian sub-continent, while Salahi- Moghaddam et al. (2010) studied VL in an endemic area of Iran using the GIS approach and Mollalo et al. (2014) focused on the relationship between the vegetation cover and occurrence of CL in Golestan Province in Iran based on satellite images.
The aims of the present study were to 1) identify the geographical distribution of CL in Isfahan Province; 2) search for hotspot areas; and 3) assess the relations between the climatic and topographic factors with CL incidence in Isfahan Province using the spatial analysis during the period 2011 to 2015.
Materials and Methods
Isfahan Province is located between 31° 43′ to 34° 22′ N and 49° 38′ to 55° 31′ E and lies in the central parts of the Iranian plateau covering an area of 107,027 km2. Iran is a mountainous country mainly situated ≥1000 m above the mean sea level. The provincial capital is the historic city of Isfahan. According to the census of 2015, it consists of 23 districts with about 4,600,000 inhabitants. The province has a moderate and dry climate on the whole, and is a wellknown endemic area of leishmaniasis (Figure 1).
Data collection and preparation
In order to conduct this study, information on the total number of CL cases, population at risk, climatic data, vegetation coverage and altitude were gathered for each district of the province.
Climatic data including mean precipitation, mean temperature and mean humidity were obtained from Tehran Meteorological Center, collected from synoptic stations in Isfahan and neighbouring provinces, including Lorestan, Kohgiloyeh VA Boyerahmad, Semnan, Chaharmahal & Bakhtiyari, Yazd, Fars, Qom and Markazi. These data were entered in ArcGIS, version 10.3 (ESRI Inc, Redlands, CA, USA) and dealt with in several steps: first, the yearly average of climatic data was calculated for 34 stations from April 2011 to March 2016; a point layer was created for 34 synoptic stations in a second step followed by ordinary Kriging was carried out for interpolation and calculation for all districts to derive a predicted value for unmeasured locations. Weights were based on the distance between the measured points, the prediction locations, and the overall spatial arrangement among the measured points (http://support.esri.com). The output for this step is a raster layer that was subsequently overlaid on the Isfahan polygon layer at the district level. Finally, the spatial analyst extension of ArcGIS, was used for calculating the average of each climatic variable. This function summarizes the values of a raster within the zones of another dataset (zonal statistics) and reports the results as a table (ArcMap 10.3, Spatial Analyst; ESRI Inc., Redlands, CA, USA).
Information on vegetation status of study area was obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) consisting of 16-day composites with 250-meter spatial resolution. These data were used to calculate the NDVI, one of the most commonly used measures to of landscape ecology and useful for the study of the epidemiology of vector-borne disease (Bavia et al., 2005). The NDVI is related to the proportion of photo-synthetically absorbed radiation and is calculated as follows:
where NIR stands for near infrared light and RED for red light (Jackson and Huete, 1991; Meijerink et al., 1994). The NDVI data for Isfahan Province covering April 2011 to March 2016 was downloaded from the United States Geological Survey (USGS) website (https://lpdaac.usgs.gov).
Digital elevation model data
A digital elevation model (DEM) is an ordered array of numbers that represents the elevation over a specified segment of the landscape (Meijerink et al., 1994). The DEM data for our study area was obtained from the Shuttle Radar Topography Mission (SRTM) through the USGS website (https://lta.cr.usgs.gov/SRTM) and overlaid on Isfahan polygon with district boundaries. Finally, zonal statistics was used for computing average of altitude for each district of Isfahan.
Population data for all 23 districts of Isfahan was obtained from National Bureau of Statistics of Iran.
Cutaneous leishmaniasis data
We obtained confirmed CL cases for Isfahan Province from the Isfahan University of medical sciences and Department of Communicable Disease Control (CDC) of the Iranian Ministry of Health in Tehran. The CL incidence rate was calculated as follows:
The total population at risk was defined as population for each district or province, which was obtained for each year.
A choropleth map was produced to show distribution of CL cases at the district level using a range of colours. Spatial data often are clustered, which means that stronger relationships may be present between proximate observations (Fotheringham and Brunsdon, 1999). Therefore, to explore spatial heterogeneity of disease distribution we used the Global Moran’s Index in ArcGIS, version 10.3 to map the clustering of CL cases across districts in Isfahan. Moran’s Index is a commonly used indicator of spatial autocorrelation (i.e. the correlation of a single variable between pairs of neighbouring observations) as well as non-random accumulation that indicates clustering. This index ranges from −1 to +1, where the value 1 means perfect positive spatial autocorrelation (high values or low values cluster together), and −1 perfect negative spatial autocorrelation (a checkerboard pattern), while 0 shows full spatial randomness. Z-scores and P values are calculated along with Moran’s Index to evaluate its significance (Moran, 1948).
Moran’s Index is a global statistic and does not show the local structure of spatial autocorrelation (Chen, 2013). Therefore, we used the Getis-Ord (Gi*) statistics (1995) to determine high-risk areas, i.e. hotspots of the disease. This statistic is suitable for distinguishing between cluster structures of high or low strength. In our case, a high value would represent a cluster of high incidence values of CL (a hotspot), while a low value would represent a lowlevel cluster (a coldspot). The Gi* statistic is usually standardized based on its sample mean and variance, which means that a Zscore can be calculated to indicate the degree of statistical significance. Positive and negative Gi* statistics with high absolute values show events created by clusters of high and low values, respectively. Conversely, a zero Gi* value implies random distribution of the spatial events observed.
To assess the correlation between DEM, NDVI, humidity, temperature, precipitation and maximum wind speed with CL incidence, we applied overlay analysis. At first, we produced three layers of information including a geographical map the study area (Isfahan), six layers representing each of the environmental factors interpolated with the Kriging method, and a layer indicating the CL distribution. The different layers were overlaid in various combinations, including a final step of all layers together.
Our study showed that incidence of CL in the Isfahan at the district level was significantly clustered (Moran’s Index=0.17, z=4.02, P<0.001). A total of 13,790 confirmed CL cases were reported in the Isfahan Province during 2011-2015 period, with a maximum incidence rate of 62 per 100,000 in 2014 (Figure 2). The choropleth map of CL incidence in the 5-year time frame is shown in Figure 3, which demonstrates that Ardestan and Natanz had the highest incidence.
The results of the Getis-Ord (Gi*) statistic, used to identify hotspots areas of CL during 2011-2015, showed that both hotspots and coldspots existed at the district level with the former seen in Ardestan and Aran va Bidgol (P<0.01) and Naein and Natanz (P<0.05), while Chadgan was identified as a coldspot area with a CL incidence of 0.05<P<0.1 (Figure 4).
Figure 5 shows the overlay of DEM on the map of CL incidence. A high incidence of CL was found at altitudes between 600 and 1,800 m. Overlaying of the NDVI and CL incidence layers showed a high incidence of CL in areas with a low level of herb coverage (0.073-0.110) (Figure 6). We also found a high incidence of CL in areas with relative humidity of 27-30% (Figure 7), mean temperature of 15-19°c (Figure 8), mean precipitation of 5-20 mm (Figure 9) and maximum wind speed of 12-16 m/s (Figure 10).
The present study investigated spatial patterns of CL incidence in an endemic area of Iran using the GIS and RS, during 2011-2015. The assessment of spatial characteristics of the CL cases by Moran’s Index and derived Z-scores indicated that CL cases, as a whole, were clustered in the study area. We identified hotspots as well as coldspots for CL incidence, which were clustered in a specific area. The fact that the highest and lowest CL incidences were found in 2011 and 2015, respectively, could be due to interventional programmes, such as reinforced training, health education, disease surveillance and strengthened vector/reservoir control interventions, which were performed these years.
As can be seen in Figure 5, CL incidence was predominantly found at altitudes between 600 and 1,800 m, which represents moderately low-laying parts of the province and, indeed, the country as a whole, which includes large areas of 2,000 m above the sea mean level. We identified several high-risk areas of CL in middle of Isfahan, mainly in Natanz and Ardestan. Other studies have also demonstrated the high incidence of CL in Ardestan (Nilforoushzadeh et al., 2014). The finding that CL was more prevalent in areas with moderately low and low altitudes is in accordance with a large number of other studies reporting an inverse relationship between altitude and CL incidence (Guernaoui et al., 2006; Bhunia et al., 2010; Holakouie-Naieni et al., 2017). Guernaoui et al. (2006) collected 2,742 specimens belonging to nine phlebotomine species and reported that P. papatasi, the vector of L. major, is more common in the lowlands. Bhunia et al. (2010) showed the significant effect of altitude on spread and outbreaks of VL and found a high prevalence at low altitudes, while most of the highlands had only few cases. A possible explanation could be that active transmission of CL increases when the population density increases as it does with decreasing altitude (Cohen and Small, 1998). Moreover, the higher temperature, lower humidity and the character of vegetation cover in areas with lower altitude may contribute to a higher rate of CL due to the positive effect of these variables on vector populations and reservoir hosts (Salahi- Moghaddam Abdoreza et al., 2010; Ali-Akbarpour et al., 2012; Mollalo et al., 2015).
NDVI is commonly used to separate three types of land cover: surfaces with sparse vegetation (NDVI<0.2), surfaces partially covered by vegetation (0.2≤NDVI≤0.5) and surfaces fully covered by vegetation (NDVI>0.5) (Momeni and Saradjian, 2007). In addition, it has been shown that an extensive vegetation cover provides shade which reduces the surface temperature, while the air temperature decreases through the process of evapotranspiration (Chaithanya et al., 2017); therefore, low-level vegetation coverage is almost always accompanied by higher temperatures and evaporation as well as less rainfall and lower relative humidity, a situation which provides favourable conditions for sand flies (Mollalo et al., 2014; Shirzadi et al., 2015). Accordingly, our study found a high incidence of CL cases in areas with low vegetation (NDVI<0.165). This finding is in line with reported findings by Bavia et al. (2005) in Brazil, by Bhunia et al. (2010) in India and by Gadisa et al. (2015) in Ethiopia.
We also found hotspot areas in semi-arid regions with moderate levels of humidity (Figure 7). These areas have also been shown as important foci of CL in Ethiopia (Gadisa et al., 2015), in India (Sudhakar et al., 2006), in Brazil (Barbosa et al., 2014) and several investigators in Iran have reported similar results (Ali- Akbarpour et al., 2012; Shirzadi et al., 2015). The contribution of humidity is, however, difficult to separate from the effect of rainfall because of the interaction of these two factors. It has been shown that peaks in rainfall can lead to reductions in sand fly numbers since excess precipitation reduces the amount of suitable resting sites for adult sand flies, limits their flight activity, and interferes with reproduction by sweeping away the eggs (Medlock et al., 2014).
Annual variation between moderately high wind speed and CL distribution showed a positive relationship in our study. A recent study has shown that P. similis, the potential vector of L. tropica in Greece, is more common in areas with low (5.9–7.7 m/s) and medium (7.7–8.6 m/s) mean wind speeds compared to very low air movement (<5.9 m/s). However, above wind speeds 8.6 m/s the risk decreased substantially (OR = 0.8; 95% CI=0.6–1.0; P= 0.002) (Ntais et al., 2013). An investigation in Darab District located in Fars Province in southern Iran showed that a wind speed above 3 m/s had a preventive effect on P. papatasi, the dominant sand fly in that area (Askari et al., 2017). Contradictory results from different researches may be explained by the dual role of wind speed in disease distribution. For example, although sand fly biting opportunities is mitigated by strong winds, the flight distance can increase (Wu et al., 2016). Moreover, the interaction between climatic factors in the epidemiology of vector-borne disease should be considered, since some climatic factors such as temperature, humidity and wind speed always operate together and interact with each other in nature (Seid et al., 2014).
There are some limitations for this study. The CL surveillance system in Iran is a passive system, so underreporting is a strong possibility, especially in the rural areas. Furthermore, some errors may occur in the surveillance system, such as unreliable diagnosis and notification, or cases acquired in areas other than where they were diagnosed and reported. In addition, climate is only one of many groups of factors influencing vector distribution, while other factors such as vector ecology and socio-economic factors vary from one area to the other and should also be considered in the study of vector ecology. However, we assessed only climatic factors, while we fully understand that comprehensive research needs to consider also other factors, such as cultural, socioeconomic, immigration, demographic, sanitation and vector diversity.
CL is a public health problem in Isfahan. Several hotspot areas were identified using spatial analysis performed by GIS and RS. Overlay analysis revealed a relationship between several climatic factors and incidence of CL in these hotspot-prone areas, the majority of which were located in semi-arid regions with low vegetation coverage. We also found fewer hotspots in lower-altitude regions with higher temperatures and less rain. In addition, a positive correlation between wind speed and hotspot areas was found. The results of the present study indicate that GIS is a feasible approach for identifying spatial disease patterns and detecting hotspots of particular infectious diseases.