Phylogenetic position of Leishmania isolates from Khyber Pakhtunkhwa province of Pakistan

Several species of the genus Leishmania are causative agents of cutaneous leishmaniasis in Pakistan. This study aimed to determine phylogenetic placement of Leishmania species causing cutaneous leishmaniasis in Khyber Pakhtunkhwa province, Pakistan (34 Leishmania tropica , 3 Leishmania infantum ), in-relation to species from other geographical areas using gene sequences encoding cytochrome b ( cytb ) and internal transcribed spacer 2 ( its2 ). Based on cytochrome b sequence analysis, L . tropica strains from Pakistan and other geographical regions were differentiated into two genotype groups, A and B. Within the province, ﬁ ve distinct L . tropica genotypes were recognized; two in group A, three in group B. Two L . infantum isolates from the province were closely associated with both Afro-Eurasian and American species of the Leishmania donovani complex, including Leishmania chagasi , L . infantum and L . donovani from Sudan and Ethiopia; while a third L . infantum isolate could not be differentiated from visceralizing Kenyan and Indian L . donovani . We observed apposite phylogenetic placement of CL-causing L . tropica and L . infantum from Khyber Pakhtunkhwa. Af ﬁ nities ascribed to Leishmania spp. From the region are valuable in tracing potential importation of leishmaniasis. © 2016 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (


Introduction
Leishmaniasis, caused by protozoan trypanosomatid parasites of genus Leishmania, affects skin (cutaneous), mucous membranes (mucocutaneous) or internal organs (visceral), primarily in humans. Recent figures suggest an approximate incidence of 0.2e0.4 million visceral leishmaniasis (VL) and 0.7 to 1.2 million cutaneous leishmaniasis (CL) cases per annum. Compared to VL, CL is more widely distributed; about one third of cases occur in each of three regions, the Americas, the Mediterranean basin and western Asia through Middle East to Central Asia (Alvar et al., 2012).
Taxonomy of Leishmania has been contentious since the genus was first described (Leishman, 1903). Lainson and Shaw (1987) classified Leishmania into sub-genera (L.) Leishmania, (L.) Viannia and (L.) Sauroleishmania on the basis of parasite development within the sandfly gut and by differences in extrinsic characters such as geographical distribution, vector preferences, clinical outcomes (VL, CL or MCL) and antigenic characteristics (Lainson and Shaw, 1987;Pratt and David, 1981). To date, multilocus enzyme electrophoresis (MLEE) remains the recognized gold standard for Leishmania species identification (Rioux et al., 1990). However, this traditional classification has been challenged using nuclear and mitochondrial molecular markers Mauricio et al., 2004Mauricio et al., , 2006Zemanova et al., 2004). There are a number of putative Leishmania species for whom species status is still under question including Leishmania panamensis, Leishmania peruviana and Leishmania lainsoni of sub-genus (L.) Viannia and L. garnhami, Leishmania pifanoi, Leishmania infantum, Leishmania chagasi, Leishmania killicki and Leishmania archibaldi of sub-genus (L.) Leishmania (Banuls et al., 2007).
Such phylogenetic studies are highly relevant in the context of disease control, since identifying discrete genetic groups can improve our understanding of virulence, transmissibility and drug susceptibility. In this study, cytochrome b (cytb) and internally transcribed region 2 (its2) were targeted to investigate the phylogenetic affinities of Leishmania tropica and L. infantum samples from Khyber Pakhtunkhwa (KP), Pakistan in-relation to species from other geographical areas.  (Khan et al., 2016 in prep). To mention briefly, localized lesions from these patients were sanitized and punctured with sterile lancets. Exudates collected using sterile pastettes were dispensed into biphasic culture medium containing rabbit blood agar and M199 medium (Sigma) supplemented with 10% heat inactivated fetal calf serum (HIFCS). DNA was extracted from pellets of positive cultures using commercial QIAGEN DNeasy Blood and Tissue Kit. For some patients DNA was also extracted from filter papers (coarse porosity) with impressions from lesions or their biopsies using resin-based Chelex ® method (Plowe et al., 1995). All samples (cultures and filter papers) were pre-identified to Leishmania species level (using KDNA primers) (Noyes et al., 1998) before they were subjected to cytb or its2 amplification.

Leishmania isolates and DNA extraction
In addition to the cultures (N ¼ 34) and filter papers (N ¼ 111) acquired from fieldwork in KP, DNA was also extracted from in vitro cultures of L. tropica (N ¼ 88) and Leishmania donovani (N ¼ 2) strains belonging to other geographical regions from the cryobank repository at the LSHTM (see Supplementary Table S1 for detail).

Cytb amplification and sequencing
Cytb, a maxicircle gene, was amplified using nested primers described by Luyo-Acero et al. (2004) Nest1 (N1) primers corresponded to COIII (cytochrome c oxidase subunit III) and MURF4 (maxi-circle unidentified reading frame 4). Nest2 (N2) primers LCBF1 and LCBR2 amplified approximately 866 bp of the internal cytb region. The cycling conditions for N1 and N2 included, an initial incubation for 1 min at 94 C, followed by 40 cycles of 94 C for 1 min, 50 C for 1 min, 72 C for 1 min and a final extension for 5 min at 74 C. DNA was extracted from gel sections using QIAquick Gel extraction Kit (QIAGEN). DNA sequencing was carried out using a set of internal primers LCBF4 and LCBR4 following manufacturer's protocol for BigDye™ terminator cycle sequencing kit V3.1.

its2 amplification and sequencing
The its2 region of chromosome number 27 was amplified (Bulle et al., 2002). Cycling conditions employed had an initial incubation for 2 min at 94 C, followed by 40 cycles of 94 C for 20 s, 53 C for 30 s, 72 C for 1 min and a final extension of 72 C for 10 min. Amplicons excised from gel were processed by QIAquick Gel extraction Kit (QIAGEN). Extracted DNA fragments were sequenced using BigDye™ sequencing kit V3.1. Three L. tropica strains (MRAT/ IQ/73/ADHANIS, MHOM/AZ/1974/SAF-K27 and MHOM/PK/1987/ FAZAL) were sequenced from clones derived using the CloneJET™ PCR cloning kit (Fermentas). We specifically undertook PCR cloning of multiple sequences per parasite isolate to circumvent the problems of differing lengths of variable regions.

Cytb and its2 phylogenetic analysis
Sequence chromatograms were analyzed and edited in Geneious v5.5.7 (Kearse et al., 2012). These sequences along with those retrieved from GenBank were used to derive Maximum likelihood (ML) trees and to perform, in parallel, Bayesian analysis. MEGA v5.05 (Tamura et al., 2011) was used to align sequences and to identify single nucleotide polymorphisms (SNPs). Trypanosoma cruzi and Herpetomonas muscarum were treated as out-groups in cytb and its2 analyses respectively. Alignments created from MEGA v5.05 were run in jMODELTEST v1.0 to select a best fit nucleotide substitution model (Posada, 2008).
For cytb the best model selected for the dataset was used to construct a Maximum Likelihood tree (ML) tree in MEGA V5.05. Bayesian analysis using the same model was carried out in TOPALi v2.5 (Milne et al., 2004) to produce posterior probabilities (abbreviated as PP) supporting the topology of ML bootstrap (abbreviated as BP) tree. For its2, both ML tree and Bayesian phylogenetic analysis was carried out in TOPALi v2.5 under same settings used for cytb. Cytb and its2 targets sequenced from this study have been submitted to GenBank (KT972143-KT972260, KT972262-KT972277).

Results and discussion
Using the diagnostic nested PCR of Noyes et al. (1998), all cultures acquired from KP were identified as L. tropica (N ¼ 34). Furthermore three L. infantum (N ¼ 3) were identified from filter paper spots, although they could not be successfully established in cultures.
Cytb genotypes from a total of 122 L. tropica cultures were analyzed along with an additional set of 41 sequences retrieved from the GenBank. No cytb product could be amplified from filter paper samples identified as L. infantum from Pakistan and thus could not be included in the analysis (Supplementary Table S1). Approximately 653 bp (with gaps) cytb gene alignment was analyzed.

Phylogenetic position of Pakistani L. tropica
Amongst all the L. tropica cytb sequences (122 from the study and 2 from Genbank) analyzed in the present study, 39 variable sites were identified in the sequence alignment giving a total of 21 genotypes. The co-occurrence of 15 single nucleotide polymorphisms (SNPs) within the L. tropica dataset discriminated the samples into 2 groups, A (12 genotypes among 101 strains) and B (9 genotypes among 22 strains). The most common genotype in the dataset was a genotype of group A (A1) observed among 91 of 124 L. tropica strains. Within KP province of Pakistan, 5 genotypes were observed (2 genotypes in group A and 3 genotypes in group B) (Supplementary Table S1). L. tropica species complex was supported by moderate bootstraps in the cytb ML tree (BP ¼ 60) (Fig. 1).
A moderately-supported L. tropica cluster (BP ¼ 95; PP ¼ 0.99), consisting of all its haplotypes, was observed in the its2 ML tree (Fig. 2). Within L. tropica, there were sub-clusters of clones with both haplotypes (H1 and H2) of its2 (including ones from the current study), although not supported by significant bootstraps (i.e. >50%). The single successfully sequenced L. tropica isolate from Northern Pakistan sampled in the present study (MHOM/PK/2010/ Fig. 1. Cytb maximum-likelihood consensus tree constructed from the Leishmania spp. typed in the current study and those retrieved from Genbank. Data for L. tropica complex samples typed in this study are expressed as genotypes (Supplementary Table S1). Strains with no species designations are expressed with standard strain codes only. Consensus tree was constructed from 1000 bootstrap replicates of cytb sequences along with Bayesian analysis (5 independent runs over 1 million generations sampling every 10 simulations with 25% burn-in period) using GTR þ G þ I (General time reversible model with gamma distributed rates and invariant sites) as the best-fit nucleotide substitution model. Only significant bootstraps (50%) and posterior probabilities (0.50) are shown in red. *Outgroup. CMH061) was placed in the sub-cluster of H2 haplotype. Overall, homogeneity was observed among the L. tropica strains from different geographical sites, including the genetically discrete Namibian strain IROS/NA/76/ROSSI-II (Schonian et al., 2001;Schwenkenbecher et al., 2006) (Fig. 2).
In this study both the phylogenetic analyses (cytb and its2) placed L. tropica as a discrete species complex within the genus, comprising L. tropica and Leishmania aethiopica as its two integrated species (Figs. 1 and 2). Using hsp70, Fraga et al. (2010), however, could not assign an individual species status to L. aethiopica. L. killicki from Tunisia, given full species rank by MLEE (Rioux et al., 1990), is indistinguishable from L. tropica by both phylogenetic topologies in the present study, a finding supported by others (Schonian et al., 2001;Schwenkenbecher et al., 2006). The cytb tree generated in the present study identified two main genotypes of L. tropica species, whereas its2 haplotypes were almost homogenous throughout the L. tropica cluster. Evidently, its2 may not be a discriminatory target for resolution of intra-specific diversity within L. tropica. This lack of congruence between our two genotypic phylogenies at species level for L. tropica may also be due to intrinsic differences in mutation rate between mitochondrial DNA as compared to nuclear DNA or evidence of historical genetic recombination (Bastos-Silveira et al., 2012;Chen et al., 2009;Gomez-Zurita and Vogler, 2003;Messenger et al., 2012;Tibayrenc, 1999).

Phylogenetic position of Pakistani L. infantum
In the its2 analysis two of the L. infantum strains from Pakistan (MHOM/PK/2010/CMH041 and MHOM/PK/2010/CMH069) were indistinctly placed within the L. donovani complex clade (BP ¼ 100; PP ¼ 1.00). A distinct sub-cluster in the main complex consisting of South Asian (India, Sri Lanka and Bhutan) and Kenyan L. donovani strains (BP ¼ 60; PP ¼ 0.97) also included a Pakistani L. infantum strain MHOM/PK/2010/CMH052 from the present study (Fig. 2). Several reports have found Indian and Kenyan L. donovani are Fig. 2. its2 maximum-likelihood consensus tree constructed from the Leishmania spp. typed in the current study and those retrieved from Genbank. Strains with no species designations are expressed with standard strain codes only. The consensus tree was constructed from 1000 bootstrap replicates of its2 sequences (656 bp). Consensus tree was constructed from 1000 bootstrap replicates of its2 sequences along with Bayesian analysis (5 independent runs over 1 million generations sampling every 10 simulations with 25% burn-in period) using SYM þ G (Symmetric model with gamma distributed rate variation among sites) as the best-fit nucleotide substitution model. Only significant bootstraps (50%) and posterior probabilities (0.50) are shown in red. *Outgroup.
genetically distinct from visceralizing agents in Sudan and Ethiopia (Ibrahim and Barker, 2001;Kuhls et al., 2005;Lewin et al., 2002;Lukes et al., 2007;Mauricio et al., 2001). Allocation of a presently genotyped Pakistani strain MHOM/PK/2010/CMH052 in this subcluster may either be a consequence of true similarity, or lack of discrimination between sequences analyzed. However, this does not seem to be the case with the two other Pakistani isolates typed as L. infantum (MHOM/PK/2010/CMH041 and MHOM/PK/2010/ CMH069), which were plausibly grouped with L. Infantum in the L. donovani complex cluster. Neither L. donovani nor CL infections caused by L. infantum have been previously reported in Pakistan. The three samples genotyped as L. infantum by KDNA from CL patients in the present study are the first reports of such cases (Khan et al., 2016 in prep).
L. chagasi (L. infantum from the Americas) could not be differentiated from the Mediterranean, Chinese and Pakistani L. infantum species. This is in agreement with previous studies based on MLEE, RAPD and RFLP analysis of nuclear targets (its, gp63, microsatellites and other protein coding genes) (Cupolillo et al., 1994;Kuhls et al., 2005;Lukes et al., 2007;Mauricio et al., 2001Mauricio et al., , 1999Mauricio et al., , 2000Ochsenreither et al., 2006). Evidently both cytb and its2 phylogenies produced herein supported a unitary status for the four defined species; L. donovani, L. archibaldi, L. infantum and L. chagasi in the L. donovani species complex (Figs. 1 and 2). Tibayrenc (2006) put forward a phylogenetic species concept that states, "A species is a monophyletic group composed of the smallest diagnosable cluster of individual organisms, within which there is a parental pattern of ancestry and descent". Based on the demarcations of this concept all the four named species in the L. donovani complex should be defined together as a singular species. Several studies support this notion of full species status (Croan et al., 1997;Lukes et al., 2007;Luyo-Acero et al., 2004;Mauricio et al., 2000;Tibayrenc, 2006).

Conclusion
Phylogenetic topologies based on cytb and its2 placed CLcausing L. tropica and L. infantum from KP, Pakistan confidently within their generally accepted respective species complexes; the latter is the first report of cutaneous disease by L. infantum in Pakistan. These observations highlight the need for further epidemiological study of parasites from the L. donovani complex in Pakistan and its neighboring countries, to further clarify the taxonomic status of strains causing cutaneous and visceral leishmaniasis in Pakistan.