Climate change and potential distribution of zoonotic cutaneous leishmaniasis in Central Iran: Horizon 2030 and 2050

Babak Shiravand, Ahmad Ali Hanafi-Bojd, Abbas Ali Dehghani Tafti✉, Mohammad Reza Abai, Ali Almodarresi, Masoud Mirzaei

1Department of Health in Disasters and Emergencies, School of Public Health, Shahid Sadoughi University of Medical Sciences, Yazd, Iran

2Department of Medical Entomology & Vector Control, School of Public Health, Tehran University of Medical Sciences, Tehran, Iran

3Engineering College, GIS & RS Department Yazd Branch, Islamic Azad University, Yazd, Iran

4Department of Epidemiology & Biostatistics, School of Public Health, Shahid Sadoughi University of Medical Sciences, Yazd, Iran

Keywords:Cutaneous leishmaniasis Climate change RCP scenario Phlebotomus papatasi Rhombomys opimus

ABSTRACTObjective: To investigate and predict the effects of climate change on the potential distribution of the main vector and reservoir hosts of the disease in Yazd province in the future.Methods: Distribution data for vector and reservoir hosts of zoonotic cutaneous leishmaniasis in Yazd province were obtained from earlier studies conducted in the area. MaxEnt ecological niche modeling was used to predict environmental suitability. BCC-CSM1-1(m) model and two climate change scenarios, RCP 4.5 and RCP 8.5 were used for horizons 2030 and 2050 climate projections. Future projections were based on data of a regional climate change model.Results: With both scenarios in 2030 and 2050, the results of jackknife test indicated that the mean temperature of wettest quarter and temperature annual range had the greatest effect on the model for the vector and the reservoir hosts, respectively.Conclusions: The climate conditions are the major determinants of zoonotic cutaneous leishmaniasis incidence rate in Yazd Province. These climate conditions provide favorable habitats for ease transmission of zoonotic cutaneous leishmaniasis in this endemic area.Habitats suitability for the vector and reservoir will be expanding in the coming years compared with the current conditions, such that, in horizon 2030 & 2050, the probability of the presence of the vector and reservoir within 38 580 and 37 949 km2, respectively, from Yazd province is above 60%. Moreover, an increase is predicted in the presence of the vector in the western parts and the reservoir in the northern and central parts of the province in the future. Understanding the role of environmental and bioclimatic factors in zoonotic cutaneous leishmaniasis occurrence can provide a guide for policy-makers in the creation and implementation of more effective policies for prevention and control.

1. Introduction

Leishmaniasis is a widespread vector-borne disease in the tropical and subtropical regions of zoonotic cutaneous leishmaniasis of the world with approximately two million new human cases reported each year[1,2]. It is estimated that more than 12 million people are infected worldwide[3]. Cutaneous leishmaniasis (CL) continues tobe a major public health problem in Iran, with many endemic foci of zoonotic cutaneous leishmaniasis (ZCL) been reported in different parts of the country[4,5]. The sand fly Phlebotomus (P.) papatasi(Diptera： Psychodidae) is the principal vector of Leishmania major(Kinetoplastida： Trypanosomatidae), the causative agent of ZCL in the old world[6]. A number of gerbils (Rodentia： Gerbillidae)are known as reservoir hosts of ZCL in the old world, including Rhombomys (R.) opimus, Meriones libycus, Tatera (T.) indica and Meriones hurrianae which are more endemic in Iran[7]. The prevalence of ZCL has been reported in different provinces of Iran ranging from 1.8% to 37.9%[8,9]. Yazd province, located in the central part of Iran, is one of the major endemic foci of ZCL in the country. From 2007 to 2016, a total of 4 229 cases have been formally reported in Yazd province[10]. Leishmania major is responsible for most of the ZCL cases reported in Yazd province,and has been isolated from various rodent reservoirs (R. opimus,Meriones libycus and T. indica) and sand fly vector (P. papatasi)[11].

Vector-borne diseases including ZCL are generally sensitive to climatic conditions[12]. It seems like that the incidence of human cases strongly correlates with the presence of sufficiently high densities of vector and reservoir species[13]. It has been suggested that climate change can cause a shift in the geographical spread of insect populations in the future through the modification of climatic conditions which affect the reproduction and length of annual activity of these insect species[14]. Different climatic conditions can cause species’ range shifts, therefore, analyses of the distributional response of vector and reservoir species to climate change scenarios can provide insight into how the spatial epidemiology of ZCL may be affected by climate change[15].

Representative concentration pathway (RCPs) scenarios are a set of greenhouse gas concentration and emission pathways designed to support research on the impacts of climate and potential policymakers’ response to climate change[16,17]. RCP scenarios provide interpolated climate layers based on the 2003 3rdassessment report of the intergovernmental panel on climate change and are publicly available[18]. These scenarios assume different climate futures based on greenhouse gases emission pathways and their associated radiative forcing to the year 2100[18]. In total, a set of four RCPs have been developed and are named according to the change in radiative forcing by the end of the century： 8.5, 6, 4.5 and 2.6 W/m2.

Ecological niche models (ENMs) based on point occurrence data,digitized environmental layers, and machine learning algorithms, which are usually all overlaid on a geographic information system (GIS) platform,provide a better understanding of the distribution of vector-borne diseases[19,20]. ENMs attempt to predict the fundamental niche of a species, defined as the set of environmental conditions in which it can maintain populations without the need for immigration[21].The fundamental niche provides the potential distributional pattern of a species when projected to a geographical space[22].Effects of geographical constraints on dispersal and negative ecological interactions such as competition may prevent a species from occupying its entire fundamental niche[23]. Thus, taking such factors into consideration defines the actual geographical distribution of a species. Practically, ENMs combine both the ecological requirements and spatial locations of species and predict the occurrence of species in an area between the potential and actual distributions[24]. Modification of these distributions using the constraining factors mentioned above generates the actual distributions. ENMs can thus provide useful framework for testing hypotheses concerning the role of different environmental variables in determining species’ distributional patterns[25]. Estimating the potential future distributions of vector and reservoir species can help to identify potential risk areas for ZCL. The aim of this study was to investigate the effect of climate change on the potential distribution of ZCL vector and reservoir species in Central Iran in horizon 2030 and 2050.

2. Material and methods

2.1. Vector and reservoir data

Existing records from all previous studies conducted on P. papatasi,the main vector of ZCL, in Yazd province were collected and entered into a database. Previous studies in the area have reported R. opimus as the main reservoir host of ZCL in the province. Coordinates of collection sites for this gerbil were also recorded in the database.

2.2. General circulation models

The results of the research conducted by National Climate Research Institute[26] were the basis for selecting General Circulation Models in this study. BCC-CSM1-1(m) model at a spatial resolution of 30 s(1 km2) was used in our analysis. The Beijing Climate Center Climate System Model version 1.1 (BCC_CSM1.1) was used to develop the BCC_CSM1.1 (m) model with a moderate resolution[27]. In the present study, two scenarios were used for modeling： RCP4.5 and RCP8.5. In the RCP4.5 emission scenario, CO2concentration is estimated to be 650 PPM by 2100, with a radiative forcing level of 4.5 watts per square meter[28,29]. Global Change Assessment Model team at the Joint Global Change Research Institute (a branch of the Pacific Northwest National Laboratory) in the United States developed the RCP 4.5 scenario[30].The RCP 4.5 scenario is considered as a stabilization scenario, without an overshoot, in which total radiative forcing is stabilized shortly after 2100[31]. RCP 8.5 scenario corresponds to the highest greenhouse gas emission trajectory compared to the scenario literature[32], and hence also to the upper bound of the RCPs. The greenhouse gas concentrations and emission trajectories in RCP 8.5 is estimated to increase considerably over time, leading to a radiative forcing level of 8.5 W/m2at the end of the century with a temperature range between 3.5-4.5 ℃[33]. The bioclimatic data based on RCP 4.5 and RCP 8.5 scenarios for the years 2030 and 2050 were downloaded from www.ccafs.cgiar.org (http：//ccafs-climate.org/data) and www.worldclim.org (http：//www.worldclim.org/cmip5_30s), respectively,at resolution 30 s.

2.3. Modeling potential occurrence of sand flies and gerbils

MaxEnt software Ver. 3.3.3 was used to predict the most appropriate ecological niches for the target species[34]. The contribution of the environmental variables was calculated by Jackknife analysis.Variables with no contribution were omitted from the final analysis.Eighty percent of the vector/reservoir collection points were used for model training and the remaining 20% were used as Test-set.MaxEnt was used for random point’s selection. ESRI ArcGIS 10.3 software was used for the preparation, management and editing of spatial and climatic data, modeling, and presentation of the model results. The management of climatic data and preparation of the expressions for modeling were assisted by Microsoft Excel 2010 software.

3. Results

3.1. Modeling for ZCL vector and reservoir hosts in Horizon 2030: RCP 4.5 scenario, BCC-CSM1-1(m) model

The estimates of relative contributions of the environmental variables to the MaxEnt model for the vector and reservoir hosts of ZCL are shown in Table 1. The variables with the highest contributions to the models for P. papatasi and R. opimus were bio8(73.5%) and bio7 (50.4%), respectively. Area under curve (AUC)for P. papatasi was calculated as 0.917. The results of jackknife test of variable importance for P. papatasi show that the environmental variable with highest gain when used in isolation is bio8, which therefore appears to have the most useful information by itself. Bio8 is the environmental variable that decreases the gain the most when it is omitted, which therefore, appears to have the most information that is not present in the other variables (Figure 1). The AUC for R.opimus was calculated as 0.851. Jackknife test of variable importance for R. opimus revealed that the environmental variable with highest gain when used in isolation is bio7, which therefore appears to have the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted is bio7, which therefore appears to have the most information that isn’t present in the other variables (Figure 2).

3.2. Modeling for ZCL vector and reservoir hosts, Horizon 2030: RCP 8.5 scenario, BCC-CSM1-1(m) model

As shown in Table 1, the variables with the highest contributions to the models for P. papatasi and R. opimus were bio8 (38.4%) and bio7(93.8%), respectively. The AUC for P. papatasi was calculated as 0.889 (Figure 3). The results of jackknife test of variable importance for P. papatasi show that the environmental variable with highest gain when used in isolation is bio8, which therefore appears to have the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted is bio8, which therefore appears to have the most information that isn’t present in the other variables (Figure 3). The AUC for R. opimus was calculated as 0.794 (Figure 4). Jackknife test of variable importance for R.opimus revealed that the environmental variable with highest gain when used in isolation is bio7, which therefore appears to have the most useful information by itself. The environmental variable les that decreases the gain the most when it is omitted is bio7, which therefore appears to have the most information that isn’t present in the other variables (Figure 4). Environmentally suitable areas for the vector and reservoir hosts in horizon 2030 using RCP 4.5 and RCP 8.5 scenarios have been presented in Figure 5.

Table 1. Contribution of variables used for the prediction of the potential distribution of Phlebotomus papatasi and Rhombomys opimus as vector and reservoir hosts of zoonotic cutaneous leishmaniasis in Yazd Province, Central Iran, horizon 2030： RCP 4.5 and 8.5 scenarios, BCC-CSM1-1(m) model (%).

Figure 1. Receiver operating characteristic curve and output of MaxEnt model and jackknife analysis for Phlebotomus papatasi in Yazd province, Central Iran, Horizon 2030： RCP 4.5 scenario, BCC-CSM1-1(m) model.

Figure 2. The receiver operating characteristic curve and output of MaxEnt model and jackknife analysis for Rhombomys opimus in Yazd province, Central Iran, Horizon 2030： RCP 4.5 scenario, BCC-CSM1-1(m) model.

Figure 3. The receiver operating characteristic curve and output of MaxEnt model and jackknife analysis for Phlebotomus papatasi in Yazd province, Central Iran, Horizon 2030： RCP 8.5 scenario, BCC-CSM1-1(m) model.

Figure 4. The receiver operating characteristic curve and output of MaxEnt model and jackknife analysis for Rhombomys opimus in Yazd province, Central Iran, Horizon 2030： RCP 8.5 scenario, BCC-CSM1-1(m) model.

Figure 5. A： Predicted environmental suitability for vector (left) and reservoir (right) of ZCL in Yazd province of Iran using RCP 4.5 scenario, BCC-CSM1-1(m) model： Horizon 2030. B： Predicted environmental suitability for vector (left) and reservoir (right) of ZCL in Yazd province of Iran using RCP 8.5 scenario, BCC-CSM1-1(m) model： Horizon 2030.

3.3. Modeling for ZCL vector and reservoir hosts, Horizon 2050: RCP 4.5 scenario, BCC-CSM1-1(m) model, P.papatasi

Table 2 shows estimates of relative contributions of the environmental variables to the MaxEnt model for the vector and reservoir hosts of ZCL. As shown below, the variables with the highest contributions to the models for P. papatasi and R. opimus were bio8 (61.1%) and bio7 (50.4%), respectively. The AUC for P. papatasi was calculated as 0.903 (Figure 6). Analysis of Jackknife test of variable importance for P. papatasi shows that the environmental variable with highest gain when used in isolation is bio8, which therefore appears to have the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted is bio8, which therefore appears to have the most information that isn’t present in the other variables (Figure 6).

The AUC for R. opimus was calculated as 0.849 (Figure 7).Jackknife test of variable importance for R. opimus revealed that the environmental variable with highest gain when used in isolation is bio7, which therefore appears to have the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted is bio2 , which therefore appears to have the most information that isn’t present in the other variables (Figure 7).

3.4. Modeling for ZCL vector and reservoir hosts, Horizon 2050: RCP 8.5 scenario, BCC-CSM1-1(m) model, P.papatasi

Estimates of relative contributions of the environmental variables to the MaxEnt model for the vector and reservoir hosts of ZCL have 1.00 0.92 0.85 0.77 0.69 0.62 0.54 0.46 0.38 0.31 0.23 0.15 0.08 0 1.00 0.92 0.85 0.77 0.69 0.62 0.54 0.46 0.38 0.31 0.23 0.15 0.08 0

been presented in Table 2. As shown in the Table 2, the variables with the most contributions to the models for P. papatasi and R.opimus were bio8 (53.0%) and bio7 (60.6%), respectively. The AUC curve calculated for P. papatasi was 0.891 (Figure 8). Analysis of jackknife test of variable importance for P. papatasi showed that the environmental variable with highest gain when used in isolation is bio8 , which therefore appears to have the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted is bio8 , which therefore appears to have the most information that isn’t present in the other variables (Figure 8).

Figure 6. The receiver operating characteristic curve and output of MaxEnt model and jackknife analysis for Phlebotomus papatasi in Yazd province, Central Iran, Horizon 2050： RCP 4.5 scenario, BCC-CSM1-1(m) model.

Figure 7. The receiver operating characteristic curve and output of MaxEnt model and jackknife analysis for Rhombomys opimus in Yazd province, Central Iran, Horizon 2050： RCP 4.5 scenario, BCC-CSM1-1(m) model.

The AUC curve calculated for R. opimus was 0.839 (Figure 9).Jackknife test of variable importance for R. opimus showed that the environmental variable with highest gain when used in isolation is bio7, which therefore appears to have the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted is bio7, which therefore appears to have the most information that isn’t present in the other variables (Figure 9).The environmentally suitable area for both vector and reservoir in horizon 2030 using RCP 4.5 and RCP 8.5 scenarios is presented in Figure 10.

Table 2. Contribution of variables used to predict the potential distribution of Phlebotomus papatasi and Rhombomys opimus as vector and main reservoir of ZCL in Yazd Province, Central Iran, horizon 2050： RCP4.5 and RCP8.5 scenarios, BCC-CSM1-1 (m) model.

3.5. Climatic suitability for vector and reservoir hosts of ZCL in Yazd province

The results of the projected (2030 & 2050) climatic suitability for the two species (P. papatasi and R. opimus) are shown in Table 3.Western and central part of Yazd province will provide more suitable climatic conditions for P. papatasi and R. opimus to the end of the horizon 2030 and 2050.

Table 3. Projected climatic suitability (%) and the area (km2) covered by different classes for Phlebotomus papatasi and Rhombomys opimus in Yazd province, center of Iran.

4. Discussion

To investigate the potential changes in the status of ZCL vector and reservoir species in the climate zone of Yazd province in the future, the probability of the presence of P. papatasi and R. opimus in 2030 and 2050 were examined. Probability of the presence of the two studied species was grouped into five classes using equal intervals in ArcGIS. Areas with more than 60% presence probability of species were considered as areas with high potential of CL transmission. In Horizon 2030, by using RCP 4.5 scenario, we found that the presence probability of the vector within 7 017 km2from Yazd province is above 60%. In this scenario, the variable with the highest contribution to the models for P. papatasi was bio8. Variation in this variable from about 70 to 140 has had a negative trend in the presence of P. papatasi, and then it does not change. The value of AUC was 0.917, which indicates that the output of the model is valid.

RCP 4.5 scenario revealed a potential expansion in the habitat of R.opimus to 10 712 km2from Yazd province. Jackknife test of variable importance for R. opimus showed that bio7 has the maximum effect on the dispersion of this species. Variation in this variable from about 395 to 445 has had a negative trend in the presence of R. opimus,and then it does not change. The value of AUC was 0.851, which indicates that the output of the model is valid.

In Horizon 2030, RCP 8.5 scenario revealed a potential dispersion of P. papatasi within 8 817 km2from Yazd province at various probabilities. The results of jackknife test show that bio8 has the highest effect on the dispersion of this species. Variation in this variable from about 0 to 160 has had a negative trend in the presence of P. papatasi, and then it does not change. The calculated AUC was 0.889, which indicates that the output of the model is valid.

In this model, we observed a potential presence of R. opimus within 11 403 km2from Yazd province. Jackknife analysis for R. opimus showed that bio7 has the highest gain among the environmental variables. Variation in this variable from about 390 to 445 has had a negative trend in the presence of R. opimus, and then it does not change. The value of AUC was 0.794 indicating that the output of the model is valid. In Horizon 2030, in general, we observed a higher expansion of habitat suitability of ZCL vector and reservoir hosts in scenario 4.5 compared with scenario 8.5.

Figure 8. The receiver operating characteristic curve and output of MaxEnt model and jackknife analysis for Phlebotomus papatasi in Yazd province, Central Iran, Horizon 2050： RCP 8.5 scenario, BCC-CSM1-1(m) model.

Figure 9. The receiver operating characteristic curve for Rhombomys opimus and output of MaxEnt model and jackknife analysis for Rhombomys opimus in Yazd province, Central Iran, Horizon 2050： RCP 8.5 scenario, BCC-CSM1-1(m) model.

Figure 10. A： Predicted environmental suitability for vector and reservoir of zoonotic cutaneous leishmaniasis in Yazd province of Iran using RCP 4.5 scenario, BCC-CSM1-1(m) model, Horizon 2050： vector (left) and reservoir (right). B： Predicted environmental suitability for vector and reservoir of zoonotic cutaneous leishmaniasis in Yazd province of Iran using RCP 8.5 scenario, BCC-CSM1-1(m) model,Horizon 2050： Vector (left) and Reservoir (right).

In Horizon 2050, RCP 4.5 scenario revealed a potential dispersion of the vector within 7 953 km2from Yazd province at various probabilities.The results of jackknife test of variable importance to the models for P.papatasi show that the environmental variable with highest gain when used in isolation is bio8. Variation in this variable from 20 to 115 has had a positive trend in the presence of P. papatasi, but in higher values up to about 175 has had a negative effect, and after that it does not change. The calculated value of AUC was 0.903, which indicates that the output of the model is valid. In this scenario, there is a potential presence of R. opimus within 11 370 km2from Yazd province. The results of jackknife test of variable importance for R. opimus show that the environmental variable with highest gain when used in isolation is bio7. Variation in this variable from about 385 to 450 has had a negative trend in the presence of R. opimus, and then it does not change. The value of AUC was 0.849 indicating that the output of the model is valid.

In Horizon 2050, by using RCP 8.5 scenario, we observed that there is a potential presence of P. papatasi within 7 440 km2of Yazd province. The results of jackknife test of variable importance for P.papatasi indicate that the environmental variable with the highest gain when used in isolation is bio8 . Variation in this variable from about 20 to 180 has had a negative trend in the presence of P.papatasi, and then it does not change. The value of AUC was 0.891,which indicates that the output of the model is valid.

In this Horizon, there is a potential dispersion of R. opimus within 11 817 km2of Yazd province. Jackknife test of variable importance 1.00 0.92 0.85 0.77 0.69 0.62 0.54 0.46 0.38 0.31 0.23 0.15 0.08 0 1.00 0.92 0.85 0.77 0.69 0.62 0.54 0.46 0.38 0.31 0.23 0.15 0.08 0 for R. opimus revealed that the environmental variable with the highest gain when used in isolation is bio7. Variation in this variable from about 395 to 445 has had a negative trend in the presence of R.opimus, and then it does not change. The calculated value of AUC was 0.839, indicating that the output of the model is valid.

According to the prediction of the model in the study area, in both 2030 and 2050 there is a decrease in the environmental suitability for P. papatasi, while the condition will increase the suitable areas for R. opimus, although the changes are not so big. Previous studies conducted in other countries revealed a significant relationship between CL and climatic factors. Toumi et al in Sidi Bouzid reported that only average temperature was not a statistically significant predictor of ZCL incidence[35]. Contrary to our findings, in another study, Temperature Annual Range & Mean Temperature of Wettest Quarter (℃) had the greatest effect on the presence of ZCL vector and reservoir, respectively. The results of our study are consistent with the findings of other researchers around the world including Roger et al in South America and Sing in Rajasthan, India. Both studies showed that the incidence of the disease increased with rise of temperature and decreased with decline of relative humidity[36,37].

In a study conducted by Yazdanpanah and Rostamianpur[38] in Ilam Province in the West of Iran, they found a positive correlation between average temperature and CL incidence, whereas this correlation was negative in another study conducted by Mozafari and Bakhshizadeh[39] in Yazd-Ardakan plain in central of Iran. In contrast to our findings, both studies conducted in Ilam and Yazd reported that the association between CL epidemic and relative humidity was not significant, whereas in our study, relative humidity appeared to be the most significant factor.

Climate factors are the major determinants of ZCL incidence rate in Yazd Province. These climate conditions provide the favorable habitats for propagation and transmission of ZCL in this endemic area. Habitat suitability for the vector and reservoir hosts in the coming years will be expanding such that in horizon 2030 and 2050 the probability presence of the vector and reservoir hosts in 38 580 and 37 949 km2, respectively, from Yazd province is above 60%. Also, we predict an increase in the presence of the vector in the western parts and the reservoir in the northern and central parts of the province in the future. In the present study, we observed a complex relationship between environmental factors, bioclimatic factors and ZCL occurrence in Yazd province. Understanding the role of environmental and bioclimatic factors in ZCL occurrence can provide a guide for policy-makers in the creation and implementation of more effective policies for the prevention and control.

Conflict of interest statement

Authors declare no conflict of interest.

Acknowledgements

The authors are grateful to the staff of Yazd Province Health Center for the kind collaboration in this study.

Foundation project

This study has been financially supported by Research Deputy,Shahid Sadoughi University of Medical Sciences, Yazd, Iran, with project No. P/17/1/19412.