INTRODUCTION
The Indian subcontinent is the most affected region for visceral leishmaniasis (VL) in the world, accounting for 60% of the global disease burden. An estimated 300 000 new cases occur annually, and 190 million people are at risk of infection (Sundar et al. Reference Sundar, Mondal, Rijal, Bhattacharya, Ghalib, Kroeger, Boelaert, Desjeux, Richter-Airijoki and Harms2008). The disease is caused by Leishmania donovani, which is transmitted by the bite of infected female sand flies. The disease epidemiology is changing rapidly due to Leishmania-HIV co-infections, the high rate of treatment failure, and rapid expansion to previously non-endemic regions (Thakur et al. Reference Thakur, Meenakshi Thakur and Thakur2009; S. Rijal, unpublished observations).
In the context of disease burden, the governments of Bangladesh, India and Nepal signed an agreement in 2005 to eliminate the disease from the region through early diagnosis and complete treatment, vector control, social mobilization of the population at risk and building partnership, networking and operational research (Joshi et al. Reference Joshi, Narain, Prasittisuk, Bhatia, Hashim, Jorge, Banjara and Kroeger2008). However, epidemiological data such as circulating strains, vector infection and spread, possible non-human reservoirs, treatment failure and drug resistance are not systematically documented in this region. Adequate tools for monitoring the (micro-) epidemiology of the disease are needed to support the elimination programme with regard to the follow-up of new epidemics, and the spread of drug-resistant clones related to treatment failure.
So far, multilocus enzyme electrophoresis (MLEE) is still being used as the reference method for Leishmania typing (Rioux et al. Reference Rioux, Lanotte, Serres, Pratlong, Bastien and Perieres1990). However, this is a time-consuming and laborious technique since it requires mass cultivation of parasites and specialized equipment. In addition, variations due to isoenzymes reflect only neutral mutations affecting electrophoretic mobility, leading to a lack of discriminatory power to fully understand limited parasite diversity.
Different PCR-based techniques with improved discriminatory power have been described for Leishmania typing (Botilde et al. Reference Botilde, Laurent, Quispe Tintaya, Chicharro, Cañavate, Cruz, Kuhls, Schönian and Dujardin2006; Reithinger and Dujardin, Reference Reithinger and Dujardin2007). These methods are easier to use, allow the study of the genetic diversity of parasite populations at different levels, and can often directly be applied on host tissue, thereby overcoming time-consuming mass cultivation and reducing the risk of selective growth in culture. The level of polymorphism and discriminatory power of PCR-based typing methods at the specific and intra-specific level essentially depends upon the genetic locus targeted as molecular marker, and on the methodology for exploring it. Polymorphisms in PCR amplicons are often revealed on the basis of their size, examples of which are microsatellites and minisatellites (Jamjoom et al. Reference Jamjoom, Ashford, Bates, Kemp and Noyes2002; Bulle et al. Reference Bulle, Million, Bart, Gállego, Gambarelli, Portús, Schnur, Jaffe, Fernandez-Barredo, Alunda and Piarroux2002; Alam et al. Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009; Ochsenreither et al. Reference Ochsenreither, Kuhls, Schaar, Presber and Schönian2006; Schwenkenbecher et al. Reference Schwenkenbecher, Wirth, Schnur, Jaffe, Schallig, Al-Jawabreh, Hamarsheh, Azmi, Pratlong and Schönian2006; Kuhls et al. Reference Kuhls, Keilonat, Ochsenreither, Schaar, Schweynoch, Presber and Schönian2007; Haralambous et al. Reference Haralambous, Antoniou, Pratlong, Dedet and Soteriadou2008; Seridi et al. Reference Seridi, Amro, Kuhls, Belkaid, Zidane, Al-Jawaberh and Schönian2008). Other methods discriminate not only by size, but also by the presence of restriction endonuclease recognition sites using a technique known as RFLP (restriction fragment length polymorphism). Examples of targets where this was applied are gp63 (glycoprotein 63 gene) (Victoir et al. Reference Victoir, Bañuls, Arévalo, Llanos-Cuentas, Hamers, Noël, De Doncker, Le Ray, Tibayrenc and Dujardin1998; Guerbouj et al. Reference Guerbouj, Victoir, Guizani, Seridi, Nuwayri-Salti, Belkaid, Ben Ismail, Le Ray and Dujardin2001; Mauricio et al. Reference Mauricio, Gaunt, Stothard and Miles2001), kinetoplast minicircle DNA (kDNA) (Morales et al. Reference Morales, Chicharro, Ares, Cañavate, Barker and Alvar2001; Cortes et al. Reference Cortes, Mauricio, Almeida, Cristovão, Pratlong, Dedet and Campino2006; Laurent et al. Reference Laurent, Rijal, Yardley, Croft, De Doncker, Decuypere, Khanal, Singh, Schönian, Kuhls, Chappuis and Dujardin2007; Nasereddin et al. Reference Nasereddin, Azmi, Jaffe, Eregat, Amro, Sawalhah, Baneth, Schönian and Abdeen2009), cpb (cysteine proteinase B gene) (Quispe Tintaya et al. Reference Quispe Tintaya, Ying, Dedet, Rijal, De Bolle and Dujardin2004), the mini-exon (Marfurt et al. Reference Marfurt, Nasereddin, Niederwieser, Jaffe, Beck and Felger2003; Mauricio et al. Reference Mauricio, Stothard and Miles2004; Pandey et al. Reference Pandey, Yanagi, Pandey, Mallik, Sherchand and Kanbara2007), and rDNA internal transcribed spacers (Mauricio et al. Reference Mauricio, Stothard and Miles2004). All these studies demonstrated the added value of PCR-based techniques to document the structure of a Leishmania population in a micro-epidemiological context.
Most reports so far have documented L. donovani populations in the Indian subcontinent to be genetically homogeneous (Alam et al. Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009), even though the structure of the parasite population has not been extensively tested with different PCR-based molecular tools in a single study. Here, we explored the discrimination performance of 6 strategies for documenting the heterogeneity of L. donovani from Nepal, namely gp63-intragenic PCR-RFLP, cpb-intragenic PCR-RFLP, cpb-intergenic PCR-RFLP, size polymorphism in haspB, multilocus microsatellite typing (MLMT), and kDNA minicircle PCR-RFLP. For kDNA minicircle PCR-RFLP, standard experimental procedures were developed to maximize repeatability and allow inter-experiment comparisons.
MATERIALS AND METHODS
Samples
Thirty-four clinical isolates (Table 1) included in this study were obtained from bone marrow aspiration of clinical VL patients who were recruited from March 2002 until September 2003 at the B. P. Koirala Institute of Health Sciences (BPKIHS, Dharan, Nepal), a tertiary medical care centre at the eastern Terai of Nepal. Informed consent was obtained from patients or guardians, and ethical clearance was granted by the institutional review boards in Belgium and Nepal. In total, the clinical isolates originated from 7 districts (Fig. 1): 6 endemic districts with a high rate of VL transmission during the past, and the non-endemic hilly district Bhojpur from which a number of cases were reported. Briefly, bone marrow aspirates were directly inoculated into Tobie's blood agar medium at ±26°C for the cultivation of parasites (Tobie et al. Reference Tobie, Von Brand and Mehlman1950). Cultured promastigotes were harvested by centrifugation, and washed 3 times using phosphate-buffered saline, pH 7·2. Total DNA was extracted with the QIAamp DNA mini kit (Qiagen, Hilden, Germany) as per the manufacturer's instructions, and the DNA concentration was estimated with the Nanodrop ND-1000 UV-Vis spectrophotometer (NanoDrop Technologies, Wilmington, USA).
Geographical distribution of nuclear (Nx) and kDNA (Kx) genotypes of the 34 studied clinical isolates according to Table 1. Villages within the districts are indicated by dots, and for each village the recovered genotypes are listed.
Origin and genotypes of the 34 isolates included in this study
a WHO code, ordered by BPK number.
b gp63 intragenic pattern 1 is indicated in Fig. 2A.
c cpb intergenic pattern 1 is indicated in Fig. 2B.
d cpb intragenic patterns are indicated by group A and B and are shown in Fig. 2C.
e haspB amplicon size is indicated by group A and B and is shown in Fig. 2D.
f Microsatellite types are detailed in Table 2.
g nDNA genotypes are derived from combining the results of all nuclear markers, and are indicated as N1 to N8.
h kDNA genotypes are as in Fig. 3.
i Combined genotypes are derived from combining the results of kDNA and nDNA genotypes.
In order to check the inter-experiment repeatability of kDNA minicircle PCR-RFLP procedures, several biological and technical replicas were used (detailed in Fig. 3). DNA from the cultured isolates and spiked blood samples were extracted according to the method described above, except for some clones maintained in modified Eagle's medium (Decuypere et al. Reference Decuypere, Vanaerschot, Rijal, Yardley, Maes, De Doncker, Chappuis and Dujardin2008) for which TELT buffer (containing TRIS, EDTA, LiCl and Triton X-100) followed by phenol-chloroform extraction was used. All samples in the study were L. donovani, which was confirmed by amplifying a partial cpbF coding sequence of L. donovani with primers cpbF2.1 and Do2.1 (Laurent et al. Reference Laurent, Van der Auwera, Hide, Mertens, Quispe-Tintaya, Deborggraeve, De Doncker, Leclipteux, Bañuls, Büscher and Dujardin2009).
gp63-intragenic genotyping
PCR of the gp63 intragenic region was done with primers SG1 and SG2 as described elsewhere (Guerbouj et al. Reference Guerbouj, Victoir, Guizani, Seridi, Nuwayri-Salti, Belkaid, Ben Ismail, Le Ray and Dujardin2001). PCR reactions were carried out in 50 μl consisting of 1 ng template DNA, 200 μm of each deoxynucleoside triphosphate (Eurogentec, Seraing, Belgium), 1× reaction buffer (Eurogentec), 1 mm of MgCl2 (Eurogentec), 2·5 units of SilverStar DNA polymerase (Eurogentec) and 400 nm of each primer (Sigma-Aldrich, Bornem, Belgium). PCR was carried out in an Eppendorf Mastercycler EP Gradient S with cycling conditions as described elsewhere (Guerbouj et al. Reference Guerbouj, Victoir, Guizani, Seridi, Nuwayri-Salti, Belkaid, Ben Ismail, Le Ray and Dujardin2001). PCR amplicons of the correct size were precipitated in 0·3 m NaAc/70% EtOH, and resuspended in 10 μl of water. About 3 μl was digested with HincII (New England Biolabs, Ipswich, MA, USA) as per the manufacturer's instructions, using 15 units of enzyme at 37°C overnight. The restriction fragments were separated in a 3% small fragment agarose gel (Eurogentec) for 5 h at 4·1 V/cm, and visualized with ethidium bromide.
cpb-intergenic and cpb-intragenic genotyping
PCRs to amplify an intra- and inter-genic region of cpb were performed with primer combinations CPBFOR-CPBREV and PIGS1A-PIGS2B respectively (Quispe Tintaya et al. Reference Quispe Tintaya, Ying, Dedet, Rijal, De Bolle and Dujardin2004). Reactions were carried out in an Eppendorf Mastercycler EP Gradient S thermocycler, using 2·5 units of SilverStar DNA polymerase in the accompanying buffer (Eurogentec), and with conditions as described by Quispe Tintaya et al. (Reference Quispe Tintaya, Ying, Dedet, Rijal, De Bolle and Dujardin2004). Specific amplicons were precipitated as described for gp63 above, and 3 μl was digested with 15 units of HaeIII (MBI Fermentas, St Leon-Rot, Germany) at 37°C overnight, following the manufacturer's instructions. The cpb-intergenic amplicons were analysed on a 3% small fragment agarose gel as for gp63, whereas the cpb-intragenic amplicons were analysed by capillary electrophoresis using the DNA 1000 lab chip kit in the Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) as described elsewhere (Quispe Tintaya et al. Reference Quispe Tintaya, Ying, Dedet, Rijal, De Bolle and Dujardin2004).
haspB genotyping
The haspB gene was amplified with primers K26F and K26R as described elsewhere (Haralambous et al. Reference Haralambous, Antoniou, Pratlong, Dedet and Soteriadou2008), except that HotStarTaq Plus DNA polymerase (Qiagen) was used in its accompanying buffer. The size of the PCR amplicons was determined on a 2% agarose gel, and some were sequenced.
MLMT genotyping
Multilocus microsatellite typing was performed with the 5 microsatellite markers Li23-41, Li22-35, LM4TA, Li45-24, and Li41-56 (Alam et al. Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009). Briefly, all the microsatellite loci were amplified in 50 μl final reaction volume containing 1 ng of template DNA, 2 mM MgCl2, 200 μM of each deoxynucleoside triphosphate, 80 nm of both the unlabelled and 5′ fluorescently labelled primer (fluorophore 6-FAM for Li23-41, VIC for Li22-35 and Li45-24, NED for LM4TA, and PET for Li41-56, Applied Biosystems, Foster City, CA, USA), and 2·5 units of HotStarTaq Plus DNA polymerase in 1× PCR buffer (Qiagen). Amplification was done in an Eppendorf Mastercycler EP gradient S as follows: initial denaturation at 95°C for 5 min; subsequently 40 cycles, each consisting of a denaturation step at 95°C for 30 sec, annealing for 30 sec at 52°C (Li23-41, Li22-35 and Li41-56) or 54°C (LM4TA, Li45-24), extension at 72°C for 1 min; and finally an extension step of 10 min at 72°C. PCR amplicon size was determined on an ABI3730XL DNA sequencer with the GENEMAPPER software and internal size standard GS500(−250)LIZ (Applied Biosystems). Nineteen Nepalese L. donovani strains that were included in the study of Alam et al. (Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009) were used to calculate the repeat numbers.
kDNA minicircle genotyping
Primers BPKDNAMINFOR (5′-CTGGGGGTTGGTGTAAAATAGGGC-3′) and BPKDNAMINREV (5′-CCCGATTTTTGGCATTTTTGG-3′) were designed on the basis of an alignment of L. donovani kDNA minicircle sequences from the Indian subcontinent retrieved from GenBank (Accession numbers AF167712, AF167713, AF167714, AF167715, AJ010084, AJ010085, AJ010086, AJ010087, X68026, X68027, Y11402, AF399822). They were used for PCR amplification in a 50 μl solution containing 1× Coralload PCR buffer (Qiagen), a total of 2 mM MgCl2, 200 μM of each deoxynucleoside triphosphate, 500 nM of each primer, 2·5 units of HotStarTaq Plus DNA polymerase, and about 1 ng of template DNA. Cycling was performed in an Eppendorf Mastercycler EP gradient S at 95°C for 5 min, followed by 40 cycles each consisting of 1 min at 94°C, 1 min at 60°C, and 1 min at 72°C, after which a final extension step of 10 min at 72°C was performed. The PCR amplicons of approximately 800 nucleotides were precipitated as described above, after which the concentration of the precipitated amplicons was estimated on agarose gel by comparison to known amounts of the Gene Ruler 100 bp plus DNA marker (MBI Fermentas). Restriction digest was done on 150 ng of precipitated PCR amplicons in a total of 20 μl of the recommended buffer with 10 units HaeIII, incubated at 37°C for 2 h, after which the reaction was stopped by adding 1 μl of 0·5 m EDTA. The fragments were analysed after electrophoresis for 5 h at 4·1 V/cm in a 3% metaphor agarose gel (Lonza, Rockland, USA), which was stained with ethidium bromide.
The obtained fingerprint patterns were analysed using the BIONUMERICS software version 5.10 (Applied Maths, Ghent, Belgium). Briefly, spectral analysis was used to correct for background signal using a gel picture without signal saturation. Smiling and other gel artifacts were filtered out using the GeneRuler 100 bp DNA ladder run at both sides and in the middle of each gel. Subsequently, RFLP fingerprint patterns were compared on a two-by-two basis using the Pearson product-moment correlation coefficient (Pearson, Reference Pearson1926), calculating overall similarities between densitometric curves. Only the region between 200 and 600 bp was taken into account. Based on the obtained similarity values, a UPGMA dendrogram was constructed. The reliability of each branch was assessed using standard deviation based upon the correlation between the dendrogram-derived and the original curve similarities. Even though UPGMA would not be the method of choice when assessing relationships between parasite isolates, our aim was merely to pragmatically group isolates based on fingerprint patterns, regardless of their origin.
The above described protocol was obtained after extensive evaluation and optimization. Initially 6 restriction enzymes (AciI, AseI, HPYCH4V, MseI, MnII, and HaeIII) were selected on the basis of in silico analysis of the kDNA alignment, and were tested on a selected set of samples using visual fingerprint evaluation. As HaeIII was retained for the final protocol, conditions for complete digestion were evaluated by incubating PCR products from 0 min to overnight, after which the obtained RFLP patterns were compared visually. The sensitivity of the kDNA PCR was tested on serial dilutions of DNA isolated from parasite cultures and from blood spiked with promastigotes.
Genetic diversity
To determine the diversity of the molecular markers, the expected heterozygosity (He) was calculated by using the equation He=[n/(n−1)] (1−∑p i2), where n is the number of samples analysed and p i is the frequency of genotype i for each marker (Koepfli et al. Reference Koepfli, Mueller, Marfurt, Goroti, Sie, Oa, Genton, Beck and Felger2009).
RESULTS
Nuclear marker genotypes
The patterns obtained with the different fingerprint methods are depicted in Fig. 2 and listed in Tables 1 and 2, and their geographical distribution is shown in Fig. 1. Gp63-intragenic PCR-RFLP and cpb-intergenic PCR-RFLP patterns were monomorphic (Fig. 2A and B) and showed the typical L. donovani pattern encountered in reference strains from India (Quispe Tintaya et al. Reference Quispe Tintaya, Ying, Dedet, Rijal, De Bolle and Dujardin2004). On the other hand, cpb-intragenic PCR-RFLP and length polymorphism of haspB showed 2 types of genetic profiles (Fig. 2C and D). Out of the 34 clinical isolates, 32 (further referred to as group A) were identical for all 4 above-mentioned markers, while the 2 remaining ones (group B) were different in the cpb-intragenic and haspB assays. Sequencing showed an haspB amplicon of 320 bp in group B, as opposed to 641 bp in group A. No sequence diversity was seen within the 2 groups.
PCR and RFLP fragments from different markers as listed in Table 1, relative to size markers M1 (GeneRuler 100 bp Plus DNA ladder) or M2 (DNA 1000 lab chip kit). (A) gp63-intragenic PCR-RFLP pattern 1 (this was the only pattern encountered). (B) cpb-intergenic PCR-RFLP pattern 1 (this was the only pattern encountered). (C) Patterns obtained from cpb-intragenic PCR-RFLP of groups A and B. (D) haspB amplicons characteristic for group A (641 bp) and B (320 bp).
The repeat number of 5 analysed microsatellite loci is shown in Table 2. While locus Li23-41 showed 7 different alleles, the other 4 loci were monomorphic within groups A and B, but different between these groups. A combination of the 5 loci uncovers 8 different profiles, 5 of which (M2, M3, M4, M7, M8) are unique, while M1 is found in 23 isolates, M6 in 4, and M5 in 2. The 2 group B isolates have a unique microsatellite profile. Of all microsatellite profiles obtained, only locus Li23-41 shows a heterozygous pattern in some strains, i.e. genotypes M5 and M7. As depicted in Table 2, all obtained profiles could be matched to genotypes reported by Alam et al. (Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009), except for M2, M3, and M7.
Microsatellite genotypes of the 34 isolates included in this study
a WHO code, ordered by BPK number.
b Each genotype represents a unique combination of alleles. Allele combinations matching to profiles in Alam et al. (Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009) are given in parentheses, whereby M1 is also compatible with Li.
c Isolates included in Alam et al. (Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009).
d The 31 repeat allele was not found in Alam et al. (Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009).
Optimization of kDNA minicircle genotyping
kDNA PCR was able to amplify a single parasite spiked in 180 μl of blood, and 0·1 pg of parasite DNA, corresponding to half a Leishmania genome. RFLP patterns obtained from 1000 pg down to 1 pg of DNA were identical, as were those obtained from spiked blood samples containing 10 000 down to 10 parasites (not shown). Out of 6 restriction enzymes, HaeIII generated the highest degree of polymorphism (not shown). Digest times of less than 2 h up to overnight yielded identical fragments, from which it was decided to use 2 h.
kDNA genotypes
kDNA restriction patterns obtained from the 34 clinical isolates were used for cluster analysis (Fig. 3). Replicas from some isolates, as well as derived parasite clones, were included in separate experiments to test the robustness of the obtained clusters. One of these clones (BPK190/0-cl3) was kept in culture for 1 year, as kDNA minicircle populations have been known to change over time. The groups diverging at the left-hand side of the dotted line could be reliably obtained (i.e. in inter-experiment replicates the same isolate is consistently linked to the same cluster), while splits occurring at the right-hand side can often be attributed to inter-experiment variations. Based on the standard deviation of the branches, and consistency of clustering, 14 kDNA fingerprint types were recognized among the 34 clinical isolates (Table 1 and Fig. 3). The 2 group B isolates branch off at the base of the tree, separate from the rest.
UPGMA dendrogram based on fingerprint densitometric curve similarities from kDNA PCR-RFLP with HaeIII, applied on the 34 clinical isolates from this study, and biological and experimental replicas. Shaded horizontal lines indicate the standard deviation of each branching position as deduced from comparing dendrogram and curve similarities. The dotted line denotes the inter-experiment repeatability threshold: only groups diverging at the left of this line define robust genotypes, while splits occurring at the right may represent experimental rather than biological variation. Inferred genotypes K1-K14 (Table 1) are indicated on the right. Groups A and B are as in Table 1 and Fig. 2. The analysis was performed using Bionumerics 5·10. The following replicas were included in the analysis. (1) Cloned parasites from isolates MHOM/NP/03/BPK206/0 and MHOM/NP/03/BPK190/0, indicated as MHOM/NP/03/BPK206/0-cl10 and MHOM/NP/03/BPK190/0-cl3. (2) MHOM/NP/03/BPK190/0-cl3-M12, isolated after 12 months of regular culturing of MHOM/NP/03/BPK190/0-cl3. (3) MHOM/NP/03/BPK282/0-cl1-S, MHOM/NP/03/BPK282/0-cl3-S, and MHOM/NP/03/BPK206/0-cl10-S, cloned parasites maintained in modified Eagle's medium (Invitrogen, Carlsbad, CA, USA) supplemented with 20% (v/v) heat-inactivated fetal calf serum (PAA laboratories GmbH, Pasching, Austria) as reported by Decuypere et al. (Reference Decuypere, Vanaerschot, Rijal, Yardley, Maes, De Doncker, Chappuis and Dujardin2008). (4) Duplicate DNA samples from isolates MHOM/NP/02/BPK173/0 (DNA extracted from a different subculture MHOM/NP/02/BPK173/0#1), MHOM/NP/03/BPK282/0 (same DNA extraction MHOM/NP/03/BPK282/0#1), and MHOM/NP/03/BPK206/0 (same DNA extraction MHOM/NP/03/BPK206/0#1 to 5).
Nuclear versus kDNA genotyping
As listed in Table 1, combining all nuclear markers showed 8 different genotypes (He=0·54), whereby 4 microsatellite markers and 2 other loci distinguished group A from group B, and only microsatellite Li23-41 (He=0·53) discriminated beyond this group level. kDNA minicircle PCR-RFLP provided 14 genotypes (He=0·92) which, combined with nuclear markers, resulted in 18 different genotypes (He=0·94) among the 34 clinical isolates.
DISCUSSION
In the present study, we have explored nuclear and kinetoplast DNA genetic markers for discriminating L. donovani strains in the VL endemic region of Nepal, which can be used to study the micro-epidemiology of parasite populations. From all studied nuclear assays, microsatellite data provided the most useful information, amounting to 8 different genotypes in the studied population, named M1-M8. When comparing with the study of Alam et al. (Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009), we found 3 extra genotypes, and results from the 19 strains common to both studies are in agreement except for 1 additional allele in MHOM/NP/02/BPK043/0. Other nuclear markers were either monomorphic, or genetically linked to the microsatellite regions, thereby not providing additional information. Seven of the 8 microsatellite genotypes were derived solely from microsatellite Li23-41, while the remaining 4 microsatellite markers discriminated the groups A and B also found with markers haspB and cpb-intragenic. Judging from the common MLMT loci analysed in our study and in that of Alam et al. (Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009), all strains of the Indian zymodeme MON2 (=LON41) belong either to MLMT genotypes M1 or M4, both present within group A. Strains of the second zymodeme circulating in India, MON38, could not be linked to any of the MLMT types here described. Likewise, the haspB amplicon size found in group A is present in MON2 (Haralambous et al. Reference Haralambous, Antoniou, Pratlong, Dedet and Soteriadou2008), but no MON38 amplicon size (515 bp) was found in our samples. The microsatellite and haspB profile of group B, to which 2 samples in our study belong, was not documented by Haralambous et al. (Reference Haralambous, Antoniou, Pratlong, Dedet and Soteriadou2008) and Alam et al. (Reference Alam, Kuhls, Schweynoch, Sundar, Rijal, Shamsuzzaman, Subba Raju, Salotra, Dujardin and Schönian2009). Given these results, it is possible that groups A and B represent different zymodemes, but this can currently not be confirmed.
The kDNA minicircle PCR-RFLP fingerprint patterns clustered into 14 reproducible kDNA genotypes K1-K14, thereby discriminating more strains as compared to nuclear data. Genotypes K3, K9, K10, and K13 were found in combination with different nuclear genotypes, while the remaining 10 kDNA types were linked to one particular nDNA type. Even though this might suggest linkage of some kDNA minicircles with particular nuclear genotypes, one should bear in mind that a kDNA fingerprint profile originates from a complex minicircle population, and therefore assigned kDNA genotypes are largely dependent on the technology used, and the defined groups are, to some extent, arbitrary as each profile is in essence unique. This is in contrast to, for instance, data from microsatellites, the repeat number of which can be determined objectively and unambiguously. When looking at kDNA, minor variants in the minicircle population of a particular isolate may remain hidden in the analysis.
When looking at nuclear DNA genotypes, N1 is the dominant genotype, found in the entire region analysed. Genotypes N2, N3, N4, N7, and N8 are unique, and seem to represent isolates that have not spread over the human population. N5 was found in 2 villages from the Sunsari district less than 30 km apart, and N6 was isolated only from 4 patients living in neighbouring villages less than 500 m apart. A striking feature of genotypes N5 and N7 is that they are characterized by a heterozygous microsatellite locus Li23-41. In addition to the widespread allele of 31 repeat units, both genotypes possess a second allele (32 and 26 repeats respectively) not found in any other sample. These heterozygous genotypes could have originated from either a hybridization event, or by a de novo generated allele.
Even though sampling is limited and biased, and hence most probably is not representative of the parasite diversity in the region, our data suggest that some nuclear genotypes are geographically confined or even unique to individual patients, which makes them especially suitable for studying the micro-epidemiology of the parasite population. These trends were observed also from kDNA PCR-RFLP. As extreme examples, genotype K13 was found in 3 patients spread from East to West over the entire study area, thereby representing a geographically widespread genotype. On the other extreme, genotype K9 was found in 6 patients from neighbouring villages that are less than 1 km apart, and was not found anywhere else. K11 is less confined, but restricted to the southern part of the two neighbouring districts Sunsari and Morang. Taken together, it is clear that the use of kDNA PCR-RFLP in combination with microsatellite Li23-41 and any marker discriminating group A from group B can be a powerful tool for parasite tracking in the region. Even though the number of genotypes and expected heterozygosity of kDNA genotypes exceeded that from nuclear data, the latter give useful information as some kDNA genotypes (K3, K9, K10, K13) were found in combination with different nuclear genotypes. As these trends were observed on the basis of a limited number of parasites isolated at a tertiary medical care centre, additional sampling is needed to substantiate these observations.
Genotypes defined on the basis of kDNA minicircle PCR-RFLP data have been used extensively throughout literature for the purpose of strain identification (Morales et al. Reference Morales, Chicharro, Ares, Cañavate, Barker and Alvar2001; Botilde et al. Reference Botilde, Laurent, Quispe Tintaya, Chicharro, Cañavate, Cruz, Kuhls, Schönian and Dujardin2006; Cortes et al. Reference Cortes, Mauricio, Almeida, Cristovão, Pratlong, Dedet and Campino2006; Laurent et al. Reference Laurent, Rijal, Yardley, Croft, De Doncker, Decuypere, Khanal, Singh, Schönian, Kuhls, Chappuis and Dujardin2007). However, the defined genotypes are extremely sensitive to experimental variations, the reason being that one starts from PCR amplification of a population of different minicircles (Singh et al. Reference Singh, Curran, Middleton and Rastogi1999). As the different minicircles compete during the PCR reaction, even slightly different reaction mixes or cycling conditions can result in a different outcome. It is worth mentioning that the use of nested PCR to increase sensitivity can generate patterns that deviate significantly from a single PCR (data not shown). For these reasons kDNA genotyping protocols have to be extremely detailed to reduce inter-experiment variations to an absolute minimum when assessing relationships among parasites, and must include detailed instructions as to the analysis of the generated fingerprints. Only with these restrictions, we were able to reliably and reproducibly assign the genotypic groups shown in Fig. 3, as determined by repeated analysis of different DNA preparations originating from the same parasite isolate. Because of these limitations, inter-laboratory comparisons of kDNA genotypes are only possible after an extensive standardization effort. As the method is also sensitive to the amount of starting material in the PCR, kDNA typing in clinical samples is not recommended so far, unless the presence of PCR inhibitors has been excluded, and the parasite load has been determined to allow standardization of the input amounts. In addition, kDNA minicircle populations may change over time in a given clone and, in that respect, also introduce a temporal dimension to parasite diversity. The fact that a clone of isolate BPK190 could still be reliably assigned to its original genotype after being kept in culture for 1 year, indicates that our kDNA genotyping protocol is sufficiently robust to account for these shifts.
In conclusion, we have shown the potential of nuclear and kDNA genotypic markers of L. donovani as a tool in studying the micro-epidemiology of the parasite population in the Terai region of Nepal. Nuclear markers were shown to be less discriminative than kDNA, but to allow inter-laboratory comparisons. However, there is a need to improve the sensitivity of some PCRs for use on clinical samples. Our tools can be applied for determining the origin of infection, distinguishing relapse from re-infection, documenting the emergence of new parasite variants, and monitoring the spread of existing ones. Finally, additional techniques such as full-genome sequencing (www.leishrisk.net/gemini), AFLP (amplified fragment length polymorphism) (Vuylsteke et al. Reference Vuylsteke, Peleman and Van Eijk2007), and application of the available methods on a large number of strains from the whole Indian subcontinent, could provide a better understanding of the dynamics and trends of parasite circulation.
ACKNOWLEDGEMENTS
We especially thank laboratory and clinical personnel of BPKIHS for sample collection and primary isolation of the parasites. We also thank Professor Gabi Schönian (Institut für Mikrobiologie und Hygiene, Berlin, Germany) for her support and suggestions to validate the MLMT protocol in our laboratory. Saskia Decuypere and Jean-Pierre Wenseleers (Institute of Tropical Medicine Antwerp) are acknowledged for their help with the figures.
FINANCIAL SUPPORT
This study was funded by the European 6th Framework Program INCO-DEV Kalanet Project (Contract no. INCO-CT 2005-01537). G.V.D.A. was supported by the 3rd Framework Program of the Belgian Development Cooperation with ITMA.
