Drug resistance analysis by next generation sequencing in Leishmania.

The use of next generation sequencing has the power to expedite the identification of drug resistance determinants and biomarkers and was applied successfully to drug resistance studies in Leishmania. This allowed the identification of modulation in gene expression, gene dosage alterations, changes in chromosome copy numbers and single nucleotide polymorphisms that correlated with resistance in Leishmania strains derived from the laboratory and from the field. An impressive heterogeneity at the population level was also observed, individual clones within populations often differing in both genotypes and phenotypes, hence complicating the elucidation of resistance mechanisms. This review summarizes the most recent highlights that whole genome sequencing brought to our understanding of Leishmania drug resistance and likely new directions.

The parasite Leishmania and chemotherapy

Leishmania are digenetic parasites that develop as promastigotes in the gut of phlebotomine sandflies and as intracellular amastigotes in the macrophages of vertebrate hosts. The Leishmania genus gathers together several species and can be further subdivided into two subgenera (the subgenus Leishmania Viannia and the subgenus Leishmania Leishmania) based on the localization of promastigotes in the alimentary tract of the insect. The species L. donovani and L. infantum (also referred to as L. chagasi in South America) are mainly responsible for visceral leishmaniasis, a condition that manifests as a gross inflammatory reaction within the viscera (spleen and liver) that is fatal if left untreated. On the other hand, the etiology of cutaneous leishmaniasis is much more diverse and manifests as self-healing skin ulcers or nodules that remain localized at the site of inoculation or metastasize to remote body sites (Reithinger et al., 2007). A more severe complication called muco-cutaneous leishmaniasis involving the destruction of the nasal mucosa can also happen and is most commonly associated with L. Viannia braziliensis in South America. Globally, 1–1.5 million new cases of leishmaniasis are estimated to occur each year with a case-fatality rate of 10% for the visceral form of the disease (Alvar et al., 2012).

There is not yet an effective vaccine registered against Leishmania but the immune protection observed upon healing of cutaneous leishmaniasis suggests that designing an effective vaccine should be achievable (reviewed in (Palatnik-de-Sousa, 2008)). Meanwhile, the control of leishmaniasis primarily relies on chemotherapy with only a limited number of registered molecules available. Pentavalent antimonials (either as sodium stibogluconate, meglumine antimoniate or generic formulations) have been the standard drug for more than 60 years and remain the primary line of treatment in many endemic regions, apart from Nepal and the Bihar state of India where antimonial formulations have been rendered almost obsolete due to widespread parasite resistance (Lira et al., 1999; Sundar et al., 2000; Rijal et al., 2003). Other recommended therapies include the polyene antibiotic amphotericin B for which a single-dose was shown to be 95% effective against visceral leishmaniasis in India (Sundar et al., 2010). Liposomal amphotericin B has become a standard treatment in many countries (Bern et al., 2006) but remains expensive even for single-course treatments (Meheus et al., 2010) and requires administration by intravenous route. There have been extensive efforts to develop new AmB formulations to replace the costly lipid-based formulations but no nonlipid-based AmB delivery systems have yet reached the clinic (Mohamed-Ahmed et al., 2012). Geographical differences in response rates to liposomal amphotericin B were also reported, with visceral leishmaniasis cases in India being more responsive than those from East Africa or South America (Berman et al., 1998). The alkyl-lysophospholipid analog miltefosine, a drug initially developed as an anti-tumoral compound, is the first effective oral drug against Leishmania (Croft et al., 1987; Jha et al., 1999). It has been successfully used for the treatment of visceral leishmaniasis since its registration in 2002 in India (Sundar et al., 2002) and was incorporated into the visceral leishmaniasis elimination programme for the Indian sub-continent. However, relapse rates of 20% recently observed in Nepal at 12 months post-treatment constitute an alarming signal and emphasize the need for careful monitoring (Rijal et al., 2013). Miltefosine might also be effective against cutaneous leishmaniasis, although regional differences in susceptibility were reported in South America (Soto et al., 2001, 2004, 2008; Gonzalez et al., 2009). The aminoglycoside paromomycin is the molecule the most recently approved for the treatment of visceral leishmaniasis. It was shown to have anti-leishmanial properties more than fifty years ago but it is only in the late 1990s that the efficacy of a parenteral formulation against visceral leishmaniasis was demonstrated by a phase III clinical trial, leading to registration for treatment of visceral leishmaniasis in India in 2006 (Jha et al., 1998; Sundar et al., 2007). Similarly to amphotericin B and miltefosine however, geographical variations also occur for the efficacy of paromomycin against visceral leishmaniasis as it appears to be less effective as monotherapy in East Africa (especially Sudan) than in India for a reason not yet understood (Hailu et al., 2010; Musa et al., 2010).

Drug combination is an established approach for the treatment of several infectious diseases (e.g. malaria, HIV-1 and tuberculosis) and is starting to be considered for the treatment of tropical diseases like visceral leishmaniasis. The appeal for the development of combination therapies is two-fold: to reduce the length of treatment in order to improve compliance (e.g. the miltefosine regimen involves a 4 weeks, twice daily drug intake), and to delay the emergence of resistance to protect the few molecules available. Nonetheless, there is still a risk that parasites could develop resistance to combination of drugs, being of particular importance those ones involving paromomycin as shown for L. donovani (Garcia-Hernandez et al., 2012). Several attempts to combine drugs have been made for the treatment of visceral leishmaniasis but as this falls outside the scope of this review the reader is referred to two recent reviews on the topic (Olliaro, 2010; van Griensven et al., 2010).

In this review, we will summarize the most recent highlights that next generation sequencing (NGS) brought to our understanding of drug resistance in the parasite Leishmania.

Next-generation sequencing technologies

NGS sequencing platforms have become widely available in recent years with associated costs now being orders of magnitude lower than classic Sanger sequencing. Different NGS platforms have reached the market during the past decade, including Illumina, Roche’s 454 pyrosequencing, Life Technologies’ Ion Torrent (PGM/Proton), Applied BioSystems’ SOLiD System, and Pacific Biosciences’ PacBio sequencing system. Each platform has specific advantages and disadvantages and it is important to carefully figure out which instrument is best suited to a given laboratory’s needs in terms of reads accuracy, reads length, throughput and cost. A thorough description of the different sequencing technologies is available in (Liu et al., 2012; Bahassi and Stambrook, 2014; Buermans and den Dunnen, 2014; van Dijk et al., 2014), for which the main highlights are summarized here. For a long time, the main advantage of 454 pyrosequencing was the long read lengths of 500–1000 bases produced in relatively fast run times. However, this technology suffers from low throughput, high cost and high error rates in homopolymer repeats. Owing to the ever increasing read lengths from competitor platforms, this technology has now become less cost efficient and will phase out in mid-2016. SOLiD sequencing provides high throughput and very high read accuracy, owing to the fact that each base is read twice during sequencing, but suffers from short read lengths and relatively long run times. Illumina has emerged as the leader in the field of NGS and is now offering the highest throughput and the lowest per-base cost. The read length of up to 300 bases also makes this technology compatible with most applications. This is the NGS platform that has been used for producing most of the Leishmania genome reported in the literature to date. The semi-conductor technology of Life Technologies’ PGM/Proton platforms benefits from fast run times because it does not rely on optical scanning of fluorescent nucleotides for sequencing. The PGM/Proton technology also generates reads whose lengths are compatible with most applications but suffers from high error rates in homopolymer repeats. For now, the sequence yield per run is less than that of Illumina’s HiSeq systems but the PGM/Proton’s platform is evolving at a very rapid pace.

The main advantage of the PacBio system compared to other platforms consists in the long reads that it generates (i.e. up to 20 kb and even longer). The PacBio system also enables to sequence individual template molecules, in contrast to the other platforms which are not sensitive enough to detect the extension of single nucleotide at the level of individual molecules and thus require a local clonal amplification of the initial template to increase the signal-to-noise ratio. PacBio sequencing suffers from high cost, high error rates and low throughput however and is thus suitable to a limited range of applications like increasing the quality of de novo genome assemblies. Indeed, owing to their small read length, typical NGS platforms are very sensitive to the presence of DNA repeats but combining high-quality short reads (either derived from Illumina or PGM/Proton platforms) with long (less accurate) PacBio reads now enables the reconstruction of long error-free reads that crosses many more DNA repeats and missing bases for producing very-high quality genome assemblies (Koren et al., 2012). This should prove especially useful for decreasing the number of gaps from genome assemblies for which DNA repeats abound as in the case of Leishmania (Ubeda et al., 2014) and for detecting chromosome rearrangements. In our experience, PacBio sequencing successfully reduced the number of L. infantum contigs by one order of magnitude, from 2000 contigs for paired-end Illumina sequencing alone to less than 175 contigs when including PacBio sequencing (N50 723,631 bp).

The road from raw sequence data to final results involves a number of analytical steps. A generalised analysis pipeline for NGS data includes processing the reads to remove adapter sequences and low-quality reads, aligning the reads to a known reference genome and/or reconstructing genome sequence by de novo assembly. Several free software packages are available to facilitate NGS data analysis, ranging from alignment tools like MAQ (Li et al., 2008a), BWA (Li and Durbin, 2009), Bowtie2 (Langmead and Salzberg, 2012), SOAP (Li et al., 2008b), SSAHA2 (Ning et al., 2001) and SMALT (www.sanger.ac.uk/resources/software/smalt/) to de novo assembler like Velvet (Zerbino and Birney, 2008), Ray (Boisvert et al., 2010), ALLPATHS-LG (Gnerre et al., 2011), SOAPdenovo2 (Luo et al., 2012), SGA (Simpson and Durbin, 2012) and ABySS (Simpson et al., 2009). Owing to the numerous tools available, selecting the software suitable for specific applications is not trivial and the reader is referred to studies that conducted comparison of the performance of major software for alignment (Fonseca et al., 2012) and de novo assembly (Salzberg et al., 2012). A comprehensive list of the tools available for the analysis of NGS data has been produced (Henry et al., 2014) and can be accessed through the OMIC tools portal (www.omictools.com). Finally, the Broad Institute’s Genome Analysis Toolkit provides a structured programming framework that offers a wide variety of tools for analyzing next-generation re-sequencing data, with a primary focus on the discovery of high quality variant sites (McKenna et al., 2010; DePristo et al., 2011).

The Leishmania genomes

Several high-quality draft or finished genomes are available for Leishmania species conferring visceral diseases (L. infantum JPCM5 and L. donovani BPK282/0cl4) (Peacock et al., 2007; Downing et al., 2011), cutaneous diseases (L. major Friedlin, L. mexicana U1103) (Ivens et al., 2005; Rogers et al., 2011) or muco-cutaneous diseases (L. braziliensis M2904) (Peacock et al., 2007). The genomes of L. tarentolae (Raymond et al., 2012) and of L. amazonensis (Real et al., 2013) have also been produced. The genomes of L. infantum, L. donovani, L. major and L. tarentolae each consist of 36 chromosomes (Ivens et al., 2005; Peacock et al., 2007; Downing et al., 2011; Raymond et al., 2012). On the other hand, the genome of L. braziliensis contains 35 chromosomes due to a fusion event between chromosomes 20 and 34 (Britto et al., 1998; Peacock et al., 2007), while the genome of L. mexicana consists of 34 chromosomes owing to two fusions events between chromosomes 8 and 29 and chromosomes 20 and 36 (Britto et al., 1998; Rogers et al., 2011). Draft genomes have also been produced for a number of other Leishmania species and these are available at TriTrypDB (www.trytripdb.org).

Globally, the genomes are characterized by a high degree of synteny between species despite an estimated 20–100 million years of separation between the L. Viannia spp. and the L. Leishmania spp. (Peacock et al., 2007). Comparative genomics also revealed a remarkably small number of species-specific genes, with less than twenty genes or paralogous groups specific to either L. infantum, L. major or L. mexicana (Britto et al., 1998; Rogers et al., 2011). The number of unique paralogous groups is slightly higher for L. braziliensis (67 genes) and L. tarentolae (95 genes), which is consistent with their outlier classification in the Viannia subgenus and the SauroLeishmania subgenus (which exclusively parasitize lizards), respectively. Moreover, L. tarentolae lacks genes associated with the intracellular life stage of pathogenic species (Raymond et al., 2012) and L. braziliensis holds features absent from Leishmania Leishmania ssp. such as functional RNA interference machinery as well as transposable elements (Peacock et al., 2007; Lye et al., 2010). This scarcity of species-specific genes suggests that disease tropism and variations in drug susceptibility between species are likely to involve differences in gene expression and gene dosage of common core genes, as previously shown for miltefosine for which the level of expression of the transporter responsible for miltefosine intake correlated with intrinsic variations in susceptibility between Leishmania species (Yardley et al., 2005; Sanchez-Canete et al., 2009).

Leishmania parasites have a peculiar mechanism of gene expression. Protein-coding genes are distributed on chromosomes as part of contiguous units on the same DNA strand and are co-transcribed as polycistronic transcription units by RNA-polymerase II before primary transcripts are matured into messenger RNAs (mRNA) by trans-splicing for the addition of a 39 nucleotides splice leader cap and poly-A tail (reviewed in (Haile and Papadopoulou, 2007)). Transcription initiation is bidirectional and takes place at strand-switch regions separating divergent transcriptional units (Martinez-Calvillo et al., 2003; Thomas et al., 2009). Conversely, transcription terminates at the level of strand-switch regions specified by convergent transcriptional units (Martinez-Calvillo et al., 2004). In the absence of promoters for regulating gene expression, the level of mRNA is controlled by RNA stability rather than by increased RNA polymerase activity (Muller et al., 2010a, b). To adapt to stressful environments (including drug pressure) in the absence of transcriptional control, Leishmania has evolved mechanisms to alter mRNA levels that include increased gene dosage through gene amplification (Ubeda et al., 2008; Leprohon et al., 2009b) or duplication (Mukherjee et al., 2011), gene deletion (Ouameur et al., 2008) and aneuploidy (Ubeda et al., 2008; Leprohon et al., 2009b). Conversely, single nucleotide polymorphisms in drug targets or key enzymes (Ritt et al., 2013) or transporters (Perez-Victoria et al., 2003; Coelho et al., 2012) can also favor adaptation without the need for altering gene expression.

Gene copy numbers

Copy number variations either through gene amplification or gene expansion enable generating substantial genome variability and Leishmania cells selected for resistance to cytotoxic compounds often amplify or delete a number of specific loci coding for either drug targets or drug transporters (reviewed in (Beverley, 1991; Ouellette and Papadopoulou, 1993)). The first example of gene amplification was the DHFR-TS gene in methotrexate resistant Leishmania (Coderre et al., 1983) but several locus have been found amplified after selection with a number of drugs (Ouellette and Papadopoulou, 1993). Amplified DNAs are generally extrachromosomal and found either as circular elements (Beverley et al., 1984; Garvey and Santi, 1986; Petrillo-Peixoto and Beverley, 1988; White et al., 1988) or as linear minichromosomes (Beverley and Coburn, 1990; Tripp et al., 1991; Papadopoulou et al., 1993; Navarro et al., 1994; Olmo et al., 1995; Grondin et al., 1998). While circles are formed by homologous recombination between direct repeated sequences (Ouellette et al., 1991; Grondin et al., 1993; Ubeda et al., 2008; Leprohon et al., 2009b), minichromosomes are generated by annealing of inverted repeats which lead to palindrome formation (Navarro et al., 1994; Olmo et al., 1995; Dubessay et al., 2001; Ubeda et al., 2008). Repeated sequences used for DNA amplification are generally non-coding and are highly conserved between different Leishmania species (Ubeda et al., 2008). Interestingly, a recent bioinformatics analysis revealed that such sequences are widespread in the Leishmania genomes and that most of the Leishmania genome is stochastically subjected to gene rearrangements at the level of these low-copy repeated sequences (Ubeda et al., 2014). Cells in the population are having a common core genome but differ in their extrachromosomal amplicon complement. This translates into an extensive heterogeneity at the population level and introduces phenotypic leaps that were shown to foster the selection of adaptive traits in response to drug pressure (Ubeda et al., 2014).

Gene amplification was first observed by comparison of digested genomic DNAs of sensitive and resistant parasites migrated on agarose gels and stained with ethidium bromide. The complexity of the DNA is low enough that it is possible to detect a 5–10 fold increase in gene copy (Ouellette et al., 1998). More recently, DNA microarrays have also allowed the detection of gene amplification given the strong correlation between RNA abundance and gene copy number (Guimond et al., 2003; Ubeda et al., 2008; Leprohon et al., 2009b; do Monte-Neto et al., 2011; Kumar et al., 2013). Now, assessing read depth coverage from whole genome sequencing data allows the precise mapping of amplified genomic loci and assists in the identification of repeated regions involved in homologous recombination. A recent analysis of the genome sequence of seventeen L. donovani field isolates from the Indian subcontinent identified gene copy number variations that could discriminate isolates according to their susceptibility to sodium stibogluconate (SSG) (Downing et al., 2011). Among these, a four-gene extrachromosomal circular amplicon coding for a mitogen-activated protein kinase derived from chromosome 36 was detected in significantly higher copies in SSG-resistant than in SSG-sensitive isolates. This difference in amplicon counts was stable as it was preserved despite the thirty in vitro passages between the collection of isolates and DNA extraction (Downing et al., 2011). While its role in antimony resistance remains to be functionally confirmed, the abundance of the four-gene amplicon could be used as a biomarker for the detection of antimony resistance in the event that its distribution proves to be widespread among SSG-resistant isolates. Amplicons derived from the H locus harboring the ABC gene MRPA were also observed in the L. donovani lines but in similar amounts between the SSG-resistant and SSG-sensitive isolates. Amplicons generated from the H locus were previously found not to be stable in the absence of drug pressure however (Leprohon et al., 2009b) and the lack of correlation with resistance might thus come from the culturing of isolates prior to DNA extraction (Downing et al., 2011). Interestingly, a locus previously reported to be amplified in L. tarentolae selected for resistance to SSG and whose abundance correlated with the level of SSG resistance (Haimeur and Ouellette, 1998) was present in more copies in the L. donovani SSG-resistant lines (Downing et al., 2011).

In addition to extrachromosomal gene amplification, altered gene dosage can also occur by gene duplication or gene deletion. Tandem gene arrays have been identified in every Leishmania species analysed. Interestingly, a third of these multicopy genes are species-specific and increased gene dosage arising from such tandem duplication may thus allow variations in transcript levels for some genes defining important species-specific features (Rogers et al., 2011). Most of the tandem arrays that are unique to each species are hypothetical proteins however and further investigation is warranted to examine these hypothetical genes for their role in the biology of Leishmania parasites (Rogers et al., 2011). Intrachromosomal segmental amplification through rearrangements at the level of direct repeats can also occur while attempting gene knockouts as shown for the subtelomeric essential gene gsh1 encoding gamma-glutamylcysteine synthetase, the rate limiting step in the biosynthesis of glutathione in Leishmania (Mukherjee et al., 2011). Chromosomal repeat expansion enabled the parasites to circumvent the sequential rounds of gene inactivation by becoming polyploid for the gsh1 locus (Mukherjee et al., 2011). Interestingly, intrachromosomal rearrangements leading to altered gene dosage were also observed in the context of antimony resistance in a series of L. major laboratory-derived mutants (Mukherjee et al., 2013). Read depth coverage from next generation sequencing data revealed that while chromosome 31 was at least hexasomic in every strain including wild-type cells, a subtelomeric region was deleted for at least four of the chromosome 31 alleles in all mutants (Mukherjee et al., 2013). The breaks mostly occurred at the level of inverted repeated sequences and differed between mutants. Most interestingly, careful analysis of sequencing reads revealed an even more heterogeneous profile at the population level where two of the three mutants were found to include subpopulations of parasites with different breakpoints for the terminal deletion of chromosome 31 (Mukherjee et al., 2013). The deleted region harbored the gene coding for the aquaglyceroporin AQP1 which is responsible for the uptake of trivalent antimony in Leishmania (Gourbal et al., 2004). Decreasing the number of AQP1 copies reduces trivalent antimony entry (Gourbal et al., 2004) and confers resistance (Gourbal et al., 2004; Marquis et al., 2005). In addition to AQP1 deletions, a subtelomeric amplification of part of chromosome 34 was also deduced from the sequencing reads in one of the L. major antimony-resistant mutants (Mukherjee et al., 2013). While an increase in copy number of a segment close to a telomere is often indicative of a linear amplicons (Ubeda et al., 2008), the chromosome profile on pulse field gel electrophoresis correlated instead with a complex intrachromosomal rearrangement (Mukherjee et al., 2013). Two genes on the amplified locus, encoding an ascorbate-dependent peroxidase and a putative glucose-6-phosphate dehydrogenase, were further linked with antimony resistance. Both genes are central to the redox defence in Leishmania and their increased dosage was shown to protect parasites against the reactive oxygen species that antimony is known to induce (Mehta and Shaha, 2006; Moreira et al., 2011; Mukherjee et al., 2013).

Chromosome copy numbers

Another layer of complexity for altering gene dosage can be found at the level of whole chromosome copy numbers. Additional chromosome copies are known to occur in various Leishmania species when attempting to knockout essential genes (Cruz et al., 1993; Mottram et al., 1996; Dumas et al., 1997; Martinez-Calvillo et al., 2005). In recent years, assessment of read depth coverage from next-generation sequencing data from reference Leishmania genomes has shown the presence of supernumerary chromosomes in every species analysed and confirmed that aneuploidy is much more widespread in Leishmania than initially anticipated, being common even in cells not submitted to any stress besides culturing (Rogers et al., 2011). The extent of aneuploidy varies between species or strains, and L. infantum/L. donovani appears to be the most aneuploid at the population level (Downing et al., 2011; Rogers et al., 2011). Also, an astonishing small number of chromosomes appear to be disomic in all species of the Leishmania subgenus (Downing et al., 2011; Rogers et al., 2011). Interestingly, the fact that many of these are enriched for gene arrays suggests that there may be a selection against whole-chromosome duplication for chromosomes harboring a higher proportion of putatively dose-sensitive genes in Leishmania species (Rogers et al., 2011). Chromosome 31, which was first shown to be supernumerary in L. major (Akopyants et al., 2009), is the only chromosome whose copy number is increased in every Leishmania species (Downing et al., 2011; Rogers et al., 2011; Mukherjee et al., 2013), even for the non-pathogenic L. tarentolae (Raymond et al., 2012). L. braziliensis is distinct from the species belonging to the Leishmania subgenus in being trisomic for most chromosomes (Rogers et al., 2011) but it remains to be seen whether this is maintained in other species from the Viannia subgenus.

Chromosome copy numbers calculated from read depth coverage can sometimes lead to a cumulative ploidy not matching with a clear-cut number of chromosomes but instead to intermediate somy values (e.g. median chromosome coverage between disomy and trisomy). Fluorescence in situ hybridization of single cells explained this phenomenon by revealing that distinct aneuploidy patterns also occur at the level of individual cells within a given population (Sterkers et al., 2011; Lachaud et al., 2014). Indeed, every chromosome analysed were show to be present in a least two ploidy states (one predominant and other minor states) in population of parasites and, similarly to the stochastic gene amplification described in section 3, this mosaic aneuploidy was proposed to generate intra-strain heterogeneity that should ultimately translate into increased adaptability at the population level (Sterkers et al., 2011).

In the context of drug resistance, aneuploidy has first been observed while performing comparative gene expression analyses by DNA microarrays on Leishmania parasites resistant to antimony (Leprohon et al., 2009b) or methotrexate (Ubeda et al., 2008), and later in L. major antimony-resistant mutants (Mukherjee et al., 2013) and in L. donovani selected for resistance to the HIV-1 protease inhibitor Nelfinavir (Kumar et al., 2013). Increased and decreased mRNA levels were observed for entire chromosomes and comparative genomic hybridization further confirmed this being related with the presence of supernumerary chromosomes and chromosome loss, respectively. In the case of L. infantum and L. major mutants resistant to antimony, this was further validated by computing chromosomal read depth coverage from whole genome sequencing data (Brotherton et al., 2013; Mukherjee et al., 2013). In all cases tested, withdrawal of drug pressure resulted in particular supernumerary chromosomes reverting to the somy observed prior to the induction of resistance, suggesting that aneuploidy can assist drug resistance (Ubeda et al., 2008; Leprohon et al., 2009b; Kumar et al., 2013). The levels of mRNA were also consistent with the altered chromosome stoichiometry but whether a specific gene or group of genes on aneuploid chromosomes is responsible for resistance has not been directly addressed, and the link between aneuploidy and drug resistance thus remains circumstantial. Quantitative proteomics by stable isotope labeling of amino acids in cell culture from the same L. infantum antimony-resistant mutant as in (Leprohon et al., 2009b) failed to correlate chromosome copy numbers with protein levels however but this is possible if the bulk of extra mRNA molecules is not translated (expect for those extra transcripts implicated in resistance) or if their protein products are degraded shortly after synthesis, as previously demonstrated in aneuploid yeast (Torres et al., 2007). Of note, resistance resulting from whole chromosome aneuploidy in fluconazole-resistant Candida albicans (Selmecki et al., 2006, 2010) and Crytptococcus neoformans (Sionov et al., 2010) was attributed to the increased expression of two specific genes located on supernumerary chromosome 5 (Selmecki et al., 2006, 2008). Similarly, adaptation of budding yeasts to the Hsp90 inhibitor radicicol induces duplication of chromosome 15 for which the synergistic effect conferred by the increased dosage of two genes encoding an Hsp90 co-chaperone and a drug pump on the chromosome was shown to be responsible for most of the resistance (Chen et al., 2012). Hence, although altered chromosome stoichiometry can lead to expression pattern changes for many genes, adaptive traits resulting from aneuploidy can be attributed to dosage fluctuations of specific genes on supernumerary chromosomes.

Analysis of read depth coverage from whole genome sequencing data is just beginning to be used for assessing aneuploidy in the context of drug resistance in natural Leishmania populations. Given the very high prevalence of antimony resistance in Nepal and India, the studies described so far mainly focused on antimony resistant isolates recovered from the Indian subcontinent. The analysis of seventeen L. donovani isolates either sensitive or resistant to SSG revealed an extensive aneuploidy among the strains (Downing et al., 2011). Again, few chromosomes were disomic in all samples (those correlated well with disomic chromosomes from (Rogers et al., 2011)) and only chromosome 31 was increased in every strain. Twenty-six chromosomes had varying ploidy across the different lines. No significant pattern of aneuploidy could be associated with the SSG resistant phenotype however and each isolate had a unique karyotype. Conversely to the laboratory-derived mutants described above, the influence of chromosomal copy numbers on gene transcripts levels was not so clear given the poor correlation between the depths of coverage for cDNA-derived and genomic DNA-derived sequencing reads (Downing et al., 2011). This analysis was only performed for one out of the seventeen L. donovani strains however and might thus have to be tested more widely.

Single nucleotide polymorphisms

Besides gene dosage fluctuations, single nucleotide mutations constitute another strategy to circumvent drug pressure. Point mutations can lead to resistance by interfering with drug-target interactions (Arrebola et al., 1994) or by altering the permeability of the cell for the drug (Vasudevan et al., 1998, 2001). While every antimony resistance mechanisms described so far implicated changes in gene copy numbers or transcripts levels, whole genome sequencing is beginning to reveal a role for point mutations in resistance. The genome sequence of a L. infantum mutant selected in vitro for resistance to antimony indeed revealed a guanine to adenine transition at position 1885 of the protein kinase LinJ.33.1810, giving rise to an E629K substitution in the protein (Brotherton et al., 2013). Overexpression of the wild-type version of the gene conferred a modest but reproducible and significant sensitization of the mutant, hence implicating the mutation in resistance (Brotherton et al., 2013). Interestingly, a mitogen-activated protein kinase 1 was recently found to be downregulated in L. donovani antimony-resistant field isolates and to correlate with resistance (Ashutosh et al., 2012). Antimony is a known protein phosphatase inhibitor (Pathak and Yi, 2001) and it is possible that mutations (or a decreased expression) in protein kinases may compensate for the inhibition of phosphatase activities. In L. guyanensis, the selection of resistance to antimony was accompanied by a reduced accumulation of the drug mainly due to genomic alterations at the level of the subtelomeric region of chromosome 31 harboring the gene coding for AQP1. For one mutant, the transport defect was correlated with the presence of a single nucleotide polymorphism at position 398 of AQP1 identified by whole genome sequencing, leading to the replacement of a glycine residue by an aspartic acid at position 133 (G133D) of the protein. The role of this missense mutation in resistance was functionally confirmed by overexpression experiments, whereby a wild-type AQP1 but not its altered AQP1G133D version was able to re-sensitize AQP1-null L. guyanensis parasites to antimony (Monte-Neto, et al., submitted).

Whole genome sequencing also revealed several mechanisms of resistance to 5-fluorouracil including point mutations while seeking for further insights into pyrimidine metabolism in Leishmania (Ritt et al., 2013). 5-fluorouracil is a pyrimidine analog specifically inhibiting thymidilate synthase that displays anti-leishmanial activity (Katakura et al., 2004). Besides amplification of the dhfr-ts target locus, point mutations in an uridine phosphorybosyl transferase, a thymidine kinase and an uridine phosphorylase, all linked to pyrimidine metabolism in kinetoplastids (Hammond and Gutteridge, 1982; Ali et al., 2013), were observed in three independent laboratory-derived L. infantum 5-fluorouracil resistant mutants (Ritt et al., 2013). All mutations were further implicated in resistance to 5-fluorouracil by functional experiments. Interestingly, these experiments also highlighted a connection between pyrimidine metabolism and folate/antifolate metabolism, where a mutation in thymidine kinase implicated in resistance to 5-fluorouracil also caused hypersensitivity to methotrexate (Ritt et al., 2013). This could prove useful for the development of specific antimetabolites against Leishmania.

The uptake of miltefosine and other alkyl-glycerophospholipids in Leishmania is mediated by a P-type ATPase named miltefosine transporter (MT) which is responsible for the translocation of phospholipids from the exoplasmic to the cytoplasmic leaflet of the plasma membrane of the parasite (Perez-Victoria et al., 2003). Miltefosine uptake further requires a protein named Ros3, the beta-subunit of the MT (Perez-Victoria et al., 2006). Targeted analyses of the MT gene in miltefosine-resistant mutants previously revealed that different inactivating mutations are able to drastically increase resistance to miltefosine (Perez-Victoria et al., 2003). Whole genome sequencing data from two laboratory-derived L. major mutants also identified several distinct MT mutations located at conserved residues and conferring resistance (Coelho et al., 2012). Interestingly, the analysis of the MT status in intermediate steps mutants indicated that mutations in MT are polyclonal, with up to five different MT genotypes identified within a given population. The number of mutated alleles also correlated with the level of resistance to miltefosine, early selection steps being more heterozygous for MT point mutations than highly resistant strains which were primarily homozygous for the mutations (Coelho et al., 2012). Interestingly, a similar stepwise increase in the number of mutated alleles was also observed for mutations in the gene LinJ30.1250 coding for a pyridoxal kinase detected in the two independent L. major mutants (Coelho et al., 2012). These mutations were also polyclonal, with at least two subpopulations differing in the number of mutated alleles being found in a given population, and were implicated in low-level resistance to miltefosine (Coelho et al., 2012). Pyridoxal kinase is involved in the scavenging of vitamin B6 precursors but the mechanism by which it confers resistance to miltefosine remains to be clarified.

Genome-wide SNPs typing by next-generation sequencing represents a powerful approach for differentiating parasite strains among natural populations and proved better than the more traditional multi-locus microsatellite typing in informing on the population history of closely related Leishmania strains (Downing et al., 2012). Principal component analysis of genome-wide nucleotide variations between several L. donovani field isolates from the Indian subcontinent displaying varying levels of antimony susceptibility supported the notion of multiple events of emergence of SSG resistance that was previously reported (Laurent et al., 2007; Downing et al., 2011, 2012). Evidence of adaptive evolution were also noted for some genes that have a role in drug resistance (Downing et al., 2011), with a high level of diversity detected in ATP-binding cassette transporters (Castanys-Munoz et al., 2008; Leprohon et al., 2009a), histone genes (Singh et al., 2007, 2010) and in the nucleoside transporter 1 (Vasudevan et al., 2001). While none of the SNPs identified was universal to the SSG-resistant isolates (Downing et al., 2011), a panel of five homozygous SNPs specifically detected in SSG-resistant isolates proved useful for classifying the L. donovani lines into eight multilocus sequencing typing haplotypes (Vanaerschot et al., 2012). Interestingly, three haplotypes were significantly more associated with non-responsive or relapsing patients compared to those that enjoyed definite cure (Vanaerschot et al., 2012). These haplotypes gathered only few isolates however and additional confirmation with a larger number of strains isolated from patients with well-defined SSG treatment outcomes are required to truly assess the potential of these haplotypes in predicting SSG treatment failure.

Natural extensions of NGS for studying the biology of Leishmania

In addition to whole genome sequencing, NGS allows for capturing the most significant changes in gene expression by means of RNA sequencing (RNA-seq), an approach that has several advantages over other transcriptional profiling approaches like DNA microarrays or serial analysis of gene expression (Nowrousian, 2010). Although RNA-seq studies addressing drug resistance in Leishmania still have to be reported, the approach has already proved valuable in providing new knowledge about the heterogeneity of spliced leader and polyadenylation sites for Leishmania mRNA, in addition to allow the identification of several new transcripts (Rastrojo et al., 2013). A dual host-parasite RNA-seq experiment also recently highlighted a potential transcriptional signature in primary L. braziliensis cutaneous lesions that may predict the long-term development of the mucosal form of the disease (Maretti-Mira et al., 2012). RNA-seq also helped deciphering the role of iron, iron superoxide dismutase and reactive oxygen species in the control of amastigote differentiation in L. amazonensis (Mittra et al., 2013). Owing to the unusual mechanisms of gene regulation in Leishmania, one should be aware however that RNA-seq experiments in Leishmania should greatly benefit from companion proteomics profiling, as exemplified by the analysis of the proteomic (by LC-MS/MS) and transcriptomic (by RNA-seq) responses of L. donovani to purine starvation that revealed that only a minority of changes detected at the protein level tracked with changes at the mRNA level (Martin et al., 2014).

Chromatin immunoprecipitation followed by hybridization on DNA microarrays (ChIP-Chip) or NGS (ChIP-Seq) allow studying the occupancy of DNA-binding proteins on the chromatin and the interplay between chromatin structure and transcription. Eukaryotic DNA is wrapped around nucleosomes composed of H2A, H2B, H3 and H4 histone proteins, with histone H1 contributing to larger order structures. By substituting histone variants or through histone modifications, eukaryotes can regulate the availability of chromatin to proteins involved in transcription. For example, acetylated histone H3 is usually found at the 5′ end of transcription start sites and is associated with increased transcription rates (Rando, 2007). While the mechanism of transcription initiation in Leishmania is poorly understood, ChIP-Chip and ChIP-Seq revealed that specific epigenetic marks are also present at putative transcriptional start and termination sites in this parasite (Thomas et al., 2009; van Luenen et al., 2012) and the related trypanosomatid Trypanosoma brucei (Siegel et al., 2009). ChIP-Seq experiments have also been conducted in our lab for mapping the precise location of the trimethylation state of lysine 4 of histone H3 (3meH3K4) in the context of drug resistance. This epigenetic mark is associated with a loosening of the chromatin and active transcription and was found to be enriched in the extrachromosomal element of drug-resistant strains (Gazanion, unpublished results).

Finally, NGS provides a unique opportunity for developing innovative strategies to discover and validate different aspects of Leishmania biology, including the identification of the mode of action and resistance mechanisms of antimicrobials. We have recently developed a multi-copy suppressor screening strategy that combines functional cloning and NGS for the detailed characterization of drugs’ mode of action and resistance determinants in Leishmania. The technique, termed Cos-seq, consists in generating genomic libraries in multi-copy cosmids (i.e. Cl-Hygro) that are subsequently transformed into drug sensitive parasites. The fitness conferred by each cosmid can then be analyzed by tracking the fate of individual transformants grown as pools in liquid media in the presence of incremental increases in drug pressure. By comparing the number of sequence reads per gene between drug increments, one can observe a gradual selection of genomic regions that bestow a selective advantage and a rapid clearance of those conferring increased sensitivity. The approach was validated by identifying known drug targets and resistance mechanisms for methotrexate and this technique was applied to a number of antileishmanial drugs and several novel resistance determinants have been identified.

Conclusions

Whole genome sequencing has become affordable for large scale studies in recent years and proved useful at detecting both known and new mechanisms of drug resistance in Leishmania. These also highlighted the complexity of drug resistance which is often polyclonal and involves several different co-existing genotypes in a given population. This was observed for episome formation, where stochastic amplification of random loci in unstressed population can favor adaptation; for gene deletion or duplication, with distinct breakpoints occurring between parasites of a given population; for aneuploidy, where the impressive variation in chromosome copy numbers between species and strains renders the link with resistance more difficult to establish; and for SNPs, where different mutations leading to resistance occur in a given population of parasite. This heterogeneity in both genotype and phenotype in resistant populations may explain the difficulties when studying field isolates. Indeed, the polyclonal nature of resistance can lead to mutations going unnoticed when analyzing parasite populations because these will constitute rare events. In the situation where no clones are dominant in the population, the analysis of several clones and looking for similar events (e.g. different mutations targeting the same gene between mutants) might thus be more revealing than population analysis for the identification of drug resistance mechanisms and for characterizing the mechanisms implicated in their genesis (e.g. to precisely map the breakpoints implicated in the deletion or the amplification/duplication of a given locus). Finally, additional studies will be required to understand more deeply the clonal variety within resistant population and how it contributes to either resistance or fitness, as recently highlighted for populations of antibiotic-resistant bacteria (Lee et al., 2010; Vega et al., 2012; El-Halfawy and Valvano, 2013).

Conflict of interest

The authors declared that there is no conflict of interest.