Population structures of Leishmania infantum and Leishmania tropica the causative agents of kala-azar in Southwest Iran

Visceral leishmaniasis (VL) is endemic in Iran and is caused predominantly by Leishmania infantum, but L. tropica is emerging as an important cause. We studied the intra-species population structure of Leishmania spp. causing VL in southwest Iran by sequence analysis of the internal transcribed spacer (ITS) 1 of DNA samples from 29 bone marrow aspiration smears. L. infantum (n = 25) and L. tropica (n = 4) were identified, consisting of 10 and three ITS1 sequence types (STs), respectively. Compared to GenBank ITS1 STs, our L. infantum parasites displayed high heterogeneity but less heterogeneity compared than northwest Iranian isolates. VL affects mostly nomadic populations in southwest Iran, and their mobility may explain partly the L. infantum heterogeneity. The VL causing L. tropica was also genetically heterogeneous but genetically indistinguishable from L. tropica strains causing anthroponotic cutaneous leishmaniasis from southwest Iran.


Introduction
Leishmaniasis is a complex of neglected tropical diseases, caused by different species of Leishmania, which are present in five continents (Zijlstra 2016). Visceral leishmaniasis (VL) is the most severe form of disease and is potentially fatal if untreated. Globally, the annual VL incidence is approximately 200,000 to 400,000 cases with 20,000 to 40,000 deaths per year (Alvar et al. 2012).
VL causing species vary geographically, Leishmania donovani in South Asia (Bern et al. 2007;Rijal et al. 2010;Singh et al. 2010) and East Africa (Elnaiem 2011;Zijlstra 2016), and Leishmania infantum (L. chagasi) in Central and South America, the Middle East, Mediterranean Basin, Central Asia and China (WHO 2010). VL caused by Leishmania tropica, a species that causes Old World cutaneous leishmaniasis (CL), has been reported sporadically in humans in South Asia and humans and dogs in the Middle East and North Africa (Mebrahtu et al. 1989;Guessous-Idrissi et al. 1997;Alborzi et al. 2006;Hajjaran et al. 2007;Khanra et al. 2012;Hosseininasab et al. 2014).
In Iran, L. infantum is the main cause of VL. However, L. tropica and L. major, the organisms causing anthroponotic and zoonotic cutaneous leishmaniasis (A and ZCL), respectively, may also cause rarely VL with L. tropica being more common than L. major (Karamian et al. 2007;Mohebali 2013). The principal VL foci in Iran are in the northwest (NW) provinces of East Azerbaijan and Ardabil and southwest (SW) provinces of Fars, Bushehr, Kohgiloye and Boyerahmad and Kerman; most VL cases are reported from Fars (Sarkari et al. 2012;Sarkari et al. 2016).
Understanding the population structure and phylogenetic divergence of Leishmania spp. at the intra-species level in time and place could provide way of tracking epidemiological changes in leishmaniasis (Ghatee et al. 2014). Multilocus isoenzyme electrophoresis (MLEE) has long been the gold standard for identifying Leishmania species and intra-species differences (Pratlong et al. 1991). However, because of its lower intra-species discriminatory power and the need for parasite mass cultivation, it has been superseded by other techniques as multilocus sequence typing (MLST), multilocus microsatellite typing (MLMT), polymerase chain reaction restriction fragment length polymorphism (PCR-RFLP) and sequencing of specific genetic markers such as hsp70, ribosomal DNA (rDNA) and the kinetoplastid DNA (kDNA) (Schönian et al. 2008Yuan et al. 2017;Ghatee et al. 2018b). Of these, MLMT and kDNA-RFLP are the most powerful techniques for studying Leishmania populations in intraspecies level , Schönian et al. 2013). However, challenges remain. MLMT and MLST have been confined in developing countries due to expensive cost. Aside from kDNA, PCR-RFLP of genes has not shown sufficiently high resolution for population structure studies, but kDNA-RFLP is associated with variable DNA yields and low reproducibility of results (Bhattarai et al. 2010;Downing et al. 2011;Ghatee et al. 2014).
ITS1 has greater polymorphism than ITS2 and shown greater discrimination of L. donovani and L. tropica in phylogenetic studies in Sudan and Iran, respectively (El Tai et al. 2001;Ghatee et al. 2014). Moreover, ITS1 is more frequently deposited in GenBank, and this allows a greater opportunity for comparing genotypes. There are no data on the genetic population structure of VL causing Leishmania spp. from southern Iran. Therefore, we conducted an ITS1 sequence analysis of VL causing Leishmania spp. in southwest Iran and compared the population structure with those from other regions of Iran.

Study area and patients
Bone marrow aspiration smears from 29 patients with microscopically confirmed kala-azar were obtained from the pathology laboratories of Shahid Faghihi (SF) hospital (n = 26), Shiraz, the capital of Fars Province in SW Iran and Afzalipour hospital (n = 3), Kerman City, the capital of Kerman Province in south-southeast Iran. The SF hospital is a referral centre and receives patients from neighbouring provinces. Patients' home addresses were retrieved from the hospital records. The VL patients were from Fars, Bushehr, Kohgiloyeh and Boyerahmad Provinces, the kala-azar endemic regions in southwest Iran and two counties in southsoutheast province of Kerman in the close proximity of Fars Province. These provinces are also endemic foci of CL. This study was approved by the Ethic Committee of Yasuj University of Medical Sciences (IR.YUMS.1.REC.1395.3), and patients' record were anonymized and de-identified prior to analysis.

Microscopic examination & DNA extraction
Bone marrow aspiration smears were re-examined microscopically to reconfirm the presence of Leishman bodies and were scratched using sterile scalpels. The sample scrapings were placed in 1.5 ml microtubes containing lysis buffer (Tris 100 mM, EDTA 10 mM, NaCl 100 mM, SDS 1%, Triton X100 2%). Ten micrograms per millilitre of proteinase K was added to the mixture and was incubated at 56°C for 1 h and then extracted with phenol/chloroform (25:24 v/v) first and then again with chloroform. DNA was precipitated with equal volumes of iso-propanol and one tenth volume of 3 M NaAc. The deposited DNA was washed with 70% ethanol, dried, and suspended in 30 μl of ultrapure water.

Kinetoplastid DNA PCR
To identify the species of the Leishmania isolates, kDNA was amplified using the primers 13Z (5′-ACT GGG GGT TGG GTG TAA AATAG-3′) and LiR (5′-TCG CAG AAC GCC CCT-3′) (Noyes et al. 1998). The PCR mixture contained 12.5 μl of 2× premix (Ampliqon, Denmark), 20 pmol of each primer, 5 μl of template DNA, and molecular grade water up to 25 μl. Reference strains and negative and inhibition controls were used for each experiment. The PCR process included a pre-denaturing step at 95°C for 5 min followed by 35 cycles at 94°C for 45 s, 55°C for 60 s and 72°C for 90 s in an Applied Biosystems thermocycler. The PCR products were electrophoresed in 1.2% agarose gel with 0.5 μg/ml ethidium bromide for 90 min at 80 Vin TBE buffer and visualised by a transilluminator. A 100-bp DNA marker was used in each experiment as the standard size. Previously defined isolates of L. major (MRHO/IR/75/ER) and L. infantum (MCAN/IR/ 07/Moheb-gh) and also L. tropica from a clinical isolate with a laboratory confirmed species finding, preserved in RPMI 1640 culture media, were used as reference strains.

ITS1-rDNA PCR
ITS1-rDNA was amplified using primers LITSR (forward, 5′-CTG GAT CAT TTT CCG ATG-3′) and L5.8S (reverse, 5_-TGA TAC CAC TTA TCG CAC TT-3′), as described by Schönian et al. (2003). The PCR mixture consisted of 25 μl of 2× premix (Ampliqon, Denmark), 20 pmol of each primer, 10 μl (50 ng) of template DNA, and molecular grade water up to 50 μl. Negative and inhibition controls were used for each PCR run. The cycling PCR conditions were 95°C for 5 min followed by 35 cycles of 94°C for 30 s, 53°C for 30 s and 72°C for 60 s, and a final elongation step at 72°C for 5 min. The PCR products were electrophoresed on 1.5% agarose gel electrophoresis with 0.5 μg/ml of ethidium bromide for 60 min at 80 V and visualised under the same conditions, as described above.

Sequence analysis
The PCR products of ITS1 DNA were excised from the gel and purified using a gel purification kit (Bioneer, South Korea), according to the manufacturer's instructions, and were sent to Bioneer Company (South Korea) for sequence analysis using an Applied Biosystems 3730 XL automated DNA sequencer. Sequencing was performed in both directions using by the same primers used for ITS1 amplification. The ITS1 sequences were deposited in the GenBank database under the accession numbers KY964622-KY964634.

Phylogenetic analysis
Species identification, based on the kDNA-PCR, was confirmed by comparison of our ITS1 sequences with those published in the GenBank database using BLAST software (http:// www.ncbi.nlm.nih.gov). The obtained sequences were analysed by Geneious Pro 5.5.6 software (http://www. geneious.com, Kearse et al. 2012). Phylogenetic trees were generated by using L. infantum and L. tropica ITS1 sequence types (ST) from this study and those from GenBank, namely L. infantum from Iran and L. tropica from Iran and west Afghanistan (the eastern neighbour of Iran). The maximum likelihood (ML) tree was constructed using MEGA 6 software (Tamura et al. 2013) after trimming all sequences at both ends. Bootstrap values for the ML method were based on 1000 replicates. Figure 1 shows the provinces from where sequence types were obtained.

Results
The patient's demographic data are shown in Table 1. Based on the species-specific size, 25 and four samples were identified as L. infantum (680 bp) and L. tropica (750 bp), respectively, by using kDNA PCR. The species identification was confirmed by comparing our ITS1 sequences with sequences that have been deposited in GenBank. Fig. 1 The provinces where ITS1 sequences of L. infantum (a) and L. tropica (b) were obtained and retrieved from GenBank are shown in yellow and orange, respectively. The green circles show provinces that our L. infantum and L. tropica were isolated from in current study. Most isolates were from Fars Province. In map A, Baft County from Kerman Province was shown Sequence types-L. infantum Ten haplotypes were found among the 25 L. infantum ITS1 sequences. ST1 was the most frequent [n = 10 (40%)], followed by ST4 (n = 4), ST2 (n = 3), ST3 (n = 2), and one each for ST5-ST10. All nucleotide variations were found in the first 120 upstream nucleotides, which occurred mostly in the first 60 nucleotides and in the microsatellite sequences, including CC, AA, and TA repeats (Fig. 2). Most similarity was found between ST1 and ST4 (99.6%), ST4 and ST7 (99.2%), and ST4 and ST3 (99.2%), and the least similarity was found between ST10 and ST8 (92.7%), and ST10 and ST6 (93%) (Fig. 3).
Sequence types-L. tropica Three STs were found among the four L. tropica parasites. STA was found in two isolates and STB and STC in the remaining two. Greater similarity was present between STA and STB in comparison to STA or B vs. STC. Deletion of CA was found in the region with repeats of C in the upstream part of STB in comparison to STA, but STC showed more differences in upstream and downstream regions that are shown in Fig. 4.

Phylogenetic trees
The L. infantum STs from different VL foci in Iran (our data and those from GenBank) were distributed in different positions in the maximum likelihood tree. Clade A-i comprised all STs from SW (Fars, Bushehr, Kohgiloyeh and Boyerahmad Provinces), SE (Kerman Province), central (Alborz, Tehran and Yazd Provinces), west (Lorestan Province) and northeast (Golestan Province) Iran and as well as some STs from northwest (NW) Iran (Ardabil Province). In this clade, STs1, 3-4, 7 and 9-10 from our SW isolates made up a subcluster. Additional sub clusters consisted of (i) all ST sin clade A-i from central, NW and NE Iran together with ST2 and a haplotype (KC570454) from SW Iran (GenBank), and (ii) ST8 and one haplotype from west Iran (Lorestan Province). STs 5 and 6 formed separate branches within clade A-i. Clade B-i consisted of only two STs from NW Iran (Ardabil Province) (Fig. 6).

Discussion
In current study, we have shown that VL in patients from southern Iran was caused mostly by L. infantum and, in a small number of patients, by L. tropica. Ten L. infantum and three L. tropica STs were identified among the 29 ITS1 sequences.
We found high heterogeneity for L. infantum from SW Iran, but these were mostly grouped together on the phylogenetic tree and included the dominant STs. NW Iranian L. infantum isolates had the highest heterogeneity, and some of these isolates formed a distinct clade (B-i), which showed the most divergence compared to other Iranian L. infantum isolates, including those from the NW in clade A-i. The least heterogeneity was found in L. infantum parasites from Central Iran. The VL-associated L. tropica isolates were also heterogeneous, but they were not genetically distinct from other L. tropica isolates that caused ACL in SW Iran.
The clade B-i, NW Iranian L. infantum isolates (GenBank: EU637915 and EU680963), have previously been detected from Phlebotomus perfilewi, the main vector of L. infantum, and dogs in Ardabil Province of NW Iran (Moshfe et al. 2008;Oshaghi et al. 2009) but, interestingly, are also identical to ITS1 from L. donovani (KT438661-2, KT4386674-9, KT438680-1). Confirmation of L. infantum and differentiation from L. donovani was also shown by cysteine protease B (CBP) gene amplification by Oshaghi et al. (2009), in the previous study. Moreover, L. donovani has also been isolated from P. perfilewi sand flies in NW Fig. 5 The similarity of three L. tropica sequence types obtained from visceral leishmaniasis patients from southwest Iran compared to GenBank sequence types from patients with anthroponotic cutaneous leishmaniasis due to L. tropica from southwest, southeast, and eastern Iran. ST1-7 (east) ) Hap A, B, and D (southeast), Hap A and C (southwest) (Ghatee et al. 2014)  Iran and Meriones persicus, the Persian jird (a rodent), from several areas of Iran (Mohebali et al. 2004;Oshaghi et al. 2009). Iran has not been known as a focus of L. donovaniassociated VL, but these data suggest that Iran could support the transmission of L. donovani.
In the Middle East, L. donovani is confined to foci in Turkey, Cyprus and Iraq (WHO 2010; Gouzelou et al. 2012). In the Cukurova region of SE Turkey, bordering NW Iran, there is a major focus of CL caused by L. donovani/ infantum hybrid strains, which are transmitted by P. tobbi, the main sandfly species in Cukurova (Seblova et al. 2015). P. tobbi has also been identified as one of the vectors of L. infantum in NW Iran (Rassi et al. 2012). Based on the results of sequence analysis of the ITS1 and Cbp loci and the epidemiological similarities between NW Iran and SE Turkey, we hypothesise that distinct strains of L. infantum from NW Iran may have evolved from the hybrid strains of L. donovani/L. infantum in neighbouring SW Turkey, and this may explain why they formed their own clade in the phylogenetic tree. More work is needed to confirm our hypothesis.
High heterogeneity was also shown in SW Iranian L. infantum haplotypes and may be explained by population movement. Indeed, VL in SW Iran may have begun because of the immigration of nomadic populations from NW Iran in the sixteenth century (Mazloumi Gavgani et al. 2011). Nomadic populations have the highest VL seropositivity rates in Iran (Mohebali et al. 2006) and their movements, especially in semi-arid regions where sand flies thrive, and the proximity of villages to nomadic travel routes are key factors affecting the distribution of VL in Iran (Ghatee et al. 2013a, b). We hypothesise that the high heterogeneity of SW Iranian L. infantum in clade A-i also may be explained by the recent L.  Fig. 6 Unrooted maximum likelihood tree inferred from the internal transcribed spacer 1 sequences of 10 representatives of L. infantum from 25 isolates from southwest Iran and those from other parts of Iran retrieved from GenBank. Only one ITS1 L. infantum from Southwest Iran (KC57454) was previously deposited in GenBank. The numbers above the branches show the bootstrap value for 1000 replicates co-evolution. Evidences of parasite-vector co-evolution for L. tropica and Ph. sergenti strains, the proven vector of L. tropica, were reported in south Iran .
VL occurs sporadically in other regions of Iran (Hamzavi et al. 2012) and is emerging in foci in central Iran (Heidari et al. 2015). The similarity of sequence haplotypes from different  Fig. 7 Unrooted maximum likelihood tree inferred from the internal transcribed spacer 1 sequences of three L. tropica isolates from patients with visceral leishmaniasis from southwest (STA also from southeast) Iran and L. tropica sequences from other regions in Iran and west Afghanistan regions of central Iran can be explained by the clonal expansion of newly established Leishmania parasites, the presence of few nomadic tribes and similar ecosystems in these foci.
CL due to L. infantum has only recently been identified in SW (Motazedian et al. 2002) and NW Iran (Badirzadeh et al. 2013), but the phylogenetic position of these isolates has not been studied until now.
L. tropica causing VL was first reported in a patient in Kenya (Mebrahtu et al. 1989) and has been followed by other human cases from American soldiers deployed in Iraq and individuals from Iran, India, and Kenya and Afghanistan (Magill et al. 1993;Sacks et al. 1995;Alborzi et al. 2006;Weiss et al. 2009;Khanra et al. 2012). India reports the highest caseload of L. tropica-related VL (Sacks et al. 1995). Canine VL due L. tropica has been reported from Iran and Morocco (Guessous-Idrissi et al. 1997;Hajjaran et al. 2007;Mohebali et al. 2011).
The first case of L. tropica VL came from SW Iran (Alborzi et al. 2006), and other cases have been reported countrywide with varying frequencies, e.g., ranging from 1.4% in 1993/4 to 14% in our study (Alborzi et al. 2006;Mohebali 2013;Hosseininasab et al. 2014). Of our four cases, two were obtained from Shiraz and Kerman City, which are well established endemic foci of L. tropica-related ACL (Ghatee et al. 2014;Izadi et al. 2016), while the other two were from areas where L. tropica has not been reported previously. These limited data hint at trend of increasing L. tropica-related VL and emergence in the new regions in Iran that may be related to population movements or global warming that promotes favourable breeding of sand flies (Ghatee et al. 2018a). More work is needed to establish the trend of L. tropica-related VL to inform future control strategies.
Homogeneity of SE and most east and central Iranian and western Afghanistan and the genetic heterogeneity of northeast and SW Iranian ACL-associated L. tropica have been described previously (Ghatee et al. 2014;Fakhar et al. 2016;Karamian et al. 2016). Viscerotropic haplotypes, STA and STB, are in the same position on the phylogenetic tree as the SW ACL isolates and that the STC haplotype clustered with other SW ACL isolates. Overall, our three L. tropica haplotypes from viscerotropic cases could not diverge from ACL isolates from SW Iran. Indian VLassociated L. tropica isolates were also not discerned from the Indian CL-associated subpopulation (Krayter et al. 2014). The visceralisation of L. tropica may be related partially to reduced host immunity, e.g., HIV positivity (Jafari et al. 2010) and/or increased pathogenicity of strains (Krayter et al. 2014). HIV and HIV/VL co-infection rate seems to increase in Iran where HIV/VL occurrence was reported up to 18% of HIV-positive patients (Shafiei et al. 2014;Gökengin et al. 2016) that may be led to higher incidence of the visceralisation of cutaneous or asymptomatic Leishmania infections. This is an area where greater understanding of L. tropica viscerotropism is needed.

Conclusion
We found 10 heterogeneous L. infantum and three heterogeneous L. tropica STs from patients with VL from SW Iran. The high heterogeneity of SW and NW L. infantum may be due to within-country population movements. Central Iran has the least heterogeneity where only small VL foci exist. L. tropica causing VL from SW Iran was heterogonous but was not distinct genetically from ACL-related L. tropica from SW Iran. Mechanisms underlying the visceralisation of L. tropica remain unknown. More research is needed to understand the dynamics of VL-related L. infantum and L. tropica and to understand the adaptation of L. tropica to visceral disease.