Leishmaniasis is widespread in Tunisia and considered to be one of the major health problems for the human population [1]. In Tataouine, southeast Tunisia, Leishmania (L.) major and L. killicki (synonymous L. tropica ) [2] are currently the causative agents of cutaneous Leishmaniasis (CL) [3]. Leishmania killicki is responsible for the chronic CL observed in villages built on the flank of the Rocky Mountains [3-5]. Leishmania major, the causative agent of zoonotic CL, probably emerged as a human pathogen in this area as a result of recent establishment of a high density of susceptible human populations at the margins of villages [3].

In Tunisia, Phlebotomus (P.) sergenti and P. papatasi are the proven vectors of L. killicki (syn. L. tropica ) and L. major species, respectively [6-8]. Recent studies on sandflies showed that P. papatasi is the most abundant in the burrows of rodents such as Meriones, which are known reservoir hosts of L. major. Indoors, P. sergenti and P. papatasi , known to be endophilic were co-dominant [8]. Phlebotomus riouxi (or P. chabaudi ), which is suspected to be linked to L. killicki transmission [9], was the dominant species in the semi-natural rocky habitats associated with the presence of Ctenodactylus (C.) gundi, described as a potential reservoir [7-10]. Concerning the reproductive strategy, the current assumption, mainly based on population genetics studies, is that Leishmania alternates between three modes of reproduction: clonality, allogamy (interspecific recombination), and endogamy (intraspecific recombination), varying according to the species and the environment [11-18].

The occurrence of genetic exchange between different Leishmania species (corresponding to allogamy) is now largely accepted. Several cases of hybridization have been described in natural populations between both genetically close or distinct species (L. braziliensis and L. panamensis/guyanensis [19]; L. panamensis and L. braziliensis [20]; L. braziliensis and L. guyanensis [21]; L. major and L. arabica [22]; multiple hybrid genotypes of Leishmania (Viannia ) species [23]; L. infantum and L. major [24]; L. donovani and L. aethiopica [25]; L. donovani, L. major and L.infantum [26]; L.braziliensis and L.peruviana ) [27]. All these hybrids were evidenced by the detection of heterozygous profiles reflecting genetic exchange between different species. Additionally, some reports also suggested that intraspecific hybridization events (corresponding to endogamy) could also occur [28-31]. Several in vivo studies demonstrated the capacity of Leishmania species to undergo intraspecific genetic exchange in sandflies [32-34]. Although genetic exchanges have been experimentally and empirically demonstrated, questions remain about its frequency in natural Leishmania populations and their impact on the evolution of these populations.

In this context, our objective was to explore the existence of hybridization patterns between L. major and L. killicki by multi-locus microsatellite typing (MLMT) in a CL focus, where the two species are spreading. A total of 64 Leishmania isolates from patients with CL in the southeast of Tunisia (Tataouine governorate) were analyzed.

Study area

The study was carried out in the governorate of Tataouine, located in far southeast Tunisia at 32 55′40.16″ North, 10 26′57.28″ East and an altitude of C. 300 m above sea level [35] (Figure 1). Tataouine is characterized by its vast area (23.7% of the country) and by its international frontiers (Libya and Algeria). This region covers an area of 38.889 km² with a population of 149.453 (2014 census) [36]. Tataouine is subdivided into three biogeographical zones: the desert, the mountain (Dhaher), and the plain (Djefara) (Republic of Tunisia Ministry of Development, Investment and International Cooperation, South Development Office). It has a Mediterranean-desert climate [35]; the average temperature is 22°C and annual rainfall varies between 88 and 157 mm [37]. The regional landscape is characterized by poor vegetation, dominated by steppe species, and rocky escarpments. The area is pastoral with olive groves and cereals cultivated around water points. Ctenodactylus gundi is the most prevalent wild rodent living in the stony mountains. However, gerbils (Meriones species) are found at the periphery of villages [38].

Figure 1: Geographical localization of samples. A: Map of Tunisia, with Tataouine governorate highlighted in green; B: Districts of Tataouine governorate; C: Localization of the samples; villages are indicated with red asterisks. Pie charts represent Leishmania species frequencies for each locality.

Parasite samples

Analysis of isolates from 64 human CL cases from the Tataouine governorate was undertaken. The medical records of the patients were collected and provided by the local government hospitals. Data related to age, address, sex, disease onset, number and location of the lesions were recorded, respecting the anonymity of the patients. Cases were mainly recorded from Farech, Ghomrassen, and Guermessa localities. Direct smear and culture growth in modified NNN medium followed by sub-culture in RPMI-1640 was performed for each case during the diagnosis of the disease as part of the standard of care.

Leishmania identification

Some isolates were previously characterized by the multi-locus enzyme electrophoresis technique performed in the World Health Organization (WHO) reference center of Leishmaniasis in Montpellier “Centre National de Référence des Leishmanioses”; they are indicated in Supplementary Table 2 with a WHO code (unpublished data). For our study, genomic DNA was extracted from the 64 parasite cultures (promastigotes) using the QIAamp^® DNA Mini Kit (Qiagen, CA) according to the manufacturer’s protocol. DNA quality was checked by agarose gel electrophoresis and the concentration was measured with a Nano Drop spectrophotometer. For identification of the Leishmania species, a nested PCR, previously designed by Noyes et al. (1998) was performed with modifications [39,40].

The presence of a product of approximately 560-bp indicated L. major and a product of approximately 750-bp indicated L. tropica . Negative control were included in each reaction mix. Two reference strains, L. tropica MHOM/SU/74/SAF-K27 and L. major MHOM/SU/ 73/5-ASKH, were included for comparison.

Microsatellite markers

A panel of 11 microsatellite markers designed for the characterization of L. tropica species or L. major was selected for their potential ability to amplify the two species (Table 1). The selected microsatellite markers were previously described in Schwenkenbecher et al. [41]. As it was recently confirmed that L. killicki belongs to the L. tropica complex with close genetic relationships [2], we used the L. tropica microsatellite markers for genetic characterization of L. killicki . In addition to microsatellite markers designed for L. tropica with positive PCR products for L. major, we selected microsatellite loci originally designed for L. major and also able to be amplified in L. tropica .

Table 1: Details of the 11 microsatellite loci selected for the characterization of L. major and L. killicki species [28].

All the selected loci were tested on three L. major strains (one reference strain and two isolates from our sample) and three L. tropica and L. killicki strains (one reference strain and two isolates from our sample). PCR conditions were optimized for the 11 markers to increase the probability of amplifying the loci for both species. Approximately 50 ng of genomic DNA was added to each PCR. A negative control was included in each experiment. Microsatellites were amplified using the following conditions: 94°C for 2 min, followed by 35 cycles of 94°C for 30s, annealing temperature for 1 min (Table 1), and 72°C for 1 min (extension), followed by a final extension at 72°C for 10 min. Genotyping was performed by capillary electrophoresis using an ABI Prism 3130 XL automated DNA sequencer (Applied Biosystems, USA). Fragment size was determined automatically using GeneMapper 4.0 software (Applied Biosystems). Genescan 500 LIZ (Applied Biosystems, France) was used as internal size standard. The loci were then selected and evaluated according to the peak quality and sizing information.

Genetic and Phylogenetic Analyses

FSTAT Version 2.9.3.2 software [42] updated from [43] was used to analyze the data. Genetic polymorphism was measured based on the allelic richness (A) and the Nei’s unbiased estimate of genetic diversity within subsamples (Hs) [44]. FST measures the relative inbreeding in subpopulations attributable to the subdivision of the total population into subpopulations of limited size. FST also measures the geneti differentiation between subpopulations. The significant departure from zero of these parameters was tested by 10.000 randomization procedures with FSTAT. The significance of these estimates was confirmed by p-values ≤ 0.05. A neighbor-joining tree, based on Nei’s minimum genetic distances was used to cluster the genotypes. Data were computed using POPULATION software to build the distance matrix (version 1.2.28; Centre National de la Recherche Scientifique, UPR9034, Langella, O) and the phylogenetic tree was generated using MEGA software, version 4.0.2 [45].

Lesion type and clinical impression

The lesions were unique in 68.75% (n=44) of cases and showed a clinical polymorphism. The most frequent pattern of eruption was seen in 78.12% of cases (n=50) patients who had a nodular ulcerative crusting aspects localized predominantly in the limbs (65.62%), followed by nine patients who had plaque-like ulcerative eruption followed by 5 patients who had a dry plaque/nodular eruption. All dry lesions were found to be Leishmania killicki (syn. L. tropica ) (see supplementary Tables 1,2 and Figure 2).

Figure 2: Clinical spectrum of the cutaneous Leishmaniases. The disease encompasses a range of clinical symptoms, including infiltrated plaque, ulcero-crusted lesions of the leg. (A: lesion caused by L. major) , small Ulcerative nodule of the hand (B: Lesion caused by L. killicki ), and non-ulcerative nodule (C: Lesion caused by L. killicki ), ulcero-crusted nodule of the ear (D: Lesion caused by L. killicki ) and ulcerative plaque on the face (E: Lesion caused by L. killicki ).

Identification of Leishmania species by nested PCR

All isolates analyzed presented a single Leishmania kinetoplast minicircle amplicon of approximately 560-bp or 750-bp, corresponding to L. major and L. tropica species, respectively. These identifications corroborated the CNRL data for L. killicki and L. major strains of the 64 samples under study, 51 isolates showed an amplicon size corresponding to L. tropica and 13 to L. major. The identification of each sample is detailed in Supplementary Table 2 and the village distribution of the two species in Supplementary Table 1 and Figure 1. Both chronic and zoonotic CL co-exist in the Tataouine governorate with L. killicki (syn. L.tropica ) as the predominant species. The majority of the villages under study contained the two species except Tlelet (seven L. killicki strains) and Oued Khil (one L. killicki strain) localities (see Figure 1).

Microsatellite genotyping

From the 11 microsatellite loci selected from Schwenkenbecher et al. [41], our choice was essentially based on the ability of primer pairs to amplify products from the two species. Eight markers were usable for our purpose including those originally designed for L. major (4GTG and 27GTG) [46] (see Table 1). Three markers (GA2, GA6, and GACA4) could not be used because of their inability to amplify both species (GA6), imperfect reproducibility (GACA4), or, in the case of GA2, under the same MLMT-PCR conditions; primer pairs amplified a 66-bp product for L. killicki and a 499-bp product for L. major , suggesting different genomic targets in the different species.

Forty of the 64 isolates yielded complete genotype data for the eight selected loci. In some cases, no PCR products were obtained in repeated PCR runs, and were treated as missing data for the statistical analyses (see supplementary Table 2). The eight microsatellite markers were polymorphic, based on number of microsatellite repeats, across the whole sample set, including the reference strains, and five of them revealed heterozygous patterns. Among the 51 L. killicki and 13 L. major isolates genotyped with the eight loci, a total of 20 different alleles were identified and allelic richness (A) was 1.12 and 1.25, respectively.

The number of alleles per locus ranged from two (GA1, GA3, GTG3, and 4GTG) to three (GA10, Mix9, GT4, and 27GTG). The level of unbiased gene diversity (Hs) across populations was very low for the two species (Hs=0.080 for L. major ; Hs=0.012 for L. killicki ) (Table 2). Overall, 20 different genotypes were detected among the 64 samples analyzed; 5 genotypes fo L. major populations of 13 isolates and 15 genotypes for L. killicki populations of 51 isolates. Nine of the 13 L.major isolates from Mdhila, Ferech, and Guermessa displayed an identical heterozygous genotype (64/68) at locus Mix9 (see Supplementary Table 2 and samples with green dots on Figure 2). Four isolates from Ghomrassen and Ksar Mourabtin differed from the first nine either by a distinct allele at locus 27GTG (GHO010) or were homozygous at the Mix9 locus (GHO022, KSM001, and KSM006). Among the L. killicki isolates, we obtained two clusters of 28 and 9 isolates showing identical genotypes in each cluster, of which 11 samples had missing data at one or more markers. Seventeen strains had unique genotypes, of which five were heterozygous at loci GA10 (79/85) or GT4 (72/78) (see Supplementary Table 2 and Figure 2). The L. tropica reference strains (SAF-K27) revealed four different heterozygous patterns at loci GA10, GTG3, GT4, and 27GTG.

Homozygous genotypes prevailed in our samples. As described above, heterozygous profiles were recorded for 14 isolates (without including L.tropica reference strain, MHOM/SU/74/SAF-K27) and on only one locus for each (MIX9 or GT4 or GA10). Ho ranged from 0 to 0.818 and he ranged from 0 to 0.506 per locus in the L. major population. For L. killicki , observed and expected values ranged from 0 to 0.063 and from 0 to 0.062, respectively (data not shown). The mean observed heterozygosity was very low for the L. killicki population and slightly higher for the L. major population (Table 2). Mean expected heterozygosity within the L. killicki population was similar to the mean observed heterozygosity and different to that of the L. major population (see Table 2)

Table 2: Genetic diversity indices estimated from microsatellite data (eight loci) for the 64 L.major and L . killicki isolates included in this study.

Furthermore, it is worth noting that all the heterozygotes displayed one allele observed in the homozygous state in the same species and a second allele absent in the homozygous state. For example, the heterozygotes observed for L. major isolates at the locus Mix9 displayed allele (68) found in the homozygous state in the L. major population (68/68) and allele (64), which is absent in the homozygous state in the whole population (L. major + L. killicki ) (see Supplementary Table 2). The same pattern is observed in the L. killicki population for the GT4 and GA10 markers. The mean FST value estimated was 0.966 (p-value ≤ 0.05), showing, as expected, a significant level of differentiation between the two species [47]. The data from the MLMT analysis on the 64 L. major and L. killicki isolates and the two reference strains were used to create matrix distances and construct a neighbor-joining tree representing the relationships among all the strains used in this study. As expected, the MLMT data classified isolates according to the species taxonomy (see supplementary Table 2). Strains with heterozygous genotypes (tagged with green dots) clustered either in L. killicki cluster or in L. major

one in agreement with their species identification. The phenetic tree did not reveal any cluster with an intermediate position between the two species (see Figure 3).

Figure 3: Clinical spectrum of the cutaneous Leishmaniases. The disease encompasses a range of clinical symptoms, including infiltrated plaque, ulcero-crusted lesions of the leg. (A: lesion caused by L. major ), small Ulcerative nodule of the hand (B: Lesion caused by L. killicki ), and non-ulcerative nodule (C: Lesion caused by L. killicki ), ulcero-crusted nodule of the ear (D: Lesion caused by L. killicki ) and ulcerative plaque on the face (E: Lesion caused by L. killicki ).

L. killicki was first detected in the micro-focus of Tataouine in southeast Tunisia in 1980, in an epidemic of CL [4,48]. The first cases of L. major in Tataouine were only recorded in 1991 [49,50], but now chronic and zoonotic CL co-exist in this governorate [3-5]. In agreement with previous epidemiological reports, our analysis found that L. major and L. killicki (syn. L.tropica ) are the only species circulating in the Tataouine governorate. The isolates were collected in different villages of the investigated region, confirming the geographical spread of Leishmania [3,7,9]. Contrary to previous studies, the most abundant species was L. killicki , accounting for 80% of our sampling. This difference between our study and previous studies was probably due to the investigation of different villages [3].

The two species showed a very low level of genetic diversity with this set of microsatellite markers. The low degree of polymorphism of L. killicki is in agreement with recent data published by Chaara et al. [2]. The authors suggested that the recent emergence and the low level of transmission in humans of this taxon can explain the low level of genetic diversity. Nevertheless, conversely to the data of Chaara et al. [2], we obtained similar observed and expected heterozygosity. This difference can be explained by the extremely low genetic diversity (Hs=0.012) in our study precluding the analysis of the L. killicki population structure. For L. major , the previous published data showed a higher genetic diversity in Gafsa, Kairouan, and Sidi Bouzid governorates (Hs=0.321, 0.306, and 0.236, respectively) compared with the results obtained in our study (Hs=0.080) [51], because L. major would have emerged simultaneously in Tataouine and Sidi Bouzid, the low heterogeneity found in our analysis is more likely due to the low number of isolates or the different sets of microsatellites between the two studies. We could not further investigate the population structure of this species in Tataouine because of the reduced sample size and low level of genetic diversity, Despite the low level of genetic diversity in each species, the marked differentiation between L.major and L.killicki allowed us to explore patterns of hybridization between the two specie. Furthermore, since the microsatellites are co-dominant and generally neutral markers [52], they have been widely used to study genetic exchange within many organisms [53-55] of which Leishmania at both intra and interspecific level (between L. (Viannia) braziliensis and L. (Viannia ) guyanensis [17]; between L.braziliensis and L. peruviana [23]; between L.guyanensis and L.panamensis [56]; between MON-1 and MON-24-80 of L. infantum [29,57]; between MON-1 and non MON-1 of L.infantum [58]; within L.tropica [59,60]; within L. donovani [31]. In our sample, neither the raw data nor the phenetic analysis suggests the occurrence of interspecific genetic exchange. Fourteenisolates (without including L.tropica MHOM/SU/74/SAF-K27 reference strain), belonging either to L.killicki or to L.major , revealed heterozygous genotypes, with no shared allele between the two species. The accumulation of mutations, consisting of gain/loss of repeat units in microsatellites [61,62] may explain these heterozygous profiles. Both species (L. major and L. killicki ) were collected in Ghomrassen, Guermessa, Farech, and Mdhila villages. The absence of interspecific genetic exchange might be explained by different epidemiological cycles at the sympatric level. These two species seem to have overlapped geographical foci but different vector species and different mammalian reservoir hosts [7,9,10,63], limiting the interactions between them. As reported by Tabbabi et al. the incriminated vectors for L. major and L. killicki in the Tataouine region would be species specific [9]. Thus, the most likely opportunity for these two species to interact is during mixed infections in humans. According to published data, cross-species genetic exchange between Leishmania strains would mainly occur in the invertebrate stage [32-34,64]. These epidemiological and biological aspects could explain the absence of interspecific genetic exchange between L. killicki and L. major in the Tataouine governorate. Nevertheless, from our data we cannot infer any hypothesis on intraspecific genetic exchange.

In our study, we detected two sympatric species responsible for CL in southeast Tunisia, Tataouine governorate: L. major and L. killicki . Each species revealed very low intraspecific diversity but high genetic differentiation. Nine L. major strains and five L. killicki strains displayed heterozygous profiles probably owing to genetic mutations and not to genetic exchange. Specific and different epidemiological cycles at the sympatric level might explain the absence of cross-species genetic exchange between the two species in the Tataouine governorate.

The authors thank the Regional Health Directorate of the Tataouine governorate (Tunisia), the “Institut Pasteur Tunis” and the “Institut de Recherche pour le Développement” (Montpellier, France), and the “Centre National de la Recherche Scientifique” (France) for this cooperative and collaborative project. We acknowledge the field staff of the Department o Medical Biology, Institut Pasteur Tunis, for strain collection and cryobank management; Adel Gharbi and Mongi Dellagi, for ensuring strain and patient data collection in accordance with current Good Clinical Practice guidelines. We gratefully acknowledge the “Centre National de Référence des Leishmanioses” for the multilocus enzyme electrophoresis typing of some strains.

The authors declare that there is no conflict of interest regarding the publication of this paper.