1. Epidemiology and Global Health
  2. Microbiology and Infectious Disease

Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent

  1. Hideo Imamura
  2. Tim Downing
  3. Frederik Van den Broeck
  4. Mandy J Sanders
  5. Suman Rijal
  6. Shyam Sundar
  7. An Mannaert
  8. Manu Vanaerschot
  9. Maya Berg
  10. Géraldine De Muylder
  11. Franck Dumetz
  12. Bart Cuypers
  13. Ilse Maes
  14. Malgorzata Domagalska
  15. Saskia Decuypere
  16. Keshav Rai
  17. Surendra Uranw
  18. Narayan Raj Bhattarai
  19. Basudha Khanal
  20. Vijay Kumar Prajapati
  21. Smriti Sharma
  22. Olivia Stark
  23. Gabriele Schönian
  24. Harry P De Koning
  25. Luca Settimo
  26. Benoit Vanhollebeke
  27. Syamal Roy
  28. Bart Ostyn
  29. Marleen Boelaert
  30. Louis Maes
  31. Matthew Berriman
  32. Jean-Claude Dujardin  Is a corresponding author
  33. James A Cotton  Is a corresponding author
  1. Institute of Tropical Medicine, Belgium
  2. Wellcome Trust Sanger Institute, United Kingdom
  3. National University of Ireland Galway, Ireland
  4. BP Koirala Institute of Health Sciences, Nepal
  5. Banaras Hindu University, India
  6. University of Western Australia, Australia
  7. Charité Universitätsmedizin Berlin, Germany
  8. University of Glasgow, United Kingdom
  9. Northeastern University, United States
  10. Université Libre de Bruxelles, Belgium
  11. Indian Institute of Chemical Biology, India
  12. University of Antwerp, Belgium
https://doi.org/10.7554/eLife.12613
Share this article
Cite this article
  1. Hideo Imamura
  2. Tim Downing
  3. Frederik Van den Broeck
  4. Mandy J Sanders
  5. Suman Rijal
  6. Shyam Sundar
  7. An Mannaert
  8. Manu Vanaerschot
  9. Maya Berg
  10. Géraldine De Muylder
  11. Franck Dumetz
  12. Bart Cuypers
  13. Ilse Maes
  14. Malgorzata Domagalska
  15. Saskia Decuypere
  16. Keshav Rai
  17. Surendra Uranw
  18. Narayan Raj Bhattarai
  19. Basudha Khanal
  20. Vijay Kumar Prajapati
  21. Smriti Sharma
  22. Olivia Stark
  23. Gabriele Schönian
  24. Harry P De Koning
  25. Luca Settimo
  26. Benoit Vanhollebeke
  27. Syamal Roy
  28. Bart Ostyn
  29. Marleen Boelaert
  30. Louis Maes
  31. Matthew Berriman
  32. Jean-Claude Dujardin
  33. James A Cotton
(2016)
Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent
eLife 5:e12613.
https://doi.org/10.7554/eLife.12613
  2. Download BibTeX
  3. Download .RIS

Abstract

Leishmania donovani causes visceral leishmaniasis (VL), the second most deadly vector-borne parasitic disease. A recent epidemic in the Indian subcontinent (ISC) caused up to 80% of global VL and over 30,000 deaths per year. Resistance against antimonial drugs has probably been a contributing factor in the persistence of this epidemic. Here we use whole genome sequences from 204 clinical isolates to track the evolution and epidemiology of L. donovani from the ISC. We identify independent radiations that have emerged since a bottleneck coincident with 1960s DDT spraying campaigns. A genetically distinct population frequently resistant to antimonials has a two base-pair insertion in the aquaglyceroporin gene LdAQP1 that prevents the transport of trivalent antimonials. We find evidence of genetic exchange between ISC populations, and show that the mutation in LdAQP1 has spread by recombination. Our results reveal the complexity of L. donovani evolution in the ISC in response to drug treatment.

https://doi.org/10.7554/eLife.12613.001

eLife digest

The parasite Leishmania donovani causes a disease called visceral leishmaniasis that affects many of the world's poorest people. Around half a million new cases develop every year, but health authorities lack safe and effective drugs to treat them. Up to 80% of these cases occur in the Indian subcontinent, where devastating epidemics have occurred in the last decades.

One reason these epidemics continue to occur is that the parasites develop genetic mutations allowing them to adapt to and resist the drugs used to kill them. As there are few existing drugs that can kill L. donovani, it is crucial to understand how drug resistance emerges and spreads among parasite populations.

Imamura, Downing, Van den Broeck et al. have now investigated the history of visceral leishmaniasis epidemics by characterising the complete genetic sequence – or genome – of 204 L. donovani parasite samples. This revealed that the majority of parasites in the Indian subcontinent first appeared in the nineteenth century, matching the first historical records of visceral leishmaniasis epidemics.

The genomes show that most of the parasites are genetically similar and can be clustered into several closely related groups. These groups first appeared in the 1960s following the end of a regional campaign to eradicate malaria. The most common parasite group is particularly resistant to drugs called antimonials, which were the main treatment for leishmaniasis until recently. These parasites have a small genetic change that scrambles most of a protein known to be involved in the uptake of antimonials.

Parasites may also be able to develop resistance to drugs through additional mechanisms that allow them to produce many copies of the same gene. These mechanisms could allow the parasites to rapidly adapt to new drugs or changes in the populations it infects. The work of Imamura et al. looks only at parasites isolated from patients then grown in the laboratory, so further research is now needed to explore how variable the Leishmania genome is in both of the parasite’s hosts: humans and sandflies.

Imamura et al.’s study reveals how L. donovani has spread throughout the Indian subcontinent in fine detail. The genome data can be used to create simple molecular tools that could form an "early warning system" to track the success of disease control programs and to determine how well the current drugs are working.

https://doi.org/10.7554/eLife.12613.002

Introduction

Parasites of the Leishmania donovani species complex cause visceral leishmaniasis (VL), the most severe presentation of leishmaniasis that is usually fatal if untreated. There are probably between 200,000 and 300,000 VL cases annually (Alvar et al., 2012), leading to as many as 50,000 deaths per year (den Boer et al., 2011; Lozano et al., 2012). VL is widespread in both the New and Old Worlds (Pigott et al., 2014), but as much as 80% of the global VL burden occurs in the Indian sub-continent (Alvar et al., 2012). Recent intensified control efforts have led to a notable decline in cases (Chowdhury et al., 2014) but the problem is not yet eliminated. VL is a key neglected tropical disease, affecting the poorest regions of the world and the poorest communities within these regions (Boelaert et al., 2009). VL was first reported in the Indian sub-continent (ISC) in the 1820s, but initially confused with malaria until the discovery of L. donovani in 1903 (Gibson, 1983). Although VL was nearly eliminated from the ISC in the 1960s (Thakur, 2007) by antimalarial spraying campaigns with DDT, it re-emerged in 1977 and has caused several subsequent major epidemics (Dye and Wolpert, 1988). Widespread chemotherapy for VL in the region has been ongoing since the 1820s, initially with quinine and other drugs, followed by extensive use of the trivalent antimonial SbIII (1915) and compounds of the less toxic pentavalent SbV (1922) such as sodium stibogluconate (SSG), and since 2005 with miltefosine (MIL) that is freely supplied through a government-subsidized control program. The parasite developed resistance to both SbIII and SbV, and after ten years of clinical use there has been a notable decline in MIL efficacy (Rijal et al., 2013; 2007; Sundar et al., 2012).

Leishmania parasites can re-shape their genome rapidly in vitro in response to stress (Leprohon et al., 2009), suggesting structural variation is an important feature by which they can rapidly adapt to changing environmental conditions and drug pressure. However, there is little data on the diversity of clinical Leishmania populations or how they evolve during treatment. While an extensive literature has made use of molecular methods to study the population genetics of Leishmania (e.g. Alam et al., 2009; Lukes et al., 2007; Mauricio et al., 2006; Schonian et al., 2008), existing genetic markers have relatively poor resolution, and in particular L. donovani within the ISC show very little genetic differentiation based on these approaches (Alam et al., 2009; Downing et al., 2012). Whole-genome sequence data has the potential to show significant population structure within the ISC, and also allows us to identify changes in genome structure.

Here we report the genome sequences of 204 L. donovani isolates (Figure 1, Supplementary file 1), obtained from VL patients between 2002 and 2011 from regions in Nepal (N=98), India (N=98) and Bangladesh (N=8) that represent the epicentre of the on-going VL epidemic in the ISC (Figure 1a).

Figure 1 with 2 supplements see all
Download asset Open asset
History and geography of Indian subcontinent L. donovani.

(a) Location of the patients from which the 204 L. donovani genomes were isolated, and of historical Kala-Azar outbreaks. Genetic groups of the parasite isolates are indicated by the colour of the dots representing them, matching those in Figure 2a,c. Sampling dates and locations are summarised in Figure 1—figure supplement 1, and detailed information about each strain including GPS coordinates are given in the source data file. Citations are to historical primary literature reviewed and cited in (Gibson, 1983). Posterior probability distributions of estimated ages for the oldest split in (b) the main population in Bihar and Nepal and (c) the ISC5 group associated with Sb resistance. Dark shading shows estimates under a strict molecular clock, light shading from relaxed molecular clock and lines show relaxed clock results with Bangladeshi and putative hybrid isolates included. (d) Estimated effective population size through time for ISC5 population (green) and the rest of the parasite population (black/grey). Lines show median of posterior distributions, dark and light shading cover 50% and 95% of the posterior density respectively. Dates for all splits on this phylogeny and other results of phylogeographic analysis are shown in Figure 1—figure supplement 2.

https://doi.org/10.7554/eLife.12613.003

Results

Calling variants against a reference genome assembly for a Nepalese L. donovani strain (BPK282/0cl4; Figure 2—figure supplement 1), we identify three divergent genetic lineages circulating in this region (Figure 2b): a core group of 191 closely related parasites found in the highly endemic lowlands of all three countries, a small population of 12 Nepalese isolates found most frequently in the highlands (ISC1) and a single divergent Nepalese isolate (BPK512/0cl9) (Downing et al., 2011). These two main groups show fixed differences at 45,743 sites (Supplementary file 2, table a), and two previously sequenced Sri Lankan L. donovani isolates (Zhang et al., 2014) were more closely related to the core population (21,546 fixed differences) than to ISC1 (45,743 fixed differences). Parasites within each group show little SNP variation with only 5,628 variable sites in ISC1 and just 2,418 sites varying within the core population (Supplementary file 2, table b) and correspondingly few SNPs in protein-coding regions (Supplementary file 2, table c). Core population isolates differ at an average of 88.3 nucleotide sites with an average nucleotide diversity of 9.7 per Mb (Supplementary file 2, table d).

Figure 2 with 4 supplements see all
Download asset Open asset
Genealogical history of L. donovani from the ISC.

(a) Maximum-likelihood tree based on SNPs called for 191 strains (see Figure 2—figure supplement 1) from the core population in the Indian subcontinent. Samples are coloured by population assignment, with putative hybrid strains not clustered in the main groups in black. Further analysis confirms the hybrid ancestry of some of these isolates (Figure 2—figure supplement 2). (b) Unrooted phylogenetic network of the L. donovani complex based on split decomposition of maximum-likelihood distances between isolates described here, reference genome isolates and two published Sri Lankan isolates (Zhang et al., 2014). (c) Model-based clustering of 191 isolates from the core population reveals six discrete monophyletic groups, and some groups and other samples of less certain ancestry. Coloured bars show the fraction of ancestry per strain assigned to a given cluster, with colours assigned to the population most closely related to each cluster. More detailed population clustering analysis shows largely congruent results (Figure 2—figure supplements 3 and 4).

https://doi.org/10.7554/eLife.12613.006

While a panel of microsatellite markers shows no variation between isolates from the core population (Downing et al., 2012), we reveal significant spatial and temporal genetic structure within this group despite this extremely low level of overall diversity (Figure 2a,c; Supplementary file 2, table e). Phylogenetic and clustering methods identify six congruent monophyletic groups (ISC2-7). Three other groups (ISC8-10) and 21 ungrouped isolates had more complex and less certain evolutionary histories (Figure 2a,c; Supplementary file 2, table f). Most of the ISC groups are present throughout our sampling window (2002–2011), and many are present in both India and Nepal (Figure 1—figure supplement 1). There are some exceptions: ISC7 represents a recent radiation (first observed in 2006) with almost no diversity (20 unique SNPs; Supplementary file 2, table g; π=1.8 per Mb) and is restricted to India, while ISC6 is an older and more diverse group restricted to Nepal (π=12.2 per Mb). We observe subsequent evolution within some groups: ISC5 is distinguished from other groups by just 32 SNP sites (Supplementary file 2, table h), but contains a subgroup with multiple novel SNPs and lower somy (Supplementary file 2, table i).

Bayesian phylogenetic models in an explicit temporal and spatial framework revealed that the core population diverged in the mid 19th century (Figure 1b), matching the dates of the earliest reports of large-scale VL outbreaks in the ISC (Gibson, 1983) and thus suggesting that modern lowland parasites descend from these early epidemics. Within the core population, the Indo-Nepalese population itself appeared around 1900 (Figure 1—figure supplement 2), almost certainly in India (0.89 posterior probability), matching the dates of the first reported outbreaks in Bihar (Gibson, 1983), more precisely in Purnea (Figure 2d). Most subsequent diversification is more recent, with many groups (ISC2 & ISC4-6) radiating from the 1960s (Figure 1b), coinciding with the end of the DDT spraying campaign. The estimated rate of migration from India to Nepal in the Core 191 group was significantly greater than that from Nepal to India, suggesting that India acts as a source population seeding the Nepalese epidemic (Figure 1—figure supplement 2).

A lack of linkage disequilibrium decay between SNP pairs with genomic distance in the core population (r2~0.33 at 5–1,400 kb) reflects a lack of detectable recombination within the six main genetic groups (ISC2-7) across the entire genome (Supplementary file 2, table j). While the low number of SNPs varying within the core population limits our power to detect recombination, we find compelling evidence of hybridisation among eight of the samples not assigned to any of the ISC groups (Figure 2—figure supplements 2–4). The identity of these isolates as hybrids and our assignment of other isolates to groups is supported by allele-frequency based methods (f-statistics), which should be robust to gene flow between groups (Supplementary file 2, tables e,f) and population structure analysis based on haplotype sharing (Supplementary files 2, tables k–n). The four-allele test also confirms that recombination is largely restricted to these hybrids (Supplementary file 2, table o). These isolates appear to result from multiple independent recent hybridizations between distinct ancestors of either ISC5 and ISC6, ISC5 and ISC7, or ISC6 and ISC7 (Figure 2—figure supplement 2).

We detect extensive variation in the structure of these L. donovani genomes. Local copy-number variants (CNVs) cover ~11% of the genome. These include sporadic gene duplication, dynamic tandem gene array sizes (Figure 3—figure supplement 1) and long sub-telomeric amplifications/deletions, the latter generally spanning whole transcription units. While structural variation in Leishmania is often considered a transient adaptation, particularly to culture conditions in vitro, we find striking conservation of many CNVs across all core population groups here. Two multigenic intra-chromosomal duplicated regions (MAPK1 and H-locus; Downing et al., 2011) are present in variable numbers in all core population isolates but are absent in ISC1 (Figure 3b,c; Figure 3—figure supplement 2). Conserved heterozygous SNPs in both of these structural variants confirm that these regions have duplicated once and been maintained throughout the evolution of this population. All known genes on these duplicated regions are associated with virulence (MAPK1, ASS, sAcP; Fernandes et al., 2013; Lakhal-Naouar et al., 2012; Wiese, 1998) or drug resistance (Brotherton et al., 2013), indicating that extensive structural variation allows these parasites to alter local copy number in response to changing environments: both aneuploidy and CNV regulate gene expression (Leprohon et al., 2009). Most isolates are aneuploid (Figure 3—figure supplement 3), even excluding the generally tetrasomic chromosome 31, and almost all chromosomes show some variation in somy (Figure 3a). Aneuploidy (r2=0.15, p=2.7x10-118), CNVs (r2=0.26, p=7.5x10-218) and indels (r2=0.30, p=2.1x10-254) are significantly correlated with SNP variation in the core isolates, suggesting that these variants have appeared gradually during the evolution of the population in the field. Most strikingly, we find two cases of recent epidemic expansions associated with major changes in aneuploidy and heterozygosity (Figure 4). Variation in somy can thus lead to changes in heterozygosity, which could allow selection to eradicate recessive deleterious mutations in the absence of recombination (Roze and Michod, 2010).

Figure 3 with 3 supplements see all
Download asset Open asset
Structural variations in ISC L. donovani.

(a) Stacked barplots per chromosome showing the proportion of ISC strains that are monosomic, disomic, trisomic, tretrasomic or pentasomic for the respective chromosome. A full breakdown of somy per strain is presented in Figure 3—figure supplement 3, and a complete catalogue of other structural variants in Figure 3—figure supplement 1. Violin plots showing the copy number of MAPK1 (b) and H-locus (c) per ISC group, except for ISC1 where these amplicons were absent. These amplicons are intra-chromosomal (Figure 3—figure supplement 2). (d) Tetrameric protein model of the transport protein aquaglyceroporin-1. The C-terminus part that is affected by the 2-nucleotide frameshift found in all ISC5 isolates is shown in magenta. Image was created using PyMOL version 1.50.04 (Schrödinger).

https://doi.org/10.7554/eLife.12613.011
Figure 4
Download asset Open asset
SNP heterozygosity and somy variation in two subclades.

Two subclades show an expansion of polysomic strains from disomic ancestors (below) and an expansion of disomic strains from polysomic ancestors (above). Somy variation per chromosome (1–36; above heatmap) and the total number of heterozygote SNPs (right to heatmap) are shown for each individual strain.

https://doi.org/10.7554/eLife.12613.015

We find no statistically significant association between any individual SNP or structural genetic variant and in vitro SSG resistance, or SSG or MIL treatment outcomes (Supplementary file 2, tables p–r), but the distribution of antimony susceptibility was uneven across different ISC populations (Supplementary file 1). 9 of 11 ISC5 samples tested were highly SbV-resistant and two out of four ISC5-ISC6/7 hybrids tested have intermediate levels of resistance. One variant – a two-bp insertion introducing a frameshift and premature stop codon in the aquaglyceroporin-1 gene (LdBPK_310030, AQP1) – is homozygous in all 52 ISC5 isolates (Table 1), and heterozygous in six hybrids between ISC5 and either ISC6 or ISC7. ISC5 isolates also share other genomic features – such as higher copy number of both the H-locus and MAPK1 amplicons (Figure 3b,c). The H-locus includes MRPA, a gene involved in the efflux of SbIII and associated with drug resistance (Leprohon et al., 2009). Other lines of evidence strongly link AQP1 with antimony resistance. While recent antimonial drugs such as SSG are compounds of pentavalent antimony (SbV), SbV is thought to act mostly as a pro-drug, being reduced to SbIII in both the macrophage phagolysosome (Frézard et al., 2001) and in the parasite itself (Denton et al., 2004; Decuypere et al., 2012). AQP1 is known to assist with SbIII uptake, both genetic and transcriptional changes at this locus have been associated with Sb resistance (Gourbal et al., 2004; Monte-Neto et al., 2015; Mukherjee et al., 2013; Uzcategui et al., 2008) and a homologous transporter is associated with drug resistance in trypanosomes (Baker et al., 2012). Recently, an AQP1 knockout line of Leishmania major was shown to be resistant to SbIII due to reduced uptake (Plourde et al., 2015). The truncated frameshift protein found in ISC5 is predicted to be incapable of forming a functional trans-membrane channel (Figure 3d). We find three other independent frameshifts in AQP1 gene in other antimony resistant isolates, including one in BPK181/12 (ISC6), an isolate taken from a patient following failure of ten months of antimony treatment that was absent in the pre-treatment isolate from the same patient (BPK181/ 0cl11, Table 1).

Table 1

Small indels. The first half of the table summarises the numbers and types of indels detected in each group. The second half of the table shows the proportion of samples within a cluster that share each group-specific coding-region indel.

https://doi.org/10.7554/eLife.12613.016
1. Number of indels
 ISC002 ISC003 ISC004 ISC005ISC006ISC007ISC008ISC009ISC010
Total number ofindels found within each group
586073796555558460
Number of group-specific Indels shared by part of the strains of that group951112708227
Number of group-specific Indels shared by all the strains of that group613300002
Number of group-specific Indels within coding regions
012340131
2. Indels within coding region
Gene IDGene productPositionType
LdBPK_310030AquaglyceroporinLd31_000777420.04
LdBPK_310030AquaglyceroporinLd31_000773521
LdBPK_310030AquaglyceroporinLd31_000766210.04
LdBPK_310030AquaglyceroporinLd31_000809920.11
LdBPK_291860Putative historie H2ALd29_0816454-20.25
LdBPK_040410Conserved hypothetical proteinLd04_015549130.08
LdBPK_070540Conserved hypothetical proteinLd07_0230487-30.12
LdBPK_190080Conserved hypothetical proteinLdl9_001515110.04
LdBPK_261790Conserved hypothetical proteinLd26_065174840.11
LdBPK_301000Conserved hypothetical proteinLd30_0311376-10.02
LdBPK_310690Conserved hypothetical proteinLd31_0241951-30.11
LdBPK_332580Conserved hypothetical proteinLd33_099596010.04
LdBPK_366590Conserved hypothetical proteinLd36_2473775-30.25
LdBPK_110650hypothetical, unknown functionLdll_0245832-30.56
LdBPK_292330hypothetical, unknown functionLd29_1008496-30.08

We propose that the AQP1 truncation is associated with antimonial resistance in the ancestor of ISC5, and has been transmitted to a group of hybrid parasites. The ISC5 lineage emerged following the end of the DDT campaigns but then proliferated quickly in the 1970s (Figure 1c,d), when antimonial dosage had to be doubled because of its declining efficacy. The persistence of this lineage beyond the era of Sb treatment perhaps reflects the increased fitness (Vanaerschot et al., 2013) of Sb-resistant parasites. These observations tally with a stronger signature of purifying selection on the ISC5 lineage, measured as a lower rate of derived allele accumulation compared to other ISC populations, which may be a consequence of higher historical exposure to drug stress (Supplementary file 2, table s). Sb resistance is also present in other genetic groups, with 4 out of 15 ISC4 lines tested in vitro being SbV-R, indicating resistance has emerged independently and recently multiple times in ISC L. donovani, and that other genetic variants responsible for SbV resistance must be present in this population. Indeed, other SbV resistance mechanisms are known in this population: previous work has shown that two resistant strains from ISC4 (BPK087 and BPK190) show significantly decreased transcription of an AQP1 locus encoding a wildtype protein sequence (Decuypere et al., 2005), and BHU764 combines a different indel mutation in AQP1 and reduced expression of MRPA, an efflux transporter of SbIII (Mukhopadhyay et al., 2011). The failure of any single resistance locus to sweep through this population may reflect the low level of gene flow and the presence of a large reservoir of untreated asymptomatic cases (Ostyn et al., 2011).

Discussion

We have shown that genomic data can retrospectively unravel the evolution and epidemiology of this parasite population, and gain new insight into possible mechanisms of drug resistance against a background of extensive variation in genome structure. We report the first analysis of the structure and history of a Leishmania population, aligned with clinical and epidemiological records, enabled by the higher resolution of genome sequence data than other genotyping approaches. These data have allowed us to describe a mechanism of resistance to one of the most ancient drugs used in the human pharmacopeia, antimonials, not only identifying a key locus, but also showing the epidemiological dynamics of a population carrying a loss-of-function variant at this locus.

Continued genetic surveillance of parasite populations is key to rapidly identify and respond to the emergence of treatment failure. In the recent emergence of artemisinin resistance in Plasmodium falciparum, genomic data has led to the identification of the major locus underlying resistance (Ariey et al., 2014; Cheeseman et al., 2012), revealed the genetic architecture of resistance (Miotto et al., 2015) and shed light on the population genetic context in which resistance is appearing (Miotto et al., 2013). Genomic surveillance is playing a key role in defining the geographic boundaries of the spreading artemisinin-resistant population. Failure of anti-Leishmania chemotherapy could become a similar public health emergency: miltefosine has shown reduced efficacy in both India (Sundar et al., 2012) and Nepal (Rijal et al., 2013). While amphotericin B is now being used against visceral leishmaniasis in ISC, few alternative treatments are available, and continued genomic surveillance will facilitate tracking the response of the Leishmania population to continued use of these drugs.

Monitoring drug resistance in clinical settings is challenging: the data set we present was generated as part of a five-year collaboration between clinicians in the endemic countries, parasitologists and genome biologists. This collaboration is critical in generating data that reflects the evolution of parasite populations in close to real time and as such is directly applicable in a public health context. The data we present here provide baseline information on the diversity of Leishmania donovani in the ISC that will contribute to future studies of drug resistance and epidemiology of this population. Our results show the promise of genomic surveillance for other Leishmania populations, where patient symptoms, the parasites involved and the main treatment modalities all differ from those in the ISC (Sundar and Chakravarty, 2015).

Materials and methods

Sample collection

Request a detailed protocol

The ethics committee of (i) the Nepal Health Research Council, Kathmandu, (ii) the Institute of Medical Sciences, Banaras Hindu University (BHU), Varanasi, India and (iii) the corresponding bodies at the Institute of Tropical Medicine of Antwerp and the Antwerp University, Belgium, reviewed and approved the study protocol. Informed written consent was obtained from each patient or his/her guardian for those <18 years of age. All the patients and caretakers/parents had the study purpose explained to them in local language.

A total of 204 parasite isolates were obtained from clinically confirmed VL patients in the high endemic regions of the Indian subcontinent (ISC) by the B.P. Koirala (BPK) Institute of Health Sciences in Dharan (Nepal, Terai, N=98), the Kala-azar Medical Research Center in Muzaffarpur (India, Bihar, N=98) and the Mymensingh Medical College in Mymensingh (Bangladesh, BD, N=8). The Indian and Nepalese isolates were collected as part of a multi-center collaborative project to investigate drug resistance in ISC and were all typed as Leishmania (Leishmania) donovani. Complete clinical and epidemiological data were available for the Indian and Nepalese isolates (Supplementary file 1).

The 204 L. donovani isolates were obtained from confirmed visceral leishmaniasis patients in previous clinical studies as described elsewhere (Rijal et al., 2013; 2007). PCR-RFLP of the cysteine proteinase gene (Quispe Tintaya et al., 2004) typed all isolates as Leishmania donovani. Strain names consisted of 2–3 letters that indicated the location of isolation (BD, BHU, BPK), 2–4 digits that indicated the patient number in that location, a forward slash followed by 1–2 digits that indicated when the sample was isolated (0: before treatment, 1: 1 month after treatment, etc) and optionally the number of the parasite clone if the strain was cloned (clone one is listed as 'cl1'; clone two is listed as 'cl2', etc). Cloning was performed using the micro-drop method (Van Meirvenne et al., 1975). Patient treatment outcome was monitored at the end of treatment and at 3, 6 and 12 months post-treatment). Treatment non-response was defined as a case with positive parasitology at the end of the treatment course. Patients who were successfully cured at the end of treatment but in whom symptoms re-emerged within the 12 month follow-up period were classified as relapse cases. Patients who were cured at the end of treatment and remained cured within the 12 month follow-up period were classified as definite cures. If patients were lost to follow-up, the last known treatment outcome was recorded. Seven pre- and post-treatment samples coming from the same patients were obtained. Patient treatment outcome after treatment with miltefosine (MIL) and pentavalent antimonials (SSG) was monitored during 12 months (at the end of treatment, month 3, month 6 and month 12 after treatment).

Sample phenotyping

Request a detailed protocol

50 strains were phenotyped for their susceptibility to SSG using a standardized in vitro susceptibility test as described elsewhere (Downing et al., 2011; Rijal et al., 2007). An SSG-susceptible reference strain (BPK206/0cl10) was included in each assay. The classification into resistance and susceptible strains was determined by calculating the activity index (AI): the ratio of the EC50 of the strain in question versus the EC50 of the susceptible reference strain. AI values clustered strongly, with most strains showing an AI≤1 (25; classified as SSG-sensitive) or ≥6 (18; classified as SSG- resistant). A few strains (7) showed AI values around 3 and were considered as showing intermediate resistance.

Genome sequencing

Request a detailed protocol

DNA isolation, sample preparation, DNA quantification and DNA library preparation were done as outlined previously (Downing et al., 2011). 100 bp paired-end sequence reads were generated (median coverage 44 per sample) with the Illumina Hiseq 2000 platform according to standard protocols. Read data are available under study ERP000140 at the European Nucleotide Archive (http://www.ebi.ac.uk/ena/data/view/ERP000140).

DNA read mapping

Request a detailed protocol

Reads were mapped to the reference L. donovani genome BPK282/0cl4 using Smalt v5.7 (www.sanger.ac.uk/resources/software/smalt/). Options for exhaustive searching for alignments and random assignment of repetitively mapped reads were used to properly estimate read coverage. Non-mapping read exclusion, read file merging, sorting and elimination of PCR duplicates were implemented with Samtools v0.1.18 and Picard v1.85.

Reference genome masking

Request a detailed protocol

The reference genome was masked at regions of the genome that were repetitive, duplicated, close to contig edges, structurally variable, or potentially mis-assembled. Five criteria masked a total of 6,358,203 bp out of the 32,444,998 bp reference genome sequence for L. donovani BPK282/0cl4, resulting in SNPs being called at 26,086,795 or 80.4% of the nuclear genome. Criteria were: 1. Manually identified repeats, commonly duplicated or deleted regions, regions with excessive rates of common SNPs and non-unique regions (Downing et al., 2011) identified with Gnuplot, the Artemis Comparison Tool, Artemis and Samtools tview (1,740,084 bp). 2. Duplicated regions determined by DNA similarity as Blast v2.2.25 (Altschul et al., 1990) hits between the two reference genome sequences for L. donovani BPK282/0cl4 and L. infantum JPCM5, with E-value less than 10e-20 (2,082,546 bp). 3. Low complexity repeat regions determined by Tantan v0.13 (www.cbrc.jp/tantan/); (2,495,070 bp). 4. 100 bp regions adjacent to each contig edge (1,641,511 bp) – initially 13.8% of candidate SNPs were in these regions. 5. The first 300 bp and last 5 kb of all chromosomes, which are more likely to contain mis-assemblies.

SNP detection using COCALL

Request a detailed protocol

SNPs were ascertained using a consensus calling approach (COCALL) that is based on the framework outlined for the 1000 Genomes project (1000 Genomes Project Consortium, 2012). COCALL applied five different variant detection approaches and combines evidence from them to calculate the support for each genotype. For complete details on the algorithm testing and development, see Appendix 1. In short, this approach avoids bias associated with systematic errors unique to each individual SNP caller by examining their consistency and identifying discordant mutations symptomatic of false positives. The five callers used were FreeBayes v0.9.5, GATK 2.0–38, Samtools Pileup v0.1.16 and Mpileup v0.1.18 based on the DNA read mapping by Smalt, and Cortex v1.0.5.13 based on its own de novo assembly and mapping. In a large population of genetically homogenous strains, superior inference power was achieved by examining the population-wide genotype at each candidate SNP site (i.e. population-based COCALL; Figure 2—figure supplement 1). Candidate SNPs with genotype qualities of 40+ across all five callers were retained. SNPs with population normalized read depth ≤0.5 or ≥1.75 or with multiple derived alleles across the five callers were excluded. Candidate SNPs in soft-masked regions were accepted where the number of callers ≥3.5; those in non-masked regions were kept where the number of callers ≥2.5. SNP sites retained in the final set of retained SNP sites were supported by a mean of 4.5 callers out of 5.

Copy number variant, somy and indel detection

Request a detailed protocol

Chromosomal read depths were computed using a trimmed median read depth (calculated as the median of read depths for sites with depths within one standard deviation of an initial, untrimmed, median read depth of each chromosome) and normalized as the depth per haploid genome as outlined previously (Downing et al., 2011). Somy levels were estimated as the median normalized chromosomal read depths (Downing et al., 2011). Local copy number variants (CNVs) were detected where the local read depth was significantly different from the median depth of approximately 60 samples from ICS4, ISC6 and ISC8 whose depth profile is similar to that of BPK282/0cl4, and were measured with respect to the haploid depth to exclude somy variability. Two CNVs in particular, the MAPK1 and H-locus, were further investigated as they show functions potentially relevant to parasite adaptation (Downing et al., 2011). A quantitative PCR assay in a subset of 46 samples was performed to confirm the copy number variation of the MAPK1 and H-locus amplicons. The nature of the amplification (extra- or intra-chromosomal) was determined by pulsed-field gel electrophoresis (PFGE) and southern blot hybridization comparing two strains that showed differential amplification of these loci (ISC6 strain BPK282/0cl4: amplification; and ISC1 strain BPK026/0cl5: no amplification). To exclude the possibility that the amplicons are a culturing artefact, PCRs using primers that enabled amplification of circular episomes or tandem duplications was also attempted directly on five bone marrow samples from VL patients. Indels were detected using a consensus calling method based on the concordance of results across four tools: Cortex, Freebayes, GATK and Samtools Mpileup. For complete details on the Somy, CNV, indel and episome detection, see Appendix 2.

Haplotype inference and linkage disequilibrium

Request a detailed protocol

Haplotypes were inferred using PHASE v2.1.1 (Stephens et al., 2001): 0.1% of genotypes in the Core 191 and 0.9% in ISC1 had confidence scores <0.95. Haplotypes were inferred with a general recombination rate model (Li and Stephens, 2003) with ten runs, each with a burn-in of 100 steps, 100 iterations and a single MC thinning step and recombination rate estimated for each chromosome. Convergence was examined for each chromosome: recombination estimates were consistent, though there was more variation between phasing runs for chromosome 16 in the core population and consequently inferred haplotypes are less certain for that chromosome. There was no correlation between the mean chromosome copy number and mean recombination rate or PHASE probability values for inferred haplotypes (r2=0.011). While variation in somy is not explicitly accounted for in the phasing process, the rapid flux in the somy levels of aneuploid chromosomes may mean this variation has no effect on haplotype inference. Of 3,567 heterozygous sites, 3,076 (86%) had a PHASE probability of exactly 100% and 437 had PHASE probabilities < 0.95: these lower-confidence haplotypes were masked. Haplotypes for BHU1087/0 were inferred along with the core population. The phased core population SNP set had 2,401 SNPs: 17 singletons were masked. The smaller sample size meant that phasing within ISC1 was less successful: phase was successfully inferred for 2,308 sites using the same (0.95) confidence score threshold: 524 were not phased and were excluded from further analysis. No correlation between phasing confidence score and trisomy or tetrasomy was apparent.

Linkage disequilibrium (LD) was inferred as the correlation in genotypes (r2 values) between SNP pairs using Bcftools v0.1.17 screened with Samtools Mpileup given SNP mapping qualities >30 and base qualities >25. These pairwise r2 values were used to examine genome-wide LD patterns and LD decay with distance. Recombination was confirmed using the four-gamete test (Hudson and Kaplan, 1985). Mean chromosomal estimates of LD in the core population did not correlate with somy level if the tetrasomic chromosome 31 was excluded (r2=0.001) but did if chromosome 31 was included (r2=0.167). Somy had little impact on the variance of LD per chromosome (r2=0.017 with chr31, r2=0.000 without chr31). Variance in somy level across chromosomes had no association with either the mean or variance of LD per chromosome. We calculated zygosity as the probability that a SNP exists at a distance d from a SNP at a site x assuming diploidy (Lynch, 2008). No differences between homozygous and heterozygous SNP clustering measured as a product of chromosomal distance was observed.

Population genomic identification of groups

Request a detailed protocol

L. infantum JPCM5 (MCAN/ES/1998/LLM-877) from Spain and LV9 (MHOM/ET/1967/HU3) from Sudan were used (Downing et al., 2012) for comparison with the L. donovani genomes generated in this study. Variants were called for these two samples using the approaches outlined above. For two additional L. donovani isolates from Sri Lanka (Zhang et al., 2014), we mapped Illumina GAII reads using Smalt v5.7 as above and called candidate SNPs at non-masked regions using Samtools Pileup v0.1.16 (Li et al., 2009), followed by screening steps as above. Two genomes were excluded in the final analyses because sequence reads were of insufficient quality (for MHOM/IN/10/BHU1087/0) or because of a suspected mixed infection (for MHOM/IN/10/BHU790/0). BHU790/0 is distantly related to the core population (most likely ISC3) and appears to be a mixed infection rather than a hybrid because its average read allele frequency of heterozygous SNPs approximates 0.17, whereas most detected hybrids had mean read allele frequencies of 0.4–0.5. Remaining data were used to construct phylogenies using the 211,536 sites containing verified SNPs in the entire sample set (ISC1/2/3/4/5/6/7/8/9/10 and ungrouped, LV9, JPCM5, BPK512/0cl9). JPCM5, LV9, the two Sri Lankan isolates and one sample from our collection (BPK512/0cl9) represented genetically distinct lineages, distinct to both the ISC1 (n=12) and core populations (n=191). Seven SNPs in the core population and ten in ISC1 had multiple derived alleles compared to reference genome sample BPK282/0cl4 (Supplementary file 2, table t). These were included in all diversity analyses, but not those involving phased haplotypes.

Genome-wide phylogenetic trees were constructed with RAxML v8.1.1 (Stamatakis, 2014) using the GTR+G substitution model and 1000 bootstrap replicates for 10 runs for the core population (881 alignment patterns), ISC1 (349 alignment patterns), and all samples including the CL and VL samples from Sri Lanka (Zhang et al., 2014) (2274 alignment patterns). The best fitting substitution model determined using MEGA v6 (Tamura et al., 2011) for the core population was GTR+G. The final phylogenies were visualised using MEGA v6 (Tamura et al., 2011) and Splitstree v4 (Huson and Bryant, 2006). Unrooted haplotype trees for the phased SNPs for each chromosome were constructed from maximum-likelihood distances for the TN93 substitution model using the package Ape (Paradis et al., 2004) v3.1–4 in R version 3.12.

Samples in the core population of 191 isolates were classified using model-based clustering as implemented in Structure v2.3.2.1 (Pritchard et al., 2000) and principal component analysis (PCA) of the allele frequencies. Given a number of genetically distinct clusters (K), samples were probabilistically assigned to a population independent of a mutation model with a prior of 1/K based on the correlation in genotypes of each sample with estimated population allele frequencies. 1≤K≤15 was examined with admixture and incomplete membership allowed to reduce overfitting. We used 105 burn-in steps before a run of 2x106 steps with three independent runs per K to confirm chain convergence. The most likely number of clusters was based on the second-order rate of change of the likelihood function (ΔK, Evanno et al., 2005). At K=4 the groups were composed of ISC2/3/9/10, ISC4, ISC5 and ISC6/7/8. Inter-population differentiation was lower for ISC2/3/9/10 (FST=0.36) compared to the others (0.85<FST<0.98). K=7 was the most probable K value (ΔK=25.8): the groups were composed of ISC2, ISC3, ISC4, ISC5, ISC6/7/8, and ISC9/10 (all FST>0.79) – the 21 ungrouped samples collectively had an FST=0. Most population membership assignments were >0.97 with few ambiguous values (range 0.80–0.97). At K=9, ISC6/7/8 split into ISC6 and ISC7/8 (both FST>0.85). At K=10, ISC7/8 segregated into ISC7 and ISC8.

Inference of historical population sizes, geographic locations and migration rates

Request a detailed protocol

Dated phylogenies, historical population sizes and migration patterns were modelled for the 191 core clinical isolates using BEAST v1.8.1 (Drummond et al., 2012). For molecular clock analyses, hybrid isolates not assigned to any of the ISC groups were removed from the dataset, as were the Bangladeshi outgroups for most analyses. Dates for each were fixed to the month of isolation, with sampling dates for those for which only isolation year data was available estimated during the MCMC but given a uniform prior on sampling ages within that year. Broadly consistent date estimates were obtained under three different models: with an uncorrelated lognormal relaxed clock model and a TVM substitution model and a Bayesian skyride model for population sizes, under the same model but with a strict clock model and finally under a GTR substitution model, with a simple constant population size coalescent model for data including the outgroups. Migration rate estimates were obtained by including a simple continuous-time Markov model of a discrete trait representing the country (Nepal/India) of isolation, so that ancestral states and rates of change in geographical location were estimated along the phylogeny. All analyses were made with a minimum of 8 independent MCMC runs, for 200 million update generations per run. Convergence was assessed by inspection in Tracer v1.6, confirming that at least 5 of the 8 runs had converged to the same stationary distribution of parameters and that this had the highest likelihood. In most analyses, seven or eight chains all converged to the same posterior distribution, but the Bayesian skyride analyses converged more slowly. ESS estimates for almost all parameters across runs was over 500, except for some skyride population size parameters. The first 20 million generations of each MCMC run were removed before combining all converged runs for inference. Historical population sizes were estimated both with the Bayesian skyride model and by transforming lineage-through-time data for all trees in the posterior probability distribution from the strict clock model above using the package Ape (Paradis et al., 2004) v3.1–4 in R version 3.12. To compare population sizes between the drug resistant clade and others, we split ISC5 from other data and removed coalescent events between the ISC groups (the oldest six) to make these comparable with the ISC5 coalescence.

Population genomic identification of admixture using allele frequency correlations

Request a detailed protocol

f-statistics describe the correlation in allele frequencies between populations (Patterson et al., 2012; Reich et al., 2009). The simplest (f2) is simply the sum-of-squares difference in allele frequency between two populations averaged across loci, and so captures the amount of divergence, or branch length between two populations. Two more complex statistics, f3 and f4 are calculated as differences between f2 statistics between groups of 3 and 4 related populations. f3(C;A,B) has the property that, for a population C derived from populations A and B, it is expected to be positive if A,B and C are related by a simple history of divergence and genetic drift, but negative if admixture from A or B has contributed to the genetic composition of population C, while being robust to the details of the relationship between the three populations. In contrast, the value of f4(A,B,C,D) does depend on the evolutionary history of populations A, B, C and D and so can be used to test a proposed relationship: if the four populations are related as ((A,B),(C,D)) the f4 statistic is expected to be zero; for ((A,C),(B,D)) it is expected to be positive and for ((A,D),(B,C)), negative. Finally, if the evolutionary history of three ancestral populations is known, the ratio of two f4 ratios is an estimate of the relative contribution of two potential parental populations to a fourth admixed population, given an outgroup.

Population genomic identification of admixture using haplotype sharing

Request a detailed protocol

Whereas groups ISC2/3/4/5/6/7 seemed clearly defined phylogenetically and by Structure, ISC8/9/10 were not and no simple relatedness among the 21 ungrouped samples was detected. Consequently, we used inferred haplotypes to test whether these represented genetically discrete populations, or whether some of those samples were mixtures of ISC3/4/5/6/7 generated by hybridisation between these groups (Lawson et al., 2012). Chromopainter v0.0.2 and FineStructure v0.0.2 inferred ancestral patterns of haplotype similarity among samples without a prior assumption of a given number of populations or of independence between mutations.

Co-ancestry matrices for the core population were computed using Chromopainter v0.0.2 as the number of segments potentially donated to or received from individual samples, using the phased haplotypes. Recombination rates between pairs of SNPs inferred by PHASE were used for each of 36 unlinked chromosomes. Groups of SNPs on a single chromosome were expected to be exchanged as blocks of different sizes, so a higher number and longer lengths of shared blocks indicate recent common ancestry. The most likely ancestral sample or population was assigned according to its similarity to corresponding segments in a set of donor isolates. Two main datasets were generated by ChromoPainter: a co-ancestry matrix where all 191 could donate to all 191 as recipients (191x191), and another where six representative samples were used as the only donors (191x6: BD09 for ISC2, BPK067/0cl2 for ISC3, BPK087/0cl11 for ISC4, BPK275/0cl18 for ISC5, BPK282/0cl4 for ISC6, BHU200/0 for ISC7). The expected number of chunks was minimised for the six representative samples, with k=80 segments and an effective number of chunks c=0.02. Reducing the number of representative strains to represent distinct groups identified by Structure with smaller K parameters resulted in smaller k and larger c values, suggesting that using six representative samples was the optimal number for discrimination within the core population. Though ISC7 was a subset of ISC6, ISC7 had a large number of fixed SNPs sufficient to differentiate it from ISC6 with Structure, so it was included. For the 191x191 comparison, k=79 segments was expected and the effective number of chunks was lower (c=0.00054) because the total diversity of the donor set per SNP had decreased.

These 191x191 and 191x6 co-ancestry matrices represented the most probable number of segments copied from each donor to each recipient, and also the relative probability of ancestry across the set of donors for each SNP for each sample. The number of donors per recipient was set to 100. 20 expectation-maximisation algorithm iterations was sufficient to maximise the recombination-scaling coefficient (Ne) and copying probabilities with 10<k<1000 iterations across different number of donor samples assuming a minimum recombination rate of 10-15 Morgans/bp. For the 191x191 matrix, the Ne=523.3 and the mutation rate (µ) was 0.000181. For the 191x6 matrix, the Ne=1015.1 and µ=0.000628: Ne and µ were higher because there were more mutations per sample.

The 191x191 matrix was clustered for 106 MCMC (Markov chain Monte Carlo) steps with a burn-in of 10,000 and a skip of 100 steps using FineStructure v0.02 to obtain aggregated expected segment sharing between samples and populations with 100 trees examined per merge step. This distinguished complex ancestral patterns of segment sharing for the strains which Structure could not fully assign to single populations.

To verify FineStructure and Structure results, the correlation in the SNP allele frequencies across samples was examined in the core population for six principal components with p<10-7 using PCA implemented by smartPCA in Eigensoft v4.2 (Price et al., 2006). The first PC separated ISC2 (10.1% of all variation), the second ISC4 (6.4%), the third ISC5 (5.8%), the fourth ISC3 (4.9%), the fifth BPK035/0cl1 and BPK043/0cl2 (4.2%) and the sixth a subset from ISC9/10 (3.9%). This was repeated for the 2353 variable sites in the core population (ISC3/4/5/6/7/8/9/10 and ungrouped samples, n=183) excluding the 8 samples from Bangladesh (ISC2). This partitioned ISC5 (PC1, 7.4%), then ISC4 (6.8%), third ISC3 (5.6%), and fourth BPK035/0cl1-BPK043/0cl2 (4.8%). Eigenstrat and FineStructure PCA results were effectively the same but with some different axis labels – PC1 in the former was PC3 in the latter. FineStructure 191x191 ancestry patterns partitioned ISC4 vs ISC6 over PC1 (16.8% of variation), and ISC5 vs ISC6 over PC2 (15.5%). The next (12.3%) differentiated ISC2, and PCs 4 (6.4%) and 5 (5.7%) separated ISC3. PC6 in FineStructure differentiated the BPK035/0cl1-BPK043/0cl2 pair.

Population genomic identification of drug-resistance elements

Request a detailed protocol

Information on in vitro SbV-resistance was available for 50/191 Core 191 isolates, from which 25 were sensitive and 25 resistant (Supplementary file 1). Links between genetic diversity (SNP, indel, CNV and somy) and in vitro SbV-resistance were assessed using the Fisher Exact test (FET), Mann Whitney U-tests (MWU) and odds ratios (ORs), implemented on 103 CNVs and 17 indels (in 14 genes) as well as 2,392 phased SNPs genotypes. SNPs were assigned to the 5’ and 3’ UTR if they were within 1 kb of the start or end of the gene (respectively). To counter bias associated with the small sample size, FET and MWU were used initially. For the FET, variants were defined as discrete variables: SNPs as 0, 1 or 2 non-reference alleles, and small indels as the diploid number of inserted or deleted basepairs. For the MWU, mutations were considered as a continuous variable such that the somy state was the haploid chromosome state, and CNVs were the haploid copy number times the somy state. The null hypothesis was that there were no significant genetic differences between SbV-R and SbV-S strains (subject to p<0.01). The FETs and MWU were limited by the partial association of different mutations with the phenotypes, so we examined ORs of the derived alleles segregating in multiple ISC populations with 6+ non-reference alleles for which the absolute difference in SbV-R and SbV-S allele frequencies >0.1 using the log-scaled EC50 values. We compared the log-scaled EC50 values of each allele pair using t-tests.

We also examined samples for which the patient was treated with SbV and was either cured or not, samples for which the patient was treated with miltefosine (MIL) and was either cured or not, and also in vitro MIL resistance levels as implemented above for SbV.

Testing for selective processes among ISC populations

Request a detailed protocol

Evidence of historical differences in selective processes on the ancestors of the major ISC populations was assessed as the rate of accumulation of derived alleles. Stronger purifying selection should purge deleterious derived alleles more quickly, detected as an excess of nonsynonymous changes relative to synonymous ones, as previously observed for ISC isolates (Downing et al., 2011). This signature should be most apparent for derived alleles, which should accumulate at a net rate dependent on the historical effective population sizes and selective coefficients. Using L. infantum JPCM5 as the outgroup, the relative abundance of derived alleles in one population that were absent in the other for each ISC population pair (ISC2-7) were determined as the statistic R (Do et al., 2015). The associated ratio R2 denoted the relative rate of homozygous derived allele accumulation between populations. R and R2 should approximate 1 assuming no difference in the strength of selection, and primarily depend on the derived allele frequency per population, so the main confounder was variance in historical effective population sizes among ISC populations. To calculate confidence intervals for these R values that take into account correlation between neighbouring sites, we used a Weighted Block Jackknife by splitting the SNPs according to chromosome (Busing et al., 1999) to counter the extensive linkage disequilibrium between SNPs (Moorjani et al., 2011): discrete chromosomal blocks may still be linked. This was adjusted for the number of SNPs per block to reflect the variability in the relative selective pressure (Kunsch, 1989). A threshold of four times the standard error of these jackknife estimates was used as a criteria for identifying comparisons deviating significantly from expected values (Do et al., 2015).

AQP1 modelling

Request a detailed protocol

A protein model of the intact Leishmania donovani AQP1 from BPK282/0cl4 was created using MODELLER 9.14 (Sali and Blundell, 1993). The template for homology modelling was the crystal structure of the aquaglyceroporin from Plasmodium falciparum in complex with glycerol (PDB code: 3c02) published by Newby and co-workers (Newby et al., 2008). The sequence identity between the target and the template was approximately 33%. PyMOL version 1.50.04 (Schrödinger) was used to generate the biological units for the aquaglyceroporin from Plasmodium falciparum (generation of symmetry mates function in pymol). The C-alpha atoms of chain A, B, C and D of the tetramer template were restrained during homology modeling using MODELLER in order to reduce the number of interatomic distances that needed to be calculated.

Appendix 1 - SNP detection procedure using the COCALL algorithm

1 Consensus variant discovery to reduce tool-specific bias

COCALL (COnsensus of SNP CALL) was designed to overcome a lack of verified SNPs for L. donovani by combining SNP information from different SNP callers using a consensus calling scheme based on the framework outlined for the 1000 Genomes project (1000 Genomes Project et al., 2012). The main premise was that no SNP caller was a priori superior, and that errors among calling tools tended to differ across sites, whereas true variants showed a concordant signal. This approach avoided systematic errors unique to each individual SNP caller by examining their consistency and identifying discordant mutations symptomatic of false positives (Weisenfeld et al., 2014). These disparities among callers arise from their differing strategies for tackling repetitive regions, ambiguous read mapping probabilities and alignments, low complexity regions, tandem repeats, inversions, translocations and regions in which genuine repeats are collapsed in the reference assembly. By combining this with population-level information to scale the prior probability of a mutation existing at a known SNP site, COCALL is optimised for variant discovery in a largely unknown population. This integrated variant discovery scheme inferred sample genotypes both dependently and independently of other samples, which may mitigate bias in rare variant discovery: dependent approaches miss true SNPs and independent ones include many false positives (Han et al., 2014). COCALL was particularly effective when samples were genetically homogenous and the sequence quality and read depth among samples was uniform, as in this study. When these conditions were not met, COCALL was less effective, and population-based methods implemented in GATK, Freebayes or Samtools Mpileup performed better because these methods identified base variations based on aggregated depth across the population, and were less affected by lower quality sequence runs. Computational time for running COCALL may become prohibitively long for samples with high genetic diversity.

This scheme applied five different mutation detection approaches and computed the support for each base variation. COCALL interpreted the results of FreeBayes v0.9.5 (Garrison and Marth, 2012), GATK 2.0-38 (McKenna et al., 2010), Samtools Pileup v0.1.16 (Li et al., 2009) and Samtools Mpileup v0.1.18 based on reads mapped with Smalt v5.7, and also Cortex v1.0.5.13 (Iqbal et al., 2012) based on its own assembly and mapping (Figure 2—figure supplement 1).

2 Attributes of SNP callers used by COCALL

1) Cortex was used as a de novo assembly haplotype-based SNP caller and offered an independent assessment of SNPs because it did not rely on read alignment and mapping by Smalt. For SNP calling, we merged SNPs with scores ≥2.5 identified on the current reference with de Bruijn graphs with k-mers of 31 and 63 bases; the latter were structurally more accurate at a cost of nucleotide level accuracy (Iqbal et al., 2012). This approach detected fewer SNPs overall than other methods, as expected (Weisenfeld et al., 2014). Cortex depended on de novo assemblies, so its accuracy was lower for repetitive regions and collapsed repeats. Cortex did not provide strand bias information, so the strand bias of Cortex SNPs was calculated from based on corresponding strand information from Pileup.

2) Freebayes applied a Bayesian strategy for SNP calling and was used to discover potential SNPs with SNP scores ≥50 using default parameters for strand bias and read depth. Unlike other SNP callers, it detected variation based on the local haplotype across the read length, and so was superior for variant calling at regions with short repeats. It required disomic chromosomes: but for this study this limitation did not affect SNP calling because no strict additional parameters for SNP calling were applied to allow it identify rare variants.

3) GATK provided two SNP calling methods: one for individual samples, and one based on a population of related samples (GATK population caller). A SNP score threshold >50 was used for both. For our sample set, GATK had comparatively lower power to find SNPs at duplicated loci, and the base quality recalibration (BQR) scheme (DePristo et al., 2011) was not applicable due to the low number of known SNPs for training and the low total number of high-frequency candidate SNPs.

4) Samtools Mpileup was applied with a SNP score ≥20, and called fewer SNPs near indels compared to GATK and Freebayes. It assumed disomy, inferred SNPs across a population of samples, and ultimately had a lower false positive rate.

5) Samtools Pileup had a lower SNP detection power at regions with low read coverage, at gaps, and near indels. However, Pileup was more accurate at high-coverage regions with lower mapping quality and produced fewer false positives. Here, it was used for candidate SNPs with a base quality >20 and >10% of reads supporting a SNP being on both the forward and reverse strands to reduce strand bias which was calculated based on pileup read depth. 11.5% of all candidate SNPs had significant strand bias such that the proportion of either forward or reverse reads was >90%, and 4% had reads only for forward or for reverse strands.

3 Screening and consensus SNP calling in populations

Candidate SNPs with mapping qualities of ≥40 across all five callers were retained: approximately 52% of candidate SNPs were excluded due to low base, SNP or mapping quality values. SNPs at sites with normalised read depth ≤0.5 or depth ≥1.75 or with multiple derived alleles across the five callers were excluded. 1.75 was used as the upper read depth limit instead of 2.0, which is often usedas a criterion for duplication, because SNP callers estimated read depth based on base and mapping quality, so even true duplicated regions tended to have depth <2.0. In addition, SNPs present in clusters of more than three per seven bp window were removed as likely to be associated with tandem repeats, small indels and other sequence artefacts. The modal read depth coverage for each candidate SNP site across the core population was 49-fold, and the minimum modal read depth in any of the core population was 14-fold: 1,986,930 genotypes were computed in total.

A panel of 90 sequencing libraries from each of three well-characterised isolates (BPK282, BPK275 and BPK087) representing three major groups from our core population that were generated as part of this study were used to calibrate the SNP calling across the population. The metric of the mean number of callers per candidate SNP site was determined for the entire set of 300 libraries. Candidate SNP sites in lower complexity regions were accepted when the mean number of callers ≥3.5; those in non-masked regions were kept where the mean number of callers ≥2.5 (Table A1.1); those in hard-masked regions were omitted. These criteria resulted in a SNP set supported by a mean of 4.5 callers per SNP. An additional screening condition on the sum of the individual SNP caller score (>380) was applied to some lower quality runs (for isolates BHU569A1, BHU572A1, BHU574A1, BHU575A1, BHU581A1, BHU592A1, BHU741A1, BHU764A1, BHU770A1, BHU777A1, BHU782A1, BHU814A1, BPK031A1, BPK091A1, BPK157A1, BPK164A1, BPK406A1 and BPK512A1).

Table A1.1

Number of putative SNPs detected by different numbers of SNP calling tools. A threshold mean number of callers ≥3 was used, which was elevated to ≥4 at soft-masked regions.

https://doi.org/10.7554/eLife.12613.019
#SNP callersCountCumulative countFreqCumulative Freq
132,07232,0720.8220.822
1.415434,4750.0040.884
1.549834,9730.0130.896
296536,1600.0250.927
2.53836,2440.0010.929
347636,7570.0120.942
3.511436,9350.0030.947
41,01938,0460.0260.975
4.517938,3470.0050.983
543839,0150.0111

The variant read allele frequency, read depth and sum of weighted SNP scores were measured for each putative SNP to evaluate SNP quality. While false positive SNPs would be expected to occur more frequently in clusters due to local structural variations or repetitive sequences: mean distance between COCALL SNPs was approximately equal between regions of the genome. Allele frequencies inferred from our genotype calls were highly correlated for the 2418 heterozygous SNPs (r2=0.999) with those derived from frequencies of reads with each allele for the 191 core isolates. The heterozygous allele read depth distributions indicated some deviations from the expected value of 0.5 due to a mix of strand bias, chromosome copy number variation, and variation in read mapping and alignment. However, few SNPs were heterozygous and few genotypes were ambiguous: 14,825 mutations had read depth allele frequencies of 0.95 or greater.

4 Verification of variants called in individual libraries

We validated 15 true heterozygous SNPs in the reference BPK292/0cl4 isolate by manually inspecting mapped reads for this isolate, and confirmed that COCALL correctly called these sites. In addition, we examined variant calls on two isolates (BPK519/12 and BPK156/0; ISC1) where chromosome 17 is inferred to be monosomic. The mean derived allele frequency at variant sites on these chromosomes was 99.3% (at 1338 SNPs in BPK519/12) and 99.9% (at 1336 SNPs in BPK156/0). On this chromosome, COCALL called 1337 homozygous SNPs and a single heterozygous one for BPK519/12 and 1318 homozygous SNPs and 18 heterozygous ones for BPK156/0. While eliminating all heterozygous calls may seem optimal, it was possible that some cells were not monosomic for this chromosome, and choosing parameters that eliminated all heterozygous SNP calls removed many valid homozygous ones

5 Comparison of COCALL with GATK population calling

We tested the relative rate of true and false SNP sites and individual genotypes for COCALL and population-based GATK (DePristo et al., 2011) using read data of three previously studied isolates (Downing et al., 2011) (BPK282/0cl4, BPK087/0cl11, BPK275/0cl18, sampled from the most frequent genetic populations ISC4, ISC5 and ISC6 respectively). We created artificial pairwise mixes of these samples in proportions of 1:9, 2:8, 3:7, 4:6, 5:5, 6:4, 7:3, 8:2, 9:1 and unmixed, with depth levels adjusted for the relative proportions of unmapped reads in each. With these data, COCALL identified more rare variant SNPs in mixed libraries for the BPK275/0cl18-BPK087/0cl11 mix (Table A1.2). BPK282/0cl4 was the reference genome strain and so should possess only heterozygous SNPs, and so its variant discovery rate was substantially lower when combined in mixed libraries. Consequently, for the mixed libraries of BPK282/0cl4-BPK087/0cl11 and BPK282/0cl4-BPK275/0cl18, once the fraction of non-reference reads <30%, the sensitivity of COCALL was superior to that for the GATK population caller.

Table A1.2

The fraction of discordant true and false SNPs between COCALL and GATK population-based variant calling. COCALL had a greater proportion of true SNPs – confirmed by visual inspection of mapping reads – both in terms of the number of variant sites, and across the population of samples.

https://doi.org/10.7554/eLife.12613.020
SNP callerSitesMutations
TrueFalseTrueFalse
COCALL10633991118
76.3%23.7%89.4%10.6%
GATK population caller33106118991
23.7%76.3%10.6%89.4%

Specifically, the population-based GATK and COCALL both called 11,389 variant genotypes in the artficially mixed samples, while GATK called 718 SNPs not detected using COCALL, and COCALL found 396 additional SNPs not in the GATK set: a total of 1109 discordant variants at 139 sites. Visual inspection of mapped reads using Artemis (Carver et al., 2008) and Samtools tview (Li et al., 2009) indicated that 89% (354) of these 396 unique to COCALL were valid. COCALL had a greater proportion of true SNPs both in terms of the number of variant sites, and across the population of samples. This illustrated COCALL’s improved power for detecting variants at novel sites, compared with GATK (Table A1.3).

Table A1.3

The rate of SNP calling for different combinations of read libraries of samples BPK282/0cl4, BPK087/0cl11 and BPK275/0cl18. BPK282/0cl4 was the reference genome strain and so should possess only heterozygous SNPs, and so its variant discovery rate was substantially lower. COCALL identified more rare variant SNPs in mixed libraries.

https://doi.org/10.7554/eLife.12613.021
Read library combinationsUnique to COCALLSharedUnique to GATK population caller
BPK282/0cl4BPK275/0cl18BPK087/0cl11
89%11%141406
80%20%71827
70%30%61917
60%40%91937
50%50%71947
40%60%51966
30%70%71916
20%80%71788
11%89%91456
89%11%10505
80%20%111125
70%30%51294
60%40%71315
50%50%71316
40%60%61326
30%70%51336
20%80%41306
11%89%41276
89%11%15643
80%20%141074
70%30%41293
60%40%41323
50%50%31334
40%60%31325
30%70%41314
20%80%31305
11%89%31266

COCALL discovered a total of 2418 SNPs (4.3% of the original candidate set) in the 191 isolates from the core population. We calculated that our calling approach had a 99% chance of finding variants segregating at a frequency ≥0.024 in the Core 191 and ≥0.173 in ISC1 (Table A1.4).

Table A1.4

The probabilities of observing alleles at a specific frequency in the core population and ISC1. We assessed whether or not the SNPs we discovered in our sample of ISC L. donovani were likely to be representative of the variation existing in the wider population of parasites circulating in the ISC. We first assessed the chance of discovering alleles at frequency f as 1-(1-f)^2n in the core population (n=191) and ISC1 (n=12) sample sets (Gutenkunst et al., 2009). As a second approach, we evaluated the chance of sampling at least one derived allele in n sampled chromosomes (as above) assuming these come from a population of size N=1000 following a hypergeometric distribution (Tennessen et al., 2012). Both approaches confirm we have good power (>99% probability) of sampling alleles with population frequencies as low as 1.2% in the core population, but are likely to miss many variants in ISC1.

https://doi.org/10.7554/eLife.12613.022
Allele FrequencyP(derived allele)=1-(1-f)2n

M0 N-Mn /Nn
P(derived allele)=1 –
Core populationISC1Core populationISC1
0.0010.3180.0240.3820.024
0.0020.5350.0470.6180.047
0.0030.6830.070.910.07
0.0040.7840.0920.9450.092
0.0050.8530.1130.9660.115
0.0060.90.1340.9790.136
0.0070.9320.1550.9870.157
0.0080.9530.1750.9920.177
0.0090.9680.1950.9970.197
0.010.9780.2140.9990.217
0.0120.990.2520.9990.254
0.0150.9970.3040.9990.307
0.0210.9990.3990.9990.403
0.0240.9990.44210.446
0.0290.9990.50710.511
0.0380.9990.60510.61
0.0390.9990.61510.619
0.0410.62510.629
0.04910.70110.705
0.06510.80110.805
0.09210.90110.904
0.17310.9910.99
0.2721111

6 Singleton SNPs and transition-transversion ratios

Singleton SNP calls (those called in only a single isolate) are expected to be more error-prone than shared variants (Han et al., 2014). We examined the fractions of pairs of singletons; doubletons (SNPs observed in two samples); and all other SNPs (≥3 samples) that co-occur within blocks of a range of sizes (10 b, 100 b, 500 b, 1 kb, 5 kb, 10 kb, 50 kb, 100 kb, 500 kb, 1 Mb, 5 Mb) to evaluate evidence of clustering: erroneous regions tend to have SNPs clustered together, and so more pairs of SNPs will co-occur within small blocks. We see this effect in SNP calls from Samtools Pileup (Table A1.5) but not for other callers or in the COCALL consensus. No difference between rates of SNPs called at coding and non-coding regions were observed. Very few sites were called as having more than two alleles (Supplementary file 2, table t).

Table A1.5

The percentage of SNP pairs for all pairs on the same chromosome for which both SNPs occur within a block, for different block sizes (kb) for Cortex, GATK, Samtools Mpileup and Pileup (Supplementary Table 6). Singletons may have clustered for Samtools Pileup, but not for Cortex, GATK or Samtools Mpileup. The mean distance was calculated the average distance between SNPs of the given type on the same chromosome. Non-rare SNPs were those found in three or more samples.

https://doi.org/10.7554/eLife.12613.023

SNP Caller

SNP type

Block size (kb)

0.01

0.1

0.5

1

5

10

50

100

500

Cortex

Non-rare

0.0

0.3

0.9

1.1

1.7

3.4

9.6

17.0

61.5

Doubletons

1.5

1.5

1.5

1.5

3.1

3.1

10.8

13.8

64.6

Singletons

0.0

0.0

0.0

0.0

3.8

3.8

23.1

30.8

80.8

All

0.0

0.3

0.9

1.1

1.7

3.3

9.6

17.0

61.5

GATK

Non-rare

0.0

0.1

0.2

0.3

1.1

2.1

9.3

17.8

64.4

Doubletons

1.3

1.3

1.3

1.3

1.3

1.3

5.3

13.3

61.3

Singletons

0.0

0.0

0.0

0.0

1.1

3.2

11.8

18.3

77.4

All

0.0

0.1

0.2

0.3

1.1

2.1

9.3

17.8

64.4

Samtools Mpileup

Non-rare

0.0

0.1

0.3

0.4

1.2

2.2

9.5

17.9

64.4

Doubletons

0.0

0.0

0.0

0.0

6.7

6.7

6.7

20.0

60.0

Singletons

0.0

0.0

0.0

0.0

9.1

9.1

9.1

9.1

81.8

All

0.0

0.1

0.3

0.4

1.2

2.2

9.5

17.9

64.4

Samtools Pileup

Non-rare

0.0

0.1

0.2

0.3

1.1

2.1

9.3

17.8

64.4

Doubletons

2.6

3.3

3.3

3.3

5.3

6.6

11.8

21.1

63.2

Singletons

0.8

4.4

5.2

5.2

5.7

6.3

13.4

24.9

71.6

All

0.0

0.1

0.2

0.3

1.1

2.1

9.3

17.8

64.5

We also looked at the transition/transversion (Ti/Tv) ratio, as transitions are expected biologically to be more common than transversions, while random (noisy) SNP calls will have Ti/Tv close to 0.5, as there are 4 possible transition mutations and 8 transversions. We calculated Ti/Tv ratios for variants in the core population of 191 (plus one additional closely related sample, BHU1087/0). We find that SNP calls supported by multiple callers had higher Ti/Tv ratios (Table A1.6), with the highest ratio for SNPs called by all five callers (2.74). Ti/Tv ratios consistently differed between singletons, doubletons and non-rare variants regardless of SNP caller: for SNPs called by all five methods, non-rare variants had higher rates (Ti/Tv=7267/1914=3.80) than doubletons (558/298=1.87) or singletons (1316/1124=1.17), as expected if singleton variants are more likely to be false positives. Ti/Tv ratios of non-rare SNPs were much higher than singleton SNPs across SNP sets called by Cortex-Mpileup (3.61 vs 0.91), Cortex-Pileup (3.08 vs 0.79), Cortex-GATK (2.26 vs 0.53), and Cortex-Freebayes (1.51 vs 0.40). While non-protozoan inter-species comparisons find a Ti/Tv ratio of 2.0–2.1 for genome-wide datasets and 3.0–3.3 for coding SNPs (Ebersberger et al., 2002; Freudenberg-Hua et al., 2003), we find higher Ti/Tv ratios at non-coding regions than coding regions, irrespective of SNP caller(s) used. Moreover, this effect was higher for high-confidence SNPs such as those called by all five callers and those present in multiple isolates. For example, for SNPs called by all five calling tools, all variants: coding Ti/Tv=4154/1736=2.39 vs non-coding 4987/1600=3.12; non-rare SNPs: coding 3285/1033=3.18 vs non-coding 3982/881=4.52; singleton SNPs: coding 617/565=1.09 vs non-coding 699/559=1.15.

Table A1.6

The transition (Ti) and transversion (Tv) mutation rates and Ti/Tv ratio across five SNP callers (Cortex, Freebayes, GATK, Samtools Mpileup, Samtools Pileup). SNP sets were selected from the intersection of the SNP calling methods marked X on the corresponding columns (for example, the first row uses the intersection of SNPs called by all five methods). The independent genome assembly and haplotype-based variant calling scheme utilised by Cortex complemented other calling methods based on Smalt assemblies well, based on the increase of the Ti/Tv ratio.

https://doi.org/10.7554/eLife.12613.024
CortexFreebayesGATKSamtools mpileupSamtools pileupTiTvTi/Tv
XXXXX914133362.74
XXXX922934112.71
XXXX915034182.68
XXX923934992.64
XXXX932338322.43
XXX941839512.38
XXX933939602.36
XX943641272.29
XXXX919444262.08
XXX922145912.01
XXX928946282.01
XX934248411.93
XXX950768991.38
XX964275361.28
XX10,39913,1010.79
X13,35128,3880.47
X1,904,0043,340,4320.57
X152226553,8890.27
X41,107137,3150.30
X42,902153,2950.28

Appendix 2 - Somy, CNV, indel and episome detection

1 Chromosome copy number estimation

The reads mapped by Smalt v5.7 were used to estimate the read depth for each sample at each chromosome and region using the coverage values from Samtools v0.1.18. Chromosomal read depths were computed from the median read depth for sites with depths within one standard deviation of the initial estimate of the depth for each chromosome, and normalised as the depth per haploid genome as outlined previously (Downing et al., 2011). The genome-wide values per sample allowed the assignment of monosomic, disomic, trisomic, tetrasomic, pentasomic and hexasomic chromosomes as continuous values to allow for intra-sample variation in somy among cells. Where the most frequent somy was trisomy rather than disomy, then chromosome copy number was normalised to reflect that trisomy was the median chromosome read depth of a given strain. The range of monosomy, disomy, trisomy, tetrasomy and pentasomy was defined to be the haploid normalized chromosome depth dch of dch < 0.65, 0.65 ≤ dch < 1.25, 1.25 ≤ dch < 1.75, 1.75 ≤ dch < 2.25 and 2.25 ≤ dch < 2.75 (Table A2.1).

Table A2.1

Somy status of 36 chromosomes in 206 strains (204 with JPCM5 and LV9). The range of monosomy, disomy, trisomy, tetrasomy and pentasomy was defined to be the haploid normalized chromosome depth dch of dch < 0.65, 0.65 ≤ dch < 1.25, 1.25 ≤ dch < 1.75, 1.75 ≤ dch < 2.25 and 2.25 ≤ dch < 2.75.

https://doi.org/10.7554/eLife.12613.025
CountPercentage
Monosomy20.03
Disomy585178.90
Trisomy122916.57
Tetrasomy3104.18
Pentasomy240.32

2 Copy number variation

Baseline coverage levels were established using median depth of over 60 strains among ISC4, ISC6, ISC7 and ISC8 so that potential CNVs were determined in comparison to the haploid normalised depth. For CNV detection, we also considered CNV length along with depth because long heterozygous deletions or duplications were observed that could be reliably identified at lower depth than shorter variants. Therefore we used three different depth thresholds: for CNVs of 2 to 5 kb, the threshold was five times the standard deviation of the chromosomal depth; for CNVs of 5 to 20 kb, the threshold was three times this value; for CNVs of >20 kb, the threshold was 1.5 times this value.

Adjacent CNV pairs in individual samples were combined when two nearby CNVs were located within 1/3 of the length of the longer of the two CNVs, because two independent adjacent CNVs would be a less parsimonious pattern. When a CNV pairwise distance was calculated, common boundaries of a given CNV among the samples were assigned to include every start and end of the CNV of each strain, and then the depth was calculated for each CNV range. Boundaries of CNVs were manually checked and CNVs overlapping large gaps were excluded. Candidate CNVs in the first 300 bp and last 5 Kb of all chromosomes were excluded. Adjacent CNVs were merged if the inter-CNV gap was smaller than the average length of the flanking CNVs.

3 Insertion and deletion discovery

Small indels, less than 10 bp, were detected by a four-caller method (Freebayes, GATK, Samtools Mpileup, Cortex). Pileup was not used because its indel detection method was not sufficiently powerful. The indel calling scheme was analogous to COCALL except that indels identified by at least three methods were accepted (instead of 2.5). To count allele frequency of indels at AQP1 more precisely, we used Smalt v7.1 that can map reads around indels more accurately by locating indels at the same base position for both forward and reverse strands. The AQP1 indel at 7,739 on chromosome 31 was visually inspected for all the strains using IGV and Samtools tview.

4 MAPK1 and H-locus amplicon detection

Read depth coverage was used to identify potential extra- or intra-chromosomal gene amplifications where levels were significantly elevated over an extended 10–50 Kb regions flanked by repetitive segments. Two of these variants: the copy number variation of the MAPK1 and H-locus amplicons was verified in a qPCR assay (SensiMix SYBR No-ROX, Bioline) performed on the LightCycler 480 (Roche Life Science). The genes MRPA (forward 5’-TGTGTTTCCGACGATTGC-3’, reverse 5’-GTGACCCGCTTTGTGGAC-3’) and MAPK1 (forward 5’-GTGGTCGCGCTGCAGAAG-3’, reverse 5’-CGGCACAACCCCTTCATTG-3’) were targeted in a subset of 46 samples. Cysteine synthase (LdBPK_363750) was used as a reference gene (forward 5’-GTCTTGGCGGTTCAGTTCG-3’, reverse 5’- GACATTGTGGTTCGTCTGCTC-3’), and BPK026/0cl5 was used for normalization of copy numbers. This strain has three and two chromosomal copies of the H-locus and MAPK1, respectively.

The nature of the amplification (extra- or intra-chromosomal) was determined by pulsed-field gel electrophoresis (PFGE) and southern blot hybridization. Intact chromosomes were prepared from logarithmic phase cultures of Leishmania promastigotes and separated by PFGE using a Bio-Rad CHEF-Mapper XA at 6V/cm, 120° separation angle and a range of separation from 100 kb to 2.5 Mb. Gels were transferred by capillarity and hybridized with [α-32P]dCTP-labelled DNA probes specific for MAPK1, HTBF (H-locus) or LinJ35.4130 (LdBPK_354130), according to standard protocols.

To exclude the possibility that the amplicons are a culturing artefact, a PCR using primers that enable amplification on circular episomes or tandem duplications was also attempted directly on five bone marrow samples from VL patients. The edges of the chromosomal loci that code for the MAPK1 and DHFR-TS genes are targeted with primers (MAPK1: forward 5’-CATGGCGCAGTGACCTTCAG-3’, reverse 5’-TCTTGGCACGGCATCAGCAG-3’; DHFR-TS: forward 5’-CCTACCCGCTTGCTTGCTTG-3’, reverse 5’- CAGCAGCACAATGGAAAGAACG-3’).

5 Correlation between SNPs, CNVs, indels and somy variation

The correlation coefficients of pairwise distances between samples based on SNP, CNV, indel and somy diversity were calculated. Some strains were clearly phylogenetically related, so correlation between these values was expected simply due to the phylogenetic structure of our samples. To alleviate this potential bias, we used Pearson’s correlation for sets of samples with different degrees of similarity based on SNP pairwise distance. The correlation coefficients of SNPs with indels and CNVs with SNPs were consistently large in all sample sets (Table A2.2).

Table A2.2

Pearson’s correlation (r2) between SNP, CNV, indels and somy diversity. The dependency of the correlation on the normalised pairwise distance is shown: this was determined from the median normalised pairwise distance based on SNP data. The Mantel correlation tests showed values that were consistent with the correlations.

https://doi.org/10.7554/eLife.12613.026
Type#samplesNormalised pairwise distance thresholdCorrelation (r2)P value
SNP-INDEL19100.4330.00E+00
570.150.3012.91E-254
330.30.2718.23E-77
191Mantel test0.4200.001
PLO-SNP19100.0730.00E+00
570.150.1522.69E-118
330.30.3049.62E-88
191Mantel test0.0600.001
PLO-INDEL19100.0243.88E-197
570.150.0402.43E-30
330.30.0721.96E-19
191Mantel test0.0160.001
CNV-SNP19100.2890.00E+00
570.150.2637.53E-218
330.30.3934.52E-120
191Mantel test0.2720.001
CNV-PLO19100.1210.00E+00
570.150.1501.96E-116
330.30.1831.59E-49
191Mantel test0.1030.001
CNV-INDEL19100.0990.00E+00
570.150.0872.34E-66
330.30.1081.04E-28
191Mantel test0.0840.001

Correlations of our values with somy were more variable, probably because somy was more prone to homoplasy as it appears to mutate faster than other genetic variants in this study (Table A2.2). Spearman's rank correlations were also calculated and showed a similar level of correlation (data not shown). We also used Mantel tests to directly test correlation between the full matrices of pairwise distances. The results of Mantel tests were in concordance with Pearson’s correlation values. Pearson’s and Spearman's correlations were calculated using python-stats-0.6 (Gary Strangman http://sourcecodebrowser.com/python-stats/0.6/stats_8py.html). Mantel tests were implemented in R package vegan (https://github.com/vegandevs/vegan).

6 Correlation between duplicated genes and somy variation

Information on duplicated genes in the reference genome was based on multi-copy gene families annotated in the L. infantum JPCM5 reference genome (Table A2.3) (Rogers et al., 2011). In the present study, each read was mapped to its most similar region in the reference, and where reads could be mapped with equal mapping scores to multiple locations, one position was selected randomly to quantify read depth, in contrast to the previous study in which read depth was under-estimated by only using uniquely-mapped reads (Rogers et al., 2011). A threshold of 1.5 times the haploid depth was used to identify duplicated genes (a more conservative threshold of 2.0 did not change results materially). This threshold of 1.5 was effective for Illumina Hiseq data, where fluctuations in read depth were smaller than those previously seen using the Illumina Genome Analyser II. To illustrate this, the normalised standard deviation (standard deviation divided by average chromosome depth) of 17 strains sequenced by Genome Analyser II in was 0.32 (Downing et al., 2011), and in the current analysis, using Illumina Hiseq 2000 data, this was 0.19. We were therefore able to quantify more subtle depth variation than previous analyses. To normalise for differences in chromosome lengths, the number of duplicated genes was expressed per unit chromosome with a length set to that of chromosome 1. Results from between-species analysis (Rogers et al., 2011) have previously suggested a link between somy variation and copy number variation. However, in the present study, the correlation of duplicated genes per chromosome with the somy level was statistically significant but the effect size was very small (n=7146 chromosomes, r2=0.027, p=7.03x10-46), suggesting no relationship or only a weak one, between somy level and copy number at the population level.

Table A2.3

Correlation between CNV and ploidy of chromosomes in 206 strains: the duplication threshold was 0.5. p values <0.02 were observed for chromosomes 1, 2, 3, 17, 19, 20, 21 and the whole genome.

https://doi.org/10.7554/eLife.12613.027
Correlation (r2)P value
chromosome0.0277.03E-46
10.0394.23E-03
20.2021.21E-11
30.0747.81E-05
40.0034.31E-01
509.30E-01
60.0043.86E-01
70.0034.68E-01
80.0016.05E-01
90.0017.32E-01
100.0157.77E-02
110.0024.88E-01
120.0053.16E-01
130.0016.57E-01
140.0017.31E-01
150.0017.42E-01
160.0091.81E-01
170.0366.06E-03
180.0053.23E-01
190.0775.41E-05
200.0642.51E-04
210.0271.87E-02
220.0091.87E-01
230.0072.42E-01
2409.68E-01
250.0091.87E-01
260.0175.97E-02
270.0091.85E-01
280.0016.80E-01
290.0091.67E-01
300.0034.28E-01
310.0043.77E-01
320.0053.28E-01
330.0111.41E-01
340.0017.61E-01
350.0111.35E-01
360.0167.38E-02

Data availability

The following data sets were generated

References

    1. Busing FMTA
    2. Meijer E
    3. Van Der Leiden R
    (1999)
    Delete-M Jackknife for Unequal M 
    Statistics and Computing 9:3–8.
    1. Dye C
    2. Wolpert DM
    (1988)
    Earthquakes, influenza and cycles of Indian kala-azar
    Transactions of the Royal Society of Tropical Medicine and Hygiene 82:843–850.
    1. Garrison E
    2. Marth G
    (2012)
    Haplotype-based variant detection from short-read sequencing
    arXiv 1207.3907.
    1. Gibson ME
    (1983)
    The identification of kala-azar and the discovery of Leishmania donovani
    Medical History 27:203–213.
    1. Hudson RR
    2. Kaplan NL
    (1985)
    Statistical properties of the number of recombination events in the history of a sample of DNA sequences
    Genetics 111:147–164.
    1. Li N
    2. Stephens M
    (2003)
    Modeling linkage disequilibrium and identifying recombination hotspots using single-nucleotide polymorphism data
    Genetics 165:2213–2233.
    1. Lozano R
    2. Naghavi M
    3. Foreman K
    4. Lim S
    5. Shibuya K
    6. Aboyans V
    7. Abraham J
    8. Adair T
    9. Aggarwal R
    10. Lim S
    11. AlMazroa MA
    12. Alvarado M
    13. Anderson HR
    14. Anderson LM
    15. Andrews KG
    16. Atkinson C
    17. Baddour LM
    18. Barker-Collo S
    19. Bartels DH
    20. Bell ML
    21. Benjamin EJ
    22. Bennett D
    23. Bhalla K
    24. Bikbov B
    25. Abdulhak AB
    26. Birbeck G
    27. Blyth F
    28. Bolliger I
    29. Lim S
    30. Bucello C
    31. Lim S
    32. Burney P
    33. Carapetis J
    34. Chen H
    35. Chou D
    36. Lim S
    37. Coffeng LE
    38. Lim S
    39. Colquhoun S
    40. Colson KE
    41. Condon J
    42. Connor MD
    43. Cooper LT
    44. Corriere M
    45. Cortinovis M
    46. de Vaccaro KC
    47. Couser W
    48. Cowie BC
    49. Criqui MH
    50. Lim S
    51. Dabhadkar KC
    52. Dahodwala N
    53. De Leo D
    54. Degenhardt L
    55. Delossantos A
    56. Denenberg J
    57. Des Jarlais DC
    58. Dharmaratne SD
    59. Dorsey ER
    60. Driscoll T
    61. Duber H
    62. Ebel B
    63. Erwin PJ
    64. Espindola P
    65. Ezzati M
    66. Feigin V
    67. Foreman K
    68. Forouzanfar MH
    69. Fowkes FGR
    70. Franklin R
    71. Fransen M
    72. Foreman K
    73. Gabriel SE
    74. Gakidou E
    75. Gaspari F
    76. Gillum RF
    77. Gonzalez-Medina D
    78. Halasa YA
    79. Haring D
    80. Harrison JE
    81. Havmoeller R
    82. Hay RJ
    83. Hoen B
    84. Hotez PJ
    85. Hoy D
    86. Jacobsen KH
    87. James SL
    88. Jasrasaria R
    89. Jayaraman S
    90. Johns N
    91. Karthikeyan G
    92. Kassebaum N
    93. Foreman K
    94. Khoo J-P
    95. Knowlton LM
    96. Kobusingye O
    97. Koranteng A
    98. Krishnamurthi R
    99. Lipnick M
    100. Lipshultz SE
    101. Ohno SL
    102. Mabweijano J
    103. MacIntyre MF
    104. Mallinger L
    105. March L
    106. Marks GB
    107. Marks R
    108. Matsumori A
    109. Matzopoulos R
    110. Mayosi BM
    111. McAnulty JH
    112. McDermott MM
    113. McGrath J
    114. Memish ZA
    115. Memish ZA
    116. Merriman TR
    117. Michaud C
    118. Miller M
    119. Miller TR
    120. Mock C
    121. Mocumbi AO
    122. Mokdad AA
    123. Foreman K
    124. Mulholland K
    125. Nair MN
    126. Naldi L
    127. Narayan KMV
    128. Foreman K
    129. Foreman K
    130. O'Donnell M
    131. Omer SB
    132. Foreman K
    133. Osborne R
    134. Ozgediz D
    135. Pahari B
    136. Pandian JD
    137. Rivero AP
    138. Padilla RP
    139. Perez-Ruiz F
    140. Perico N
    141. Phillips D
    142. Pierce K
    143. Pope CA
    144. Porrini E
    145. Pourmalek F
    146. Raju M
    147. Ranganathan D
    148. Lim S
    149. Rein DB
    150. Remuzzi G
    151. Rivara FP
    152. Roberts T
    153. De León FR
    154. Rosenfeld LC
    155. Rushton L
    156. Sacco RL
    157. Salomon JA
    158. Sampson U
    159. Foreman K
    160. Schwebel DC
    161. Segui-Gomez M
    162. Shepard DS
    163. Singh D
    164. Singleton J
    165. Sliwa K
    166. Smith E
    167. Steer A
    168. Taylor JA
    169. Thomas B
    170. Tleyjeh IM
    171. Towbin JA
    172. Truelsen T
    173. Undurraga EA
    174. Venketasubramanian N
    175. Vijayakumar L
    176. Lim S
    177. Wagner GR
    178. Naghavi M
    179. Wang W
    180. Watt K
    181. Weinstock MA
    182. Weintraub R
    183. Wilkinson JD
    184. Woolf AD
    185. Wulf S
    186. Yeh P-H
    187. Yip P
    188. Zabetian A
    189. Zheng Z-J
    190. Lopez AD
    191. Murray CJL
    192. Shibuya K
    193. Aboyans V
    194. Memish ZA
    (2012) Global and regional mortality from 235 causes of death for 20 age groups in 1990 and 2010: a systematic analysis for the Global Burden of Disease Study 2010
    The Lancet 380:2095–2128.
    https://doi.org/10.1016/S0140-6736(12)61728-0
    1. Paradis E
    2. Claude J
    3. Strimmer K
    (2004)
    APE: Analyses of Phylogenetics and Evolution in R language
    Bioinformatics 20:289–290.
    1. Pritchard JK
    2. Stephens M
    3. Donnelly P
    (2000)
    Inference of population structure using multilocus genotype data
    Genetics 155:945–959.
    1. Thakur CP
    (2007)
    A new strategy for elimination of kala-azar from rural Bihar
    The Indian Journal of Medical Research 126:447–451.
    1. Van Meirvenne N
    2. Janssens PG
    3. Magnus E
    (1975)
    Antigenic variation in syringe passaged populations of Trypanosoma (Trypanozoon) brucei. 1. Rationalization of the experimental approach
    Annales De La Société Belge De Médecine Tropicale 55:1–23.

Article and author information

Author details

  1. Hideo Imamura

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    HI, Developed and verified SNP and indel detection scheme COCALL, Aligned reads and determined read depth and performed genetic variant detection and correlation analysis among genetic variants, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Tim Downing and Frederik Van den Broeck
    Competing interests
    The authors declare that no competing interests exist.

  2. Tim Downing

    1. Wellcome Trust Sanger Institute, Hinxton, United Kingdom
    2. School of Maths, Applied Maths and Statistics, National University of Ireland Galway, Galway, Ireland
    Present address
    School of Biotechnology, Dublin City University, Dublin, Ireland
    Contribution
    TD, Developed and verified SNP and indel detection scheme COCALL, Analysed population structure and performed phylogenetic analyses, Examined the correspondence of mutation with drug phenotype variation, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Hideo Imamura and Frederik Van den Broeck
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8385-6730

  3. Frederik Van den Broeck

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    FVdB, Analysed correlations of aneuploidy with heterozygosity, Examined the correspondence of mutation with drug phenotype variation, Performed additional sequence data analyses, Analysis and interpretation of data, Drafting or revising the article
    Contributed equally with
    Hideo Imamura and Tim Downing
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2542-5585

  4. Mandy J Sanders

    Wellcome Trust Sanger Institute, Hinxton, United Kingdom
    Contribution
    MJS, Coordinated sequence materials and sequencing. Reviewed the final manuscript., Acquisition of data
    Competing interests
    The authors declare that no competing interests exist.

  5. Suman Rijal

    BP Koirala Institute of Health Sciences, Dharan, Nepal
    Contribution
    SRi, Coordinated sample collection and provided and documented samples and clinical data, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  6. Shyam Sundar

    Department of Medicine, Institute of Medical Sciences, Banaras Hindu University, Varanasi, India
    Contribution
    SSu, Coordinated sample collection and provided and documented samples and clinical data, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  7. An Mannaert

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    AM, Performed molecular-clock and phylogeographic analyses and assessed epidemiological and phylogenetic history, Performed additional sequence data analysis and additional experimental analysis and verification, Reviewed the final manuscript., Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.

  8. Manu Vanaerschot

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Present address
    Department of Microbiology and Immunology, Columbia University College of Physicians and Surgeons, New York, United States
    Contribution
    MV, Coordinated sample collection and provided and documented samples and clinical data, Coordinated sequence materials and sequencing, Performed additional sequence data analysis, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  9. Maya Berg

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    MBe, Performed additional sequence data analysis, Reviewed the final manuscript, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.

  10. Géraldine De Muylder

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    GDM, Performed amplicon verification, Performed AQP1 modelling and analysis, Performed additional experimental analysis and verification, Reviewed the final manuscript., Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.

  11. Franck Dumetz

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    FD, Performed additional experimental analysis and verification, Reviewed the final manuscript, Acquisition of data
    Competing interests
    The authors declare that no competing interests exist.

  12. Bart Cuypers

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    BC, Performed additional experimental analysis and verification, Reviewed the final manuscript, Acquisition of data
    Competing interests
    The authors declare that no competing interests exist.

  13. Ilse Maes

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    IM, Conducted drug susceptibility assays, Coordinated sequence materials and sequencing, Performed additional experimental analysis and verification, Reviewed the final manuscript., Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  14. Malgorzata Domagalska

    Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    MD, Performed analyses on aneuploidy. Reviewed the final manuscript
    Competing interests
    The authors declare that no competing interests exist.

  15. Saskia Decuypere

    1. Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    2. Telethon Kids Institute, University of Western Australia, Perth, Australia
    Contribution
    SD, Co-ordinated sample collection and provided and documented samples and clinical data, Co-ordinated sequence materials and sequencing. Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  16. Keshav Rai

    BP Koirala Institute of Health Sciences, Dharan, Nepal
    Present address
    West Bengal State University, Kolkata, India
    Contribution
    KR, Conducted drug susceptibility assays, Performed additional sequence data analysis, Performed additional experimental analysis and verification, Reviewed the final manuscript, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9747-3431

  17. Surendra Uranw

    BP Koirala Institute of Health Sciences, Dharan, Nepal
    Contribution
    SU, Coordinated sample collection and provided and documented samples and clinical data, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  18. Narayan Raj Bhattarai

    BP Koirala Institute of Health Sciences, Dharan, Nepal
    Contribution
    NRB, Conducted drug susceptibility assays, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  19. Basudha Khanal

    BP Koirala Institute of Health Sciences, Dharan, Nepal
    Contribution
    BK, Co-ordinated sample collection and provided and documented samples and clinical data, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  20. Vijay Kumar Prajapati

    Department of Medicine, Institute of Medical Sciences, Banaras Hindu University, Varanasi, India
    Present address
    Department of Biochemistry, Central University of Rajasthan, Ajmer, India
    Contribution
    VKP, Coordinated sample collection and provided and documented samples and clinical data, Conducted drug susceptibility assays, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  21. Smriti Sharma

    Department of Medicine, Institute of Medical Sciences, Banaras Hindu University, Varanasi, India
    Contribution
    SSh, Co-ordinated sample collection and provided and documented samples and clinical data, reviewed the final manuscript, acquired data, contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  22. Olivia Stark

    Institut für Mikrobiologie und Hygiene, Charité Universitätsmedizin Berlin, Berlin, Germany
    Contribution
    OS, Co-ordinated sample collection and provided and documented samples and clinical data, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  23. Gabriele Schönian

    Institut für Mikrobiologie und Hygiene, Charité Universitätsmedizin Berlin, Berlin, Germany
    Contribution
    GS, Co-ordinated sample collection and provided and documented samples and clinical data, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  24. Harry P De Koning

    Institute of Infection, Immunity and Inflammation, College of Medical, Veterinary and Life Sciences, University of Glasgow, Glasgow, United Kingdom
    Contribution
    HPDK, Performed AQP1 modelling and analysis, Reviewed the final manuscript, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9963-1827

  25. Luca Settimo

    1. Institute of Infection, Immunity and Inflammation, College of Medical, Veterinary and Life Sciences, University of Glasgow, Glasgow, United Kingdom
    2. Department of Chemistry and Chemical Biology, Northeastern University, Boston, United States
    Contribution
    LS, Performed AQP1 modelling and analysis, Reviewed the final manuscript, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.

  26. Benoit Vanhollebeke

    Laboratory of Molecular Parasitology, Université Libre de Bruxelles, Gosselies, Belgium
    Contribution
    BV, Performed amplicon verification, Performed additional experimental analysis and verification, Reviewed the final manuscript, Acquisition of data, Analysis and interpretation of data
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0353-365X

  27. Syamal Roy

    Department of Infectious Diseases and Immunology, Council of Scientific and Industrial Research, Indian Institute of Chemical Biology, Kolkata, India
    Contribution
    SRo, Co-ordinated sample collection and provided and documented samples and clinical data, Conducted drug susceptibility assays, Performed additional experimental analysis and verification, Reviewed the final manuscript, Acquisition of data, Analysis and interpretation of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  28. Bart Ostyn

    Department of Public Health, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    BO, Co-ordinated sample collection and provided and documented samples and clinical data, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  29. Marleen Boelaert

    Department of Public Health, Institute of Tropical Medicine, Antwerp, Belgium
    Contribution
    MBo, Co-ordinated sample collection and provided and documented samples and clinical data, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8051-6776

  30. Louis Maes

    Department of Biomedical Sciences, Faculty of Pharmaceutical, Biomedical and Veterinary Sciences, University of Antwerp, Antwerp, Belgium
    Contribution
    LM, Conducted drug susceptibility assays, Reviewed the final manuscript, Acquisition of data, Contributed unpublished essential data or reagents
    Competing interests
    The authors declare that no competing interests exist.

  31. Matthew Berriman

    Wellcome Trust Sanger Institute, Hinxton, United Kingdom
    Contribution
    MB, Conception and design, Drafting or revising the article
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9581-0377

  32. Jean-Claude Dujardin

    1. Department of Biomedical Sciences, Institute of Tropical Medicine, Antwerp, Belgium
    2. Department of Biomedical Sciences, Faculty of Pharmaceutical, Biomedical and Veterinary Sciences, University of Antwerp, Antwerp, Belgium
    Contribution
    J-CD, Co-ordinated sample collection and provided and documented samples and clinical data, Performed molecular-clock and phylogeographic analyses and assessed epidemiological and phylogenetic history, Performed additional sequence data analysis, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article, Contributed unpublished essential data or reagents
    For correspondence
    JCDujardin@itg.be
    Competing interests
    The authors declare that no competing interests exist.

  33. James A Cotton

    Wellcome Trust Sanger Institute, Hinxton, United Kingdom
    Contribution
    JAC, Analysed population structure and performed phylogenetic analyses, Performed molecular-clock and phylogeographic analyses and assessed epidemiological and phylogenetic history, Performed additional sequence data analysis, Analysis and interpretation of data, Drafting or revising the article
    For correspondence
    james.cotton@sanger.ac.uk
    Competing interests
    The authors declare that no competing interests exist.
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5475-3583

Funding

European Commission (EU framework program FP7-222895)

  • Hideo Imamura
  • Tim Downing
  • Frederik Van den Broeck
  • Mandy J Sanders
  • Suman Rijal
  • Shyam Sundar
  • An Mannaert
  • Manu Vanaerschot
  • Maya Berg
  • Géraldine De Muylder
  • Franck Dumetz
  • Ilse Maes
  • Saskia Decuypere
  • Keshav Rai
  • Surendra Uranw
  • Narayan Raj Bhattarai
  • Basudha Khanal
  • Vijay Kumar Prajapati
  • Smriti Sharma
  • Olivia Stark
  • Gabriele Schönian
  • Syamal Roy
  • Bart Ostyn
  • Marleen Boelaert
  • Matthew Berriman
  • Jean-Claude Dujardin
  • James A Cotton

Belgian Science Policy Office (TRIT, P7/41)

  • Hideo Imamura
  • Frederik Van den Broeck
  • An Mannaert
  • Manu Vanaerschot
  • Maya Berg
  • Géraldine De Muylder
  • Franck Dumetz
  • Ilse Maes
  • Benoit Vanhollebeke
  • Jean-Claude Dujardin

Flemish Fund for Scientific Research (G.0.B81.12)

  • Hideo Imamura
  • Frederik Van den Broeck
  • An Mannaert
  • Manu Vanaerschot
  • Maya Berg
  • Géraldine De Muylder
  • Franck Dumetz
  • Ilse Maes
  • Jean-Claude Dujardin

Department of Economy, Science and Innovation in Flanders (ITM-SOFIB)

  • Hideo Imamura
  • Manu Vanaerschot
  • Maya Berg
  • Saskia Decuypere
  • Jean-Claude Dujardin
  • Malgorzata Domagalska

Wellcome Trust (098051)

  • Mandy J Sanders
  • Matthew Berriman
  • James A Cotton

Belgian Development Cooperation (FA3 II VL control and FA3 project 95502)

  • Suman Rijal
  • Keshav Rai
  • Surendra Uranw
  • Narayan Raj Bhattarai
  • Basudha Khanal
  • Bart Ostyn
  • Marleen Boelaert
  • Jean-Claude Dujardin

INBEV-Baillet Latour Foundation

  • Maya Berg

Flemish Fund for Scientific Research (11O1614N)

  • Bart Cuypers

JC Bose National Fellowship, DST, Government of India (SB/S2/JCB-65/2014)

  • Syamal Roy

Council of Scientific and Industrial Research (BSC0114)

  • Syamal Roy

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors gratefully acknowledge financial support for the Kaladrug-R consortium from the European Union framework program (FP7-222895). JCD, BO, MBo are supported by the Belgian Development Cooperation (FA3 II VL control and FA3 project 95502), HI, FVdB, AM, MV, MBe, GD, FD, IM, and JCD by the Belgian Science Policy Office (TRIT, P7/41) and the Flemish Fund for Scientific Research (G.0.B81.12) and MBe by the INBEV-Baillet Latour foundation. JAC, MB, MJS are supported by the Wellcome Trust through their core funding of the Wellcome Trust Sanger Institute (grant 098051) and thank members of the DNA pipelines team at WTSI for generating the sequencing libraries.

Copyright

© 2016, Imamura et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 4,948
    views
  • 755
    downloads
  • 134
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Hideo Imamura
  2. Tim Downing
  3. Frederik Van den Broeck
  4. Mandy J Sanders
  5. Suman Rijal
  6. Shyam Sundar
  7. An Mannaert
  8. Manu Vanaerschot
  9. Maya Berg
  10. Géraldine De Muylder
  11. Franck Dumetz
  12. Bart Cuypers
  13. Ilse Maes
  14. Malgorzata Domagalska
  15. Saskia Decuypere
  16. Keshav Rai
  17. Surendra Uranw
  18. Narayan Raj Bhattarai
  19. Basudha Khanal
  20. Vijay Kumar Prajapati
  21. Smriti Sharma
  22. Olivia Stark
  23. Gabriele Schönian
  24. Harry P De Koning
  25. Luca Settimo
  26. Benoit Vanhollebeke
  27. Syamal Roy
  28. Bart Ostyn
  29. Marleen Boelaert
  30. Louis Maes
  31. Matthew Berriman
  32. Jean-Claude Dujardin
  33. James A Cotton
(2016)
Evolutionary genomics of epidemic visceral leishmaniasis in the Indian subcontinent
eLife 5:e12613.
https://doi.org/10.7554/eLife.12613

Share this article

https://doi.org/10.7554/eLife.12613

Categories and tags

Further reading

    1. Medicine
    2. Microbiology and Infectious Disease
    3. Epidemiology and Global Health
    4. Immunology and Inflammation

    Infectious Diseases: A Collection of Translational and Clinical Research Articles

    Edited by Jos WM van der Meer et al.
    Collection

    eLife has published articles on a wide range of infectious diseases, including COVID-19, influenza, tuberculosis, HIV/AIDS, malaria and typhoid fever.

    1. Epidemiology and Global Health
    2. Microbiology and Infectious Disease

    Livestock abortion surveillance in Tanzania reveals disease priorities and importance of timely collection of vaginal swab samples for attribution

    Felix Lankester, Tito J Kibona ... Sarah Cleaveland
    Research Article

    Lack of data on the aetiology of livestock diseases constrains effective interventions to improve livelihoods, food security and public health. Livestock abortion is an important disease syndrome affecting productivity and public health. Several pathogens are associated with livestock abortions but across Africa surveillance data rarely include information from abortions, little is known about aetiology and impacts, and data are not available to inform interventions. This paper describes outcomes from a surveillance platform established in Tanzania spanning pastoral, agropastoral and smallholder systems to investigate causes and impacts of livestock abortion. Abortion events were reported by farmers to livestock field officers (LFO) and on to investigation teams. Events were included if the research team or LFO could attend within 72 hr. If so, samples and questionnaire data were collected to investigate (a) determinants of attribution; (b) patterns of events, including species and breed, previous abortion history, and seasonality; (c) determinants of reporting, investigation and attribution; (d) cases involving zoonotic pathogens. Between 2017–2019, 215 events in cattle (n=71), sheep (n=44), and goats (n=100) were investigated. Attribution, achieved for 19.5% of cases, was significantly affected by delays in obtaining samples. Histopathology proved less useful than PCR due to rapid deterioration of samples. Vaginal swabs provided practical and sensitive material for pathogen detection. Livestock abortion surveillance, even at a small scale, can generate valuable information on causes of disease outbreaks, reproductive losses and can identify pathogens not easily captured through other forms of livestock disease surveillance. This study demonstrated the feasibility of establishing a surveillance system, achieved through engagement of community-based field officers, establishment of practical sample collection and application of molecular diagnostic platforms.

    1. Epidemiology and Global Health
    2. Genetics and Genomics

    Systematic evaluation of multifactorial causal associations for Alzheimer’s disease and an interactive platform MRAD developed based on Mendelian randomization analysis

    Tianyu Zhao, Hui Li ... Li Chen
    Research Article

    Alzheimer’s disease (AD) is a complex degenerative disease of the central nervous system, and elucidating its pathogenesis remains challenging. In this study, we used the inverse-variance weighted (IVW) model as the major analysis method to perform hypothesis-free Mendelian randomization (MR) analysis on the data from MRC IEU OpenGWAS (18,097 exposure traits and 16 AD outcome traits), and conducted sensitivity analysis with six models, to assess the robustness of the IVW results, to identify various classes of risk or protective factors for AD, early-onset AD, and late-onset AD. We generated 400,274 data entries in total, among which the major analysis method of the IVW model consists of 73,129 records with 4840 exposure traits, which fall into 10 categories: Disease, Medical laboratory science, Imaging, Anthropometric, Treatment, Molecular trait, Gut microbiota, Past history, Family history, and Lifestyle trait. More importantly, a freely accessed online platform called MRAD (https://gwasmrad.com/mrad/) has been developed using the Shiny package with MR analysis results. Additionally, novel potential AD therapeutic targets (CD33, TBCA, VPS29, GNAI3, PSME1) are identified, among which CD33 was positively associated with the main outcome traits of AD, as well as with both EOAD and LOAD. TBCA and VPS29 were negatively associated with the main outcome traits of AD, as well as with both EOAD and LOAD. GNAI3 and PSME1 were negatively associated with the main outcome traits of AD, as well as with LOAD, but had no significant causal association with EOAD. The findings of our research advance our understanding of the etiology of AD.