Evolutionary history of Leishmania killicki (synonymous Leishmania tropica) and taxonomic implications

The taxonomic status of Leishmania (L.) killicki, a parasite that causes chronic cutaneous leishmaniasis, is not well defined yet. Indeed, some researchers suggested that this taxon could be included in the L. tropica complex, whereas others considered it as a distinct phylogenetic complex. To try to solve this taxonomic issue we carried out a detailed study on the evolutionary history of L. killicki relative to L. tropica. Thirty-five L. killicki and 25 L. tropica strains isolated from humans and originating from several countries were characterized using the MultiLocus Enzyme Electrophoresis (MLEE) and the MultiLocus Sequence Typing (MLST) approaches. The results of the genetic and phylogenetic analyses strongly support the hypothesis that L. killicki belongs to the L. tropica complex. Our data suggest that L. killicki emerged from a single founder event and that it evolved independently from L. tropica. However, they do not validate the hypothesis that L. killicki is a distinct complex. Therefore, we suggest naming this taxon L. killicki (synonymous L. tropica) until further epidemiological and phylogenetic studies justify the L. killicki denomination. This study provides taxonomic and phylogenetic information on L. killicki and improves our knowledge on the evolutionary history of this taxon.


Background
Leishmaniases are neglected tropical diseases caused by Leishmania parasites and transmitted to mammals through bites by infected Phlebotomine sandflies of the genus Phlebotomus [1]. In humans, these diseases can have cutaneous (CL), muco-cutaneous (MCL) or visceral (VL) clinical manifestations.
Since the first description of the genus Leishmania Ross, 1903, the classification methods have considerably evolved. Indeed, between 1916 and 1987, Leishmania taxonomy followed the Linnaean classification system, mainly based on extrinsic features, such as clinical manifestations, geographical distribution, epidemiological cycles and behaviour in sandfly vectors. This method has led to the subdivision of the Leishmania genus in the two sub-genera Leishmania and Viannia [2,3].
In the eighties, the biochemical classification based on the study of the parasite isoenzymatic patterns started to be developed. This approach has evolved from the classical Adansonian to the numerical cladistic classification method that uses isoenzymes as evolutionary markers [4][5][6][7][8]. The description of several Leishmania complexes in the Old and New World is based on these analyses. Specifically, by using numerical phenetic and phylogenetic approaches, Rioux et al. [9] identified four main Leishmania groups in the Old World, while Thomaz et al. and Cupolillo et al. [10,11] defined eight complexes and two Leishmania groups in the New World.
Currently, the numerical taxonomic approach based on isoenzyme analysis is considered as the gold standard for the classification of the genus Leishmania and is routinely used for classification updates and for epidemiological studies [12,13]. The drawbacks of this approach are the need of bulk cultures of Leishmania parasites and its relatively poor discriminatory power. It is also time-consuming. Therefore, DNA-based techniques represent valuable alternatives for the identification and the classification of these parasites.
MLST is one of the most appropriate approaches for taxonomic studies because it provides data on the genetic variations of housekeeping genes. This approach has been increasingly used for phylogenetic investigations to understand the epidemiological and transmission features of many Leishmania complexes [20,[29][30][31][32][33]. However, because of the complexity of this genus and the lack of studies, several taxa need to be detailed further [34].
Leishmania killicki is a recently described taxon that causes CL in Tunisia [35], Libya [36] and Algeria [37]. L. killicki taxonomic status and evolutionary history relative to L. tropica are based on very few studies and samples. The numerical taxonomic analysis using the Multilocus Enzyme Electrophoresis (MLEE) approach first included this parasite in the L. tropica complex [9,38]. However, after the revision of the Leishmania genus classification, it was considered as a separate phylogenetic complex [39]. Recently, an update study by Pratlong et al. [12] confirmed the inclusion of L. killicki within the L. tropica complex. Phenetic and phylogenetic studies using MLMT [40], PCR-sequencing [41] and MLST [31] also classified L. killicki within the L. tropica complex and suggested a closer genetic link with L. tropica from Morocco. However, these data were obtained using only seven L. killicki strains: two strains were analyzed by Schwenkenbecher et al. [40], two by Chaouch et al. [41] and three by El Baidouri et al. [31]. Therefore, the present study wanted to analyze by MLST a large number of L. killicki and L. tropica strains in order to precisely determine the evolutionary history and the taxonomic status of L. killicki.

Origin of strains
For this study, strains of L. killicki (n = 35), L. tropica (n = 25), L. major (n = 1) and L. infantum (n = 1) from different geographic areas and with various zymodeme patterns were included (total = 62 strains). These strains were from human cutaneous lesions, except the L. infantum strain that was isolated from a patient with VL. Most strains (n = 53) were selected from the Cryobank of the Centre National de Référence des Leishmanioses (CNRL) (Montpellier, France) and nine L. killicki strains were collected by the team of the Laboratoire de Parasitologie -Mycologie Médicale et Moléculaire (Monastir, Tunisia) during epidemiological investigations.
Forty-eight strains, among which 34 L. killicki strains (six from Algeria, one from Libya and 27 from Tunisia) and 14 L. tropica strains from Morocco were analyzed by MLST for the first time during this study. The eleven remaining L. tropica strains were from several countries (one from Egypt, one from Greece, two from Israel, two from Jordan, three from Kenya and two from Yemen) and were previously typed by MLST. Their sequences were published in Genbank under the following accession numbers: KC158621, KC158637, KC158643, KC158677, KC158682, KC158683, KC158690, KC158696, KC158711, KC158722 and KC158761 (see [31]). One L. killicki strain (LEM163) MHOM/TN/80/LEM163 had also already been analyzed by MLST (Genbank accession number KC158820 (see [31]).

Isoenzymatic identification
All studied strains were identified by MLEE, according to Rioux

DNA extraction
Genomic DNA from cultured parasites was extracted using the QIAamp DNA Mini Kit (Qiagen, Germany) following the manufacturer's recommendations and eluted in 150 μl.

Analysis by Multilocus sequence typing (MLST)
The L. killicki (n = 34) and L. tropica (n = 14) strains that had not been previously assessed by MLST were typed using the MLST approach based on the analysis of seven loci coding for single-copy housekeeping genes that was developed and optimized by El Baidouri et al. [31]. Genomic DNA was amplified by real-time PCR using the SYBR Green method (Light cycler 480 II, Roche). The amplified products were sequenced on both strands (Eurofins MWG Operon, Germany) and the obtained sequences were aligned and checked in both directions using the CodonCode Aligner software, v.4.0.1 (Codon Code Co., USA). For each strain, polymorphic sites (PS) and ambiguous positions corresponding to heterozygous sites (HS) were identified in each locus using the same software. The DnaSP software v.5 [42] was used to calculate the number of haplotypes from the concatenated sequences.  Phylogenetic relationships were inferred using a Bayesian approach implemented with the MrBayes software v. 3.2.3 [43]. The concatenated duplicated sequence alignments of the seven loci for the 32 Leishmania strains representing all the identified haplotypes and the two outgroup strains (n = 34 in total) were used to run two independent chains for 10,000,000 generations each and trees sampled every 1000 generations. The burn-in period was set to 200,000 generations to fit the first 20% of the analyses. Analyses were conducted using the General time reversible model of substitution with a proportion of invariable sites and gamma distribution estimated by the program (GTR + G + I).
The chain convergence was assessed using the average standard deviation of split frequencies (ASDSF). If two runs converge onto the stationary distribution, the ASDSF is expected to approach zero, reflecting the fact that the two tree samples become increasingly similar. An average standard deviation below 0.01 is thought to be a very good indication of convergence (below 0.004 in our analysis). The consensus tree was constructed using 1000 trees sampled from the stationary phase. The MEGA 5.10 software [44] was used to identify amino acid variations between L. killicki and L. tropica.

Isoenzymatic identification of Leishmania strains
Among the 62 strains under study, 53 had been previously characterized by MLEE at the Centre National de Référence des Leishmanioses. The nine strains collected by the team of the Laboratoire de Parasitologie -Mycologie Médicale et Moléculaire (Monastir, Tunisia) were identified for the first time in this study using the same technique [12,[35][36][37][38][45][46][47]. Nevertheless, all the strains were analyzed again by MLEE at the Centre National de Référence des Leishmanioses (Montpellier, France).
Seventeen zymodemes were identified: three for L. killicki, 12 for L. tropica and a single zymodeme for each L. major and L. infantum strain (Table 1) [35]. On the other hand, the MDH, ME, GOT1, GOT2 and FH profiles were different in the zymodemes MON-317 and MON-301 [37], and the MDH, GOT1, GOT2 and FH profiles allowed discriminating between MON-317 and MON-306 (a zymodeme described in Algeria, but not included in our sample collection) [48] (Table 2). For L. tropica, all the identified zymodemes were already known [12,45] (Table 1).

Sequence analysis
The sequences of the L. killicki (n = 34) and L. tropica (n = 14) strains were submitted to GenBank (accession numbers from KM085998 to KM086333). The sizes of the seven loci under study were identical to those reported by El Baidouri et al. [31], except for locus 12.0010 (only 579 pb instead of 714 pb), leading to a total length of 4542 pb for the concatenated sequences (Table 3). All chromatograms were clearly readable. Polymorphic sites (PS) and heterozygous sites (HS), which corresponded to ambiguous positions with two peaks, were easily identified. No tri-allelic site was detected.
Genetic polymorphisms in L. killicki and in L. tropica  Table 4).
Assessment of the presence of mutations in the seven loci under study in the L. killicki and L. tropica (heterozygous mutations were excluded from the analysis)   identified 55 mutations of which 29 were silent substitutions and 26 resulted in altered amino acid residues ( Table 5). All L. killicki mutations corresponded to a single amino acid change. Conversely, in the L. tropica strains, mutations could lead to more than one amino acid change.
Phylogenetic analysis of L. killicki In total, 32 different haplotypes were identified: 10 for the 35 L. killicki strains and 22 for the 25 L. tropica strains. Twenty-six haplotypes were unique (eight for L. killicki and 18 for L. tropica) and the two taxa did not share any haplotype. The L. killicki MON-317 (strain LEM6173) had its own haplotype ( Table 6). The Bayesian consensus tree using 32 strains representing all the identified haplotypes was constructed based on the concatenated sequences and duplicated nucleotide sites to avoid the loss of genetic information in ambiguous positions ( Figure 1). The phylogenetic tree showed that L. killicki formed a separate group, although it belonged to the L. tropica complex. The L. killicki cluster showed low structuring and low polymorphism (see Figure 1). In contrast, L. tropica was highly polymorphic with strong structuring supported by high bootstrap values and some links with the country of origin, especially for strains from Kenya and Yemen. The larger and main clade was composed by all the Moroccan strains with the addition of other strains from other countries.

Discussion
Previous studies using a small number of strains and different molecular tools and analytic methods [9,12,31,38,40] included L. killicki in the L. tropica complex, except the study by Rioux and Lanotte [39] in which it was considered as a separate phylogenetic complex. The present study wanted to improve the knowledge on L. killicki phylogeny and its evolutionary history relative to L. tropica by using a larger sample of L. killicki strains from different countries. The phylogenetic analyses performed in this study confirm the position of this taxon within L. tropica in agreement with the previous biochemical and genetic findings. The close phylogenetic relationships between these taxa were also confirmed by the low number of polymorphic sites compared to those found between various Leishmania species [30,31]. The phylogenetic tree shows that L. killicki creates an independent group within L. tropica with high bootstrap value and no common haplotypes between them. Nevertheless, this taxon is included in the L. tropica complex and our data indicate that the species status of L. killicki is not justified. Furthermore, based on the L. tropica complex diversity and the multiple monophyletic branches in this complex, if L. killicki were to be considered as a species, the L. tropica complex would be composed of many species. Therefore, we suggest calling this taxon L. killicki (synonymous L. tropica) as it was previously done before for L. chagasi (synonymous L. infantum) [9,11,49,50]. Further epidemiological and clinical studies in the different countries where this taxon has been reported will say whether the L. killicki denomination should be maintained. From an evolutionary point of view, these data strongly suggest that L. killicki descends from L. tropica following only one founder event. This hypothesis is supported by the structure of the phylogenetic tree and by biochemical and genetic data. Indeed, the isoenzymatic characterization showed a low number of L. killicki zymodemes compared to those of L. tropica. This low polymorphism in L. killicki was confirmed by the low numbers of PS, HS and haplotypes and amino acid variations in the sequence of the different strains. The analysis of the phylogenetic tree suggests that L. killicki could have originated from an L. tropica ancestor from the Middle East. This ancestor would have separated into L. tropica in Morocco and other countries and into L. killicki in several other countries.
Finally, the lack of shared haplotypes and the identification of the new zymodeme MON-317 and its own haplotype suggest that L. killicki is now evolving independently from L. tropica, probably due to their different transmission cycles (zoonotic for L. killicki [51,52] and both anthroponotic and zoonotic for L. tropica [45,53,54]).
As the L. killicki strains showed low structuring and low polymorphism, we could not determine the precise evolutionary history of this taxon and particularly the country in which it emerged for the first time. Based on the epidemiological data, the higher genetic diversity and especially the relatively high number of described cases in Tunisia compared to the other countries [35,[55][56][57], it is likely that this taxon has emerged for the first time in Tunisia and then has spread in other North-African countries. Nevertheless, this should be further investigated.

Conclusion
The present work brings new insights into the evolutionary history of L. killicki and its taxonomic classification relative to L. tropica. However, more investigations need to be carried out on this model and particularly a detailed population genetics analysis to better understand the epidemiology and population dynamics of this parasite in comparison to L. tropica. revision of the manuscript; ALB participated in the analysis, interpretation of data and has contributed to the draft and revision of the manuscript; NH has been involved in the revision of the manuscript; PL and LT have participated in the technical experiments; FEB has contributed to data analysis; KJ and ZH have participated in sample collection; JPD was involved in the revision of the manuscript; HB directed the study; FP directed the study, revised and approved the manuscript. All authors read and approved the final manuscript.