Genomic and transcriptomic alterations in Leishmania donovani lines experimentally resistant to antileishmanial drugs.

Leishmaniasis is a serious medical issue in many countries around the World, but it remains largely neglected in terms of research investment for developing new control and treatment measures. No vaccines exist for human use, and the chemotherapeutic agents currently used are scanty. Furthermore, for some drugs, resistance and treatment failure are increasing to alarming levels. The aim of this work was to identify genomic and trancriptomic alterations associated with experimental resistance against the common drugs used against VL: trivalent antimony (SbIII, S line), amphotericin B (AmB, A line), miltefosine (MIL, M line) and paromomycin (PMM, P line). A total of 1006 differentially expressed transcripts were identified in the S line, 379 in the A line, 146 in the M line, and 129 in the P line. Also, changes in ploidy of chromosomes and amplification/deletion of particular regions were observed in the resistant lines regarding the parental one. A series of genes were identified as possible drivers of the resistance phenotype and were validated in both promastigotes and amastigotes from Leishmania donovani, Leishmania infantum and Leishmania major species. Remarkably, a deletion of the gene LinJ.36.2510 (coding for 24-sterol methyltransferase, SMT) was found to be associated with AmB-resistance in the A line. In the P line, a dramatic overexpression of the transcripts LinJ.27.T1940 and LinJ.27.T1950 that results from a massive amplification of the collinear genes was suggested as one of the mechanisms of PMM resistance. This conclusion was reinforced after transfection experiments in which significant PMM-resistance was generated in WT parasites over-expressing either gene LinJ.27.1940 (coding for a D-lactate dehydrogenase-like protein, D-LDH) or gene LinJ.27.1950 (coding for an aminotransferase of branched-chain amino acids, BCAT). This work allowed to identify new drivers, like SMT, the deletion of which being associated with resistance to AmB, and the tandem D-LDH-BCAT, the amplification of which being related to PMM resistance.

Introduction

Leishmaniasis is a complex of diseases caused by protists of the genus Leishmania. Many species of this genus are human pathogens and are responsible for different clinical forms of the disease, among others cutaneous (CL), mucosal and visceral leishmaniasis (VL). These infectious diseases are endemic worldwide, in tropical and sub-tropical regions, including also the Mediterranean basin (Dujardin et al., 2008). Recent estimates indicate that there are approximately 0.4 million new VL cases and 1.2 million CL cases every year (Alvar et al., 2012). Moreover, the prevalence of VL co-infection with human immunodeficiency virus (HIV) is increasing, and HIV infection can increase the risk of patients to develop VL, reducing their ability to respond to chemotherapy, and increasing the probability of disease relapse (van Griensven et al., 2014). Unfortunately, there are currently no effective vaccines against this parasite (Singh and Sundar, 2012). Thus, control of leishmaniasis relies primarily on chemotherapy and vector control (Singh et al., 2012). However, available drugs are limited in number, and toxicity together with emerging resistance might further jeopardize their use (Ouellette et al., 2004). Understanding drug resistance is important to protect the few existing compounds, but also for developing new ones (Hefnawy et al., 2017).

Drug resistance phenotype can be the result of different strategies followed by Leishmania (reviewed in (Ouellette et al., 2004)), including: decrease in drug uptake, efflux of drugs, inactivation of drugs (once inside the parasite), loss of drug activation pathways, alteration of drug targets, and finally the parasite may become more proficient in repairing drug damage. In addition, molecular adaptations related to drug resistance may provide some advantages to Leishmania donovani parasites that may be counterintuitive, such as a higher capacity to survive in stress conditions (Garcia-Hernandez et al., 2015), and an increased fitness (Vanaerschot et al., 2014).

Antimonials (Glucantime and Pentostam, SbV) were the first drugs used for treatment of leishmaniasis, almost eighty years ago and they continue being a standard first-line treatment in many places of the world. However, the appearance of resistant parasites has been described frequently in many endemic regions (Croft et al., 2006; Singh et al., 2012; Hefnawy et al., 2017). Several membrane transporters of the ABC superfamily have been involved in SbV resistance. PGPA/MRPA (ABCC3), was the first ABC transporter described to be responsible for clinical resistance to antimonials in L. donovani (Mittal et al., 2007; Mukherjee et al., 2007; Mukhopadhyay et al., 2011) and experimental resistance to trivalent antimony (SbIII) or the related metal AsIII in L. tarentolae (Legare et al., 2001). The second ABC protein involved in antimony resistance is PRP1 (ABCC7) (Coelho et al., 2003; Leprohon et al., 2009). Recently, two other ABC transporters, ABCI4 and ABCG2 (Manzano et al., 2013; Perea et al., 2016), have been involved in antimony resistance in L. major. These ABC transporters confer resistance to antimony due to a significant decrease in drug accumulation as a consequence of an efflux of metals, probably conjugated with thiols. Other mechanisms involved in SbV resistance in L. donovani include decreased drug uptake through inactivation of AQP1 (aquaglyceroporin) transporter, increased tolerance to oxidative stress through upregulation of antioxidant pathways (increase in thiol levels) and modifications in membrane fluidity (Sundar, 2001; Decuypere et al., 2005; Croft et al., 2006; Mittal et al., 2007; Imamura et al., 2016).

In recent years, a few other drugs have been used for treatment of leishmaniasis, including amphotericin B liposome formulation (AmBisome, AmB), oral miltefosine (MIL) and a parenteral formulation of aminosidine (paromomycin, PMM) (reviewed in (Croft et al., 2006; Monge-Maillo and Lopez-Velez, 2013; Sundar and Singh, 2016)). However, resistance was also reported for these drugs: (i) experimental resistance in L. donovani to AmB (Mbongo et al., 1998), MIL (Perez-Victoria et al., 2003b; Seifert et al., 2003) and PMM (Maarouf et al., 1998; Jhingran et al., 2009), and (ii) clinical resistance in L. donovani to AmB (Purkait et al., 2012, 2015) and in L. infantum and L. donovani to MIL (Mondelaers et al., 2016; Srivastava et al., 2017). Molecular mechanisms of AmB resistance of one clinical isolate of L. donovani have been described including an altered membrane composition (decreased ergosterol content) due to a modified expression of the two transcripts coding for the enzyme S-adenosyl-L-methionine: C24-sterol methyltranferase (SMT); the parasites show a decreased AmB affinity and uptake (Purkait et al., 2012). Additionally, the involvement of SMT in AmB susceptibility has been validated in yeast (Konecna et al., 2016). A recent publication has demonstrated that a mutation in the sterol 14-α-demethylase, another enzyme of the sterol biosynthetic pathway, leads to AmB resistance in Leishmania mexicana (Mwenechanya et al., 2017), supporting the relevance of sterol biosynthetic pathway in AmB susceptibility. Also, the expression of ABCB (ABCB4, MDR1) transporter was found to be higher in AmB resistant parasites, suggesting a high efflux of AmB (Purkait et al., 2012). In addition, upregulation of thiol metabolism and reduced intracellular thiol levels have been described in clinical isolates of L. donovani resistant to AmB (Purkait et al., 2012; Suman et al., 2016). MIL is the latest and unique oral antileishmanial drug; however, after a decade of use, recent reports indicate a significant decrease in its efficacy against VL (Rijal et al., 2013). Experimental resistance to MIL was easily obtained (Seifert et al., 2003; Perez-Victoria et al., 2003a). L. donovani resistant parasites showed changes in the length and levels of unsaturation of fatty acids, reduction in ergosterol levels (Rakotomanga et al., 2005) and a significant reduction in drug internalization due to mutations in the MIL transporter genes LdMT and/or LdRos3 (Perez-Victoria et al., 2003b; Sanchez-Cañete et al., 2009). Only two clinical isolates with demonstrated MIL resistance have been described, corresponding to L. infantum lines isolated from HIV-coinfected patients in France that show a MIL-resistant phenotype due to a frameshift mutation in the MIL transporter gene LiRos3 and/or SNPs in the LiMT gene (Mondelaers et al., 2016). PMM, a natural aminoglycoside antibiotic, has a high promising efficiency against VL and CL. Although rapid development of experimental resistance to PMM in L. donovani and L. infantum has been described (Maarouf et al., 1998; Jhingran et al., 2009; Hendrickx et al., 2015), no clinical resistant isolates to PMM have been reported, probably by its low routine use in the field. An increased membrane fluidity associated with a decreased intracellular PMM accumulation have been described in experimental PMM resistant L. donovani lines (Jhingran et al., 2009; Bhandari et al., 2014). Additionally, an increased expression of ABC transporters (MDR1 and MRPA) and protein phosphatase 2A was observed, together with an increased drug efflux and a greater tolerance of L. donovani parasites to host defense mechanisms including nitrosative, oxidative and complement-mediated stresses (Bhandari et al., 2014).

Next-generation sequencing (NGS) technologies are being used for expediting the untargeted identification of drug resistance markers in many organisms, and also in Leishmania (reviewed in (Leprohon et al., 2015)). In the absence of transcriptional regulation, a common strategy of the parasite is to modulate gene dosage (Requena, 2011; Berg et al., 2013), either locally by gene copy number variation (Inga et al., 1998) and/or by generating episomal amplicons (Beverley, 1991) or at whole chromosome level through aneuploidy (Sterkers et al., 2012; Imamura et al., 2016). Recently, a combination of genome-wide cosmid-based functional screening with NGS, termed as Cosmid Sequencing (Cos-Seq) has been successfully used for identification of drug target and resistance mechanisms in Leishmania (Gazanion et al., 2016). However, this method can only identify gain-of-function and individual loci involved in drug resistance. Drug resistance selection by progressive increase of exposure to the compound followed by NGS of both whole DNA and mRNA remains the most adequate method for identifying multi-genic adaptations in a broad context, including loss-of-function. In present study, this method was applied on a same L. donovani line allowing to study experimental resistance to either SbIII, AmB, MIL or PMM, within a same genetic background. Chromosomal reorganizations and substantial changes in gene expression were found to be associated with the different drug resistances. We further validated the contribution of selected drug resistance markers, allowing us to identify drivers of resistance to some drugs currently in use to treat leishmaniasis.

Materials and methods

Chemicals

Trivalent antimony (SbIII), paromomycin (PMM), amphotericin B (AmB), Triton X-100, paraformaldehyde, MTT [3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide], rhodamine 123 (Rh123), phorbol 12-myristate 13-acetate (PMA), Trizol (TRI) Reagent, were obtained from Sigma-Aldrich (St. Louis, MO). Miltefosine (MIL) was purchased from Zentaris GmbH (Frankfurt am Main, Germany). Hygromycin B (HyB), CellTracker, dichlorodihydrofluorescein diacetate (H2DCFDA), MitoSOX red, and 4’, 6-diamidino-2-phenylindole dihydrochloride (DAPI) were purchased from Invitrogen (Carlsbad, CA). Geneticin G418 Sulphate, L-glutamine, and penicillin/streptomycin were obtained from Gibco. Nourseothricin was obtained from Werner Biogents. All chemicals were of the highest quality available.

Leishmania culture conditions

Promastigotes of Leishmania donovani (MHOM/ET/67/HU3 and MHOM/IND/80/Dd8), Leishmania major (MHOM/JL/80/Friedlin), Leishmania infantum JPCM5 (MCAN/ES/98/LLM-877), and derivative lines used in this study were grown at 28 °C in RPMI 1640-modified medium (Invitrogen) supplemented with 20% heat-inactivated fetal bovine serum (HIFBS, Invitrogen) and stored in cryobanks. Leishmania lines overexpressing D-lactate dehydrogenase-like (D-LDH), branched-chain amino acid aminotransferase (BCAT), protein associated with differentiation (PAD), and sterol 24-c-methyltransferase (SMT) were obtained in this study and grown in the same medium in the presence of 100 μg/ml of HyB. In a previous work (Garcia-Hernandez et al., 2012), the lines A, M, P, and S (that are resistant to AmB, MIL, PMM, and SbIII, respectively) were selected in vitro from L. donovani HU3 promastigotes by a stepwise adaptation process till 80 μM SbIII, 0.1 μM AmB, 8 μM MIL and 20 μM PMM. Wild-type (WT) and resistant lines were grown during 30 passages under drug pressure before proceeding with extraction of both DNA and RNA for NGS sequencing. However, for the rest of purposes, the lines were cultured without drugs just before being subjected to the different assays.

Isolation of DNA and RNA

Around 4 × 108 promastigotes in the late logarithmic phase were harvested by centrifugation, and the pellet suspended in 1 ml of Trizol (TRI Reagent, Sigma-Aldrich, product No. T9424). Manufacturer's instructions were followed. Samples were kept at −70 °C for a week before proceeding with the phase separation. After thawing, 0.2 ml of chloroform was added, and the mixtures were shaken vigorously for 15 s. After centrifugation, three phases were observed: a red organic phase (containing protein), an interphase (containing DNA), and a colorless upper aqueous phase (containing RNA). Both interphase and aqueous phases were processed separately for isolation of DNA and RNA, respectively. RNA samples were suspended in DEPC-treated water, and their concentrations were determined using the Nanodrop ND-1000 (Thermo Scientific); all samples showed A260/A280 ratios higher than 2.0. In addition, RNA integrity was checked in a bioanalyzer (Agilent 2100). The DNA samples were also quantified by absorbance at 260 nm using the Nanodrop, and the integrity by electrophoresis in an 0.8% agarose gel.

Next-generation sequencing

Library preparations (mRNA and DNA) and sequencing were performed by Centro de Análisis Genómico (CNAG, Spain). Sequencing was done using Illumina HiSeq 2000 technology. From RNA libraries, paired reads of 75 nucleotides were obtained, whereas 2 × 100 reads were obtained from DNA libraries. Sequence quality metrics were assessed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

Mapping of RNA-Seq and DNA-Seq reads to reference genomes

Reads were mapped to the reference genomes (v.9; TriTrypDB.org) of L. infantum (JPCM5 strain (Peacock et al., 2007)) and L. donovani (MHOM/NP/03/BPK282/0; clone BPK282/0cl4 (Downing et al., 2011)) genomes using Bowtie2 with default parameters (Langmead and Salzberg, 2012). This aligner was chosen because the end-to-end mode has a better performance when working with a genome containing a high level of repeated regions such as tandemly repeated genes. Other aligners, by default, trim the reads to try to map as much reads as possible, and consequently produce an over-representation of repeated regions. As percentages of mapped reads were higher against the L. infantum genome than against the L. donovani one, the former was used as reference for subsequent analyses.

Assignment of RNA-Seq to specific transcripts and determination of differentially expressed transcripts

L. infantum (JPCM5 strain) transcriptome was delineated using RNA-Seq data derived from poly-A enriched RNA of log-phase promastigotes following the procedure described elsewhere (Rastrojo et al., 2013). Afterwards, the HTSeq program (v0.6.1 (Anders et al., 2015)) was used for counting the reads overlapping with the transcripts (mode: intersection-nonempty). These sets of values (3 replicates for each sample) were used for differential expression analysis by the DESeq2 program (v1.8.1 (Love et al., 2014)).

Chromosome ploidy analysis and determination of chromosomal alterations

The Illumina DNA-Seq reads were aligned to the L. infantum reference genome by Bowtie2 (parameters: --np 0 --n-ceil L,0,0.02 --rdg 0,6 --rfg 0,6 --mp 6,2 --score-min L,0,-0.06 –no-unal). SAMtools view (Li et al., 2009) (parameters: -Sh -f 2 -F 256) was used to discard multi-hit reads, and an in-house script was used to discard orphan reads. Again, SAMtools view was used to obtain the alignments in standard BAM format and to count the number of reads (parameters: -F 0x4). Afterwards, the data in BAM format were analyzed by genomeCoverageBed (BEDTools (Quinlan and Hall, 2010)) using the parameter –d to obtain the read depth per position along each chromosome and the genome. This procedure was repeated for each DNA-Seq sample (WT, S, A, M and P lines). A normalization of coverage values was done taking into account the total number of reads for each sample. In the P line, due to the over-amplification of the genomic region 857320–865727 of chromosome 27, this region was not considered for the calculation of the mean coverage of this chromosome. The somy of each chromosome was calculated with the following formula 2*di/dm (Mondelaers et al., 2016), where di is the median read depth of each chromosome and dm is the median depth of the genome.

For analyzing the read coverage in particular genomic regions, genomeCoverageBed (parameter –d) was again used to obtain the coverage by position in the selected genomic regions. Afterwards, an in-house Perl script was used to smooth read-depth values using overlapping windows of 5-Kb in length. For comparison between different samples, a further normalization considering the total number of reads was applied. Finally, graphs were generated from the coverage values by the tool plotly included in the R package (https://CRAN.R-project.org). Figure drawing was done using the software Inkscape (https://www.inkscape.org).

Gene cloning and transfection experiments

DNA fragments coding for D-LDH (LinJ.27.1940), BCAT (LinJ.27.1950), PAD (LinJ.13.1430), and SMT (LinJ.36.2510) were PCR amplified from genomic DNA of the L. donovani HU3 strain using the following primer pairs:

LDHBamHIF 5′-TATGGATCCA TGACGTACCC ACTGAAACA-3′ and LDHXbaIR 5′-ATATTCTAGA TTACAGGTGT GCCGTCGGGT-3′ for D-LDH; BCABglIIF; 5′-ATATAGATCT ATGCTGCTGA GCCGACGCTG-3′ and BCAXbaIR 5′-ATATTCTAGA AGCCTCAACC TTCACTGACC-3′ for BCAT; PADBglIIF 5′-ATATAGATCT ATGTTTGACA TTCGCGAGAA-3′ and PADXbaIR 5′-ATATTCTAGA CTACAGTGTC GTCTCCGGAG-3′ for PAD and, SMTBglIIF 5′-ATATAGATCT ATGTCCGCCG GTGGCCGTGA-3′ and SMTXbaIR 5′- ATATTCTAGA CTAAGCCTGC TTGGACGGCT-3′ for SMT. Restriction sites (in bold) were added for subsequent cloning. The PCR products were subcloned and sequenced into the pLEXSY-hyg2 expression vector (Jenabioscience, Jena, Germany) and pIR1SAT expression vector (from Prof S.M. Beverley, Washington University, USA), to generate the pLEXSY-DLDH, pLEXSY-BCAT, pLEXSY-PAD, pLEXSY-SMT, pIR1SAT-DLDH, and pIR1SAT-BCAT constructs, respectively. These plasmids were subsequently used for transfection of Leishmania promastigotes by electroporation in 2-mm gap cuvettes at 450 V and 500 μF (BTX Electro Cell Manipulator 600). Transfectants were selected in medium containing Hyg (100 μg/ml) for pLEXSY-hyg2 or nourseothricin (100 μg/ml) for pIR1SAT. The A line transfected with pLEXSY-hyg2 empty and WT line transfected with pLEXSY-hyg2 empty and co-transfected with both pLEXSY-hyg2 and pIR1SAT empty (controls) were named as A + Lexsy, WT + Lexsy and WT + Lexsy + pIRSAT, respectively, while WT line transfected with pLEXSY-DLDH, pLEXSY-BCAT and pLEXSY-PAD, and co-transfected with pLEXSY-DLDH + pIR1SAT-BCAT, and pLEXSY-BCAT + pIR1SAT-DLDH were named as WT + D-LDH, WT + BCAT, WT + PAD, WT + DLDH + BCAT, and WT + BCAT + DLDH, respectively. On the other hand, the A line transfected with pLEXSY-SMT was named as A + SMT.

For PCR amplification of the locus expanding genes LinJ.36.2510 and LinJ.36.2520, the following primers were used: gSMT-Fw, 5′- AAAGTTGGTG GCAACCTCC-3’; and gSMT-Rv, 5′- CCCAGGCAAG ACCACCTTAC-3’. As template, 25 ng of genomic DNA from either WT or A lines were employed, and the cycling conditions were: 98 °C for 30 s, followed by 30 cycles of 98 °C for 10 s, 65 °C for 20s and 72 °C for 3 min, and a final extension step of 10 min at 72 °C. Phusion Hot Start II High Fidelity DNA Polymerase (Thermo Scientific) was used. The 1.6 kb amplicon obtained in the A line was cloned into pNZY28 vector, using the NZY-blunt PCR cloning kit (Nzytech) following manufacturer's instructions. The sequence of the resulting clone was determined in the Sequencing Unit of Parque Científico de Madrid.

RT-qPCR gene expression

Total RNA was extracted from the indicated L. donovani lines during the early log phase of growth using the RNeasy Plus Mini Kit (Qiagen). Total RNA (400 ng) was retrotranscribed to cDNA with the qScript cDNA Synthesis Kit (Quanta Biosciences, Inc) following the manufacturer's instructions. The lack of genomic DNA contamination was checked by PCR amplification of RNA samples without prior cDNA synthesis. The specific primer pairs, designed by using the Primer3 software (Untergasser et al., 2007) and used to amplify cDNA, were: 5′-GTGGCAGTGG TGAAAGTGTG-3′ and 5′-TCGAAGTTCA TGCTGGTGAG-3′ for D-LDH, 5′-AGCGCCCGTA CTACATGTCT-3′ and 5′-CGATACGCAA TCATGCTCAC-3′ for BCAT, 5′CCACACAACA ACTCCAATGA TG-3′ and 5′-CACAGATGGA CGGAGAGGAG-3′ for PAD, 5′-CATGAGCTTA GCCGACAACA-3′ and 5′-GCACCACTCG TACAGGACAA AG-3′ for SMT and 5′-TTCTCGGCTT CACCAACGAC-3′ and 5′-CCATGTAGCG CACCAGGTC-3′ for GAPDH. Standard curves for each primer pair were generated with 10-fold serial dilutions of L. donovani genomic DNA to determine primer efficiency. Quantitative PCR was performed in a CFX96 cycler (BioRad); each 10 μl reaction were set up containing 5 μl of PerfeCta SYBR Green SuperMix (Quanta Biosciences), 500 nM of each primer and 2 μl of a 1:5 dilution of the synthetized cDNA. All reactions were performed in duplicate and the specificity of the amplification was verified by melting curve analysis. Gene expression data were normalized to the expression of the reference gene GAPDH, and relative to the control sample, using the CFX Manager software with ΔΔCt method (Livak and Schmittgen, 2001; Pfaffl, 2001).

Drug susceptibility in promastigotes and intracellular amastigotes of Leishmania lines

The susceptibility of Leishmania promastigotes to SbIII, AmB, PMM, geneticin, and nourseothricin was determined using the MTT colorimetric assay, after incubation for 72 h at 28 °C in the presence of increasing concentrations of the drug, as described previously (Kennedy et al., 2001). In this way, it was determined the 50% effective concentration (EC50), which is the concentration of drug that decreases the rate of cell growth by 50%, whereas the resistance index (RI) was calculated by dividing the EC50 for each resistant line or lines overexpressing D-LDH, BCAT, PAD or SMT, by that for the corresponding Leishmania control line (WT for A, P and S lines, and WT + Lexsy or A + Lexsy for the transfected lines).

To determine the susceptibility of intracellular Leishmania amastigotes to PMM, macrophage-differentiated-THP-1 cells were used. THP-1 cells were plated at a density of 3 × 105 macrophages/well in RPMI 1640 medium supplemented with 10% HIFBS, 2 mM glutamate, penicillin (100 U/ml), and streptomycin (100 mg/ml) in 24-well tissue culture chamber slides and, differentiated to macrophages with 20 ng/ml of PMA treatment for 48 h followed by 24 h of culture in fresh medium (Gomez-Perez et al., 2014). Late-stage promastigotes from L. donovani lines were used to infect macrophages at a macrophage/parasite ratio of 1:10. Twenty-four h after infection at 35 °C and 5% CO2, extracellular parasites were removed by washing with phosphate buffered saline (PBS; 1.2 mM KH2PO4, 8.1 mM Na2HPO4, 130 mM NaCl, and 2.6 mM KCl, pH 7). Infected macrophage cultures were maintained in RPMI 1640 medium plus 10% HIFBS at 37 °C with 5% CO2 at different PMM concentrations for 72 h. Afterwards, macrophages were fixed for 30 min at 4 °C with 2.5% paraformaldehyde in PBS and permeabilized with 0.1% Triton X-100 in PBS for 30 min. Intracellular parasites were detected by nuclear staining with Prolong Gold antifade reagent plus DAPI (Invitrogen). Drug activity was determined from the percentage of infected cells and the number of amastigotes per cell in drug-treated versus non-treated cultures (Seifert et al., 2007).

Measurement of mitochondrial membrane potential (ΔΨm) in Leishmania lines

ΔΨm was measured by flow cytometry using Rh123 accumulation as described previously (Garcia-Hernandez et al., 2012). Depolarization of ΔΨm is indicated by a reduced accumulation of the dye. Parasites (107 promastigotes/ml) were washed twice with PBS and incubated in HBS buffer (21 mM HEPES, 0.7 mM Na2HPO4, 137 mM NaCl, 5 mM KCl, and 6 mM glucose, pH 7) with 0.5 μM Rh123 for 15 min at 28 °C. Then they were washed twice with PBS, resuspended in PBS, and analyzed by flow cytometry in a FACScan flow cytometer (Becton-Dickinson, San Jose, CA) equipped with an argon laser operating at 488 nm. Fluorescence emission between 515 and 545 nm was quantified using the Cell Quest software.

Measurement of reactive oxygen species (ROS) production

The generation of intracellular ROS was measured using the cell-permeable nonfluorescent probes dichlorodihydrofluorescein diacetate (H2DCFDA) and MitoSOX red, as described previously (Berg et al., 2015). Briefly, parasites (107/ml) of the different L. donovani HU3 lines were washed twice with PBS, and either directly analyzed (control untreated) or further incubated in PBS for 2 h at 28 °C (nutritional stress) before analyzing. Next, the promastigotes were incubated in the same buffer supplemented with 4 μM H2DCFDA or 5 μM of MitoSOX for 30 min at 28 °C. The fluorescent dichlorofluorescein (DCF) product obtained after esterase cleavage and oxidation was measured by flow cytometry using a FACScan flow cytometer (Becton-Dickinson, San Jose).

Statistical analysis

Statistical comparisons between groups were performed using Student's t-test. Differences were considered significant at a level of p < 0.05. Also, when it cannot be assumed that data followed a normal distribution, the non-parametric Wilcoxon–Mann–Whitney test was used. A value of p < 0.05 was considered statistically significant.

Accession numbers

RNA and DNA sequencing data have been deposited in the European Nucleotide Archive with the accession number PRJEB20187.

Results and discussion

Experimental design, global transcriptomics changes and chromosomal alterations

In a previous work using a step-wise adaptation process (Garcia-Hernandez et al., 2012), we established lines derived from L. donovani (HU3 strain) with increased resistance to the four clinically relevant antileishmanial drugs for VL: SbIII, AmB, MIL and PMM. Since the establishment of the different lines (Garcia-Hernandez et al., 2012), they were stored in cryobanks. For this study, the lines were thawed and grown in the presence of the appropriate drug concentrations during 30 passages before proceeding with nucleic acids extraction. For each passage, promastigote cultures were diluted in fresh medium at 1 × 106 parasites/ml and grown for 4 days until the cultures arose the stationary phase (cell density: 3-3.5 × 107 parasites/ml). The concentrations of drugs used during culturing of the resistant lines were: 80 μM SbIII (line S), 0.1 μM AmB (line A), 8 μM MIL (line M) and 20 μM PMM (line P). All resistant lines present a similar growth curves to the parental strain HU3 (WT), supporting that drug pressure was not affecting the viability of parasites (Garcia-Hernandez et al., 2012).

Three independent biological replicates of each resistant line and the WT strain, each originating from a separate growth, were processed for RNA and DNA isolation. Illumina RNA-Seq data were independently obtained for the poly-A+ enriched RNA fractions derived from each one of the three replicates. In-depth genomic DNA sequencing was carried out for replicate one from WT strain and the four resistant lines (A, M, P and S). Table 1 shows the sequence reads that were produced across the 20 samples, and the percentages of reads mapped to the reference genomes for L. donovani and L. infantum. In all samples, it was found a slightly larger number of reads mapping to the L. infantum (JPCM5) reference genome than to the L. donovani (BPK) genome. These results are not surprising given that L. infantum and L. donovani species have diverged very recently and are genetically almost indistinguishable (Lukes et al., 2007). Moreover, even though both reference genomes are not complete, the L. infantum reference genome (Peacock et al., 2007) is currently more complete than the L. donovani BPK reference genome (Downing et al., 2011), and therefore a larger number of reads mapped to the L. infantum reference genome. Additionally, the transcriptome of L. infantum has been delineated (version 1; http://leish-esp.cbm.uam.es/l_infantum_downloads.html). Hence, we decided to use the L. infantum JPCM5 genome as the reference for the subsequent analyses.

Table 1

DNA- and RNA-seq data and mapping summary.

Sample	Number of readsa	Aligned vs LdBPKb	Aligned vs LinJPCb
RNA_HU3_wt_1	18,658,148	94.97	95.92
RNA_HU3_A_1	23,474,747	94.89	95.94
RNA_HU3_M_1	15,473,695	94.87	95.79
RNA_HU3_P_1	19,785,684	94.87	95.84
RNA_HU3_S_1	18,837,447	95.12	96.14
RNA_HU3_wt_2	18,573,421	94.90	95.86
RNA_HU3_A_2	16,205,678	94.38	95.39
RNA_HU3_M_2	18,979,494	94.93	95.86
RNA_HU3_P_2	19,885,687	94.96	95.91
RNA_HU3_S_2	17,406,090	94.77	95.76
RNA_HU3_wt_3	19,274,987	95.22	96.18
RNA_HU3_A_3	18,511,557	94.94	95.97
RNA_HU3_M_3	19,042,397	95.11	96.03
RNA_HU3_P_3	21,501,077	94.91	95.87
RNA_HU3_S_3	19,868,776	95.12	96.12
DNA_HU3_wt_1	16,980,871	66.70	68.92
DNA_HU3_A_1	15,550,019	58.82	61.05
DNA_HU3_M_1	13,465,664	61.04	63.37
DNA_HU3_P_1	13,562,485	60.76	62.87
DNA_HU3_S_1	13,026,200	68.89	71.29

The size of reads were: 2 × 76 bp for RNA samples and 2 × 101 bp for DNA ones.

Reference genomes for L. donovani (strain BPK282A1) and L. infantum (strain JPCM5) were retrieved from TriTrypDB database (www.tritrypdb.org, release 9.0).

RNA-Seq reads mapping to the L. infantum reference genome were assigned to L. infantum transcripts and quantified using HTseq (0.6.1 (Anders et al., 2015);) and the data analyzed by DESeq2 (Love et al., 2014) as detailed in Methods section. Transcripts showing a q-value <0.01 and FC differences higher than 1.5 regarding the WT line were categorized as differentially expressed for each one of the drug-resistant lines (S, A, M and P). A total of 1006 differentially expressed transcripts were identified in the S line, 379 in the A line, 146 in the M line, and 129 in the P line. A detailed list of the transcripts differentially expressed in each one of the resistant lines is provided in the Supplementary File 1. When the FC differences between the resistant lines and the reference strain is fixed as 2 or above (Fig. 1), the numbers of differentially expressed transcripts decrease substantially: 194 in the S line, 68 in the A line, 18 in the M line and 13 in the P line. Most differentially expressed transcripts corresponded to transcript upregulated in the resistant lines regarding the WT line. Venn diagrams with the numbers of transcripts having significantly altered levels (either up- or down-regulated) showed the existence of a remarkable number of concurring transcripts when the A and S lines are compared, particularly among the up-regulated transcripts: 38 out of the 64 transcripts up-regulated in the A line were found up-regulated also in the S line. However, as discussed below, many of coincidental data might be derived from somy changes for particular chromosomes observed in these resistant lines regarding the parental strain. In fact, when Euclidean distance heatmap analyses were used to visualize the relationship between samples based on the relative expression of all the transcripts (Fig. 1B), a clear separation among the different lines was observed; in particular, the lines S and A are clearly separated.

Fig. 1

General analysis of differentially expressed transcripts in the amphotericin B (A), miltefosine (M), paromomycin (P) and antimonial (S) lines. A) Venn diagram showing transcripts down- and up-regulated more than 2-fold in the resistant lines regarding the expression levels found in the parental WT line. The Venn diagram was constructed using the online software Venn diagrams, available at Web site of the Bioinformatics and Evolutionary Genomics group of the University of Gent and the VIB institute (http://bioinformatics.psb.ugent.be/webtools/Venn/). B) Heatmap of a hierarchical clustering analysis based on the data from three individual replicates. Euclidean distances were calculated over rlog transformed counts using DESeq2 and plotted with pHeatmap R package (https://cran.r-project.org). C) Volcano plots build using DESeq2 results. Dots in red represent differentially expressed genes with abs(log2(FC)) > 1 and adjusted p-value (Q value) < 0.01. To avoid calculation errors, we have assigned the minimum adjusted p-value higher than 0 to those transcripts with adjusted p-value equal to 0. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

On the other hand, we used alignments of Illumina genomic reads against the reference L. infantum genome to assess differences in chromosomal copy number among the different resistant lines and the WT strain (Fig. 2). After normalization, taking into account the total number of reads obtained for each line (see Materials and Methods for further details), the median read depth (coverage) of the chromosomal set for each line was considered as the diploid status (solid line in Fig. 2). In all cases (WT and resistant lines), chromosome 31 showed the greatest reads coverage, and the data suggested that this chromosome would be pentasomic in all samples here studied. It has been reported previously that this chromosome is at least tetrasomic in L. major (Akopyants et al., 2009) and other Leishmania species (Imamura et al., 2016). Furthermore, chromosome copy number variations compared to the WT somy were observed in the different resistant lines. Thus, in the S line, two extra copies for the chromosomes 1 and an extra copy for chromosome 29 were observed, whereas the somy of chromosome 23 decreased from three (WT) to two (or lover) in the S line (Fig. 2A–B). In the A line, the copy number of chromosome 29 was found to be one-fold higher than that in the WT strain (Fig. 2C). In the M line, the chromosome 5 achieved the status of pentasomic, even though this chromosome is at least trisomic in the WT strain (Fig. 2D). Finally, in the P line the somy for chromosomes 5 and 8 was found to be one-fold higher than in the WT strain (Fig. 2E). With the exception of chromosome 29, whose amplification was found to be partial (see below), the amplification of the other chromosomes seemed to be total as the read depth coverage was even along the whole molecules.

Fig. 2

Chromosome copy number variation in The somy estimation was done using a 2-loop method (Mondelaers et al., 2016). The median coverage of the genome is shown by a solid line, and it was assigned as 2, taking into account that diploid is considered the major ploidy status in Leishmania. The dotted lines indicated the estimated values for other somies. Graphs were generated from the median coverage values by barplot function of R package (https://cran.r-project.org).

Differentially expressed transcripts in the L. donovani S line

Given the wide use of antimonials for leishmaniasis treatment, great efforts are being done in order to understand their damaging effects on the parasites and the underlying mechanisms of resistance developed by Leishmania. Results support the existence of many cellular targets for these drugs, and a variety of mechanisms of resistance. Thus, several metabolic pathways and membrane transporters have been implicated in the resistance phenotype to antimonials (Laffitte et al., 2016b). In agreement with this broad spectrum of responses elicited by this drug, our study showed a tremendous alteration at the transcriptomic level in the S line. Thus, the levels of 194 transcripts were found to be significantly altered (with FC higher than 2) in the resistant line regarding the parental L. donovani HU3 strain (Fig. 1). Nevertheless, it is likely that a fraction of these changes are the result of the chromosomal alterations observed in the resistant line (Fig. 2).

Active efflux of the drug has been suggested as a resistance mechanism to antimonials (Dey et al., 1994). In this regards, it is noteworthy the up-regulation in the S line of two transcripts (LinJ.31.T3190 and LinJ.28.T2050) that are annotated as coding for putative transporters for iron/zinc and zinc, respectively (Table 2). Another up-regulated transcript coding also for a transporter was LinJ.06.T1320, which encodes for a pteridine transporter.

Table 2

Top 30 differentially expressed transcripts in the L. donovani line resistant to trivalent antimony (S line).

Transcript	Log2FC (WT/A)	Producta
LinJ.35.T2640	1.65	hypothetical protein (ubiquitin supergroup?)
LinJ.36.T5110	1.62	hypothetical protein, conserved
LinJ.36.T5100	1.53	hypothetical protein, conserved
LinJ.01.T0850	1.48	peptidyl dipeptidase (pseudogenic transcript)
LinJ.29.T3030	−1.38	amastin, putative
LinJ.34.T1150	−1.38	amastin-like surface protein, putative
LinJ.36.T4910	−1.40	cytochrome b5 1, putative
LinJ.29.T2010	−1.41	hypothetical protein, conserved
LinJ.06.T0890	−1.42	dihydrofolate reductase-thymidylate synthase (DHFR-TS)
LinJ.29.T2020	−1.43	hypothetical protein, conserved
LinJ.31.T3190	−1.45	iron/zinc transporter protein-like protein
LinJ.10.T0590	−1.47	hypothetical protein
LinJ.29.T1020	−1.53	A-1 protein, putative
LinJ.30.T2590	−1.55	formate--tetrahydrofolate ligase (FTHS)
LinJ.18.T0430	−1.61	hypothetical protein
LinJ.26.T2600	−1.62	protein kinase, putative
LinJ.29.T1450	−1.62	amastin-like protein
LinJ.29.T1210	−1.79	hypothetical protein, conserved
LinJ.06.T1320	−1.87	pteridine transporter, putative
LinJ.07.T1320	−1.88	peptide methionine sulfoxide reductase-like, putative
LinJ.29.T1610	−1.95	hypothetical protein, conserved
LinJ.31.T0340	−1.97	hypothetical protein, conserved
LinJ.14.T0470	−1.99	cystathionine beta-lyase-like protein
LinJ.29.T3000	−2.02	amastin, putative
LinJ.31.T0010	−2.27	5-methyltetrahydropteroyltriglutamate-homocysteine S-methyltransferase, putative
LinJ.36.T6670	−2.29	methylenetetrahydrofolate reductase, putative
LinJ.29.T3010	−2.34	amastin, putative
LinJ.29.T0930	−2.40	hypothetical protein, conserved
LinJ.28.T2050	−2.50	zinc transporter, putative
LinJ.13.T1430	−4.11	Protein Associated with Differentiation, putative

Hypothetical protein: predicted bioinformatically; conserved: predicted protein presents also in Trypanosoma brucei and/or T. cruzi (GeneDB).

In addition to those transcripts shown in Table 2, there are other differentially expressed transcripts that match with data found in previous reports, reinforcing their possible association with resistance to antimonials. Thus, Singh and co-workers demonstrated that over-expression of histone H2A increases resistance of L. donovani parasites to sodium antimony gluconate (SAG) (Singh et al., 2010). In agreement with that finding, we found that one of the H2A-coding transcripts (i.e. LinJ.29.T1870) showed a 2-fold over-expression in the S line. Also, the Salotra's group showed that over-expression of PSA-2 in L. donovani results in development of resistance to SAG (Bhandari et al., 2013). In agreement, we found that four transcripts coding for PSA-2 (LinJ.12.T0662, LinJ.12.T0665, LinJ.12.T0666, LinJ.12.T0690) were about 2-fold more abundant in the resistant line (see Supplementary File1). Gamma-glutamylcysteine synthetase (γGCS, encoded by gene GSH1) is the rate-limiting enzyme in the biosynthesis of thiols glutathione and trypanothione in Leishmania, and these thiols are considered to have a critical role in neutralizing the oxidative stress induced by antimonials. In fact, overexpression of the gene GSH1 in Leishmania leads to an increased resistance to arsenite (Grondin et al., 1997). Conversely, Leishmania heterozygous mutants, having only one GSH1-allele, synthesized less glutathione and trypanothione, and consequently showed a significant decreased survival in the presence of SbV compared with control cells (Mukherjee et al., 2009). In a recent article, it was shown that L. guyanensis parasites overexpressing GSH1 presented an increase of 4-fold in SbIII-resistance index when compared with the wild-type line (Fonseca et al., 2017). Interestingly, our data showed that transcript LinJ.18.T1660 (encoding for γGCS) was about 2.5-fold more abundant in the resistant S line than in the WT parasites. Furthermore, other transcripts encoding enzymes that are part of the biosynthesis pathway of trypanothione were found to be significantly up-regulated in the resistant line: (i) LinJ.30.T3580, which codes for the S-adenosylmethionine synthetase (METK2), was 1.8-fold more abundant in the S line; (ii) transcript LinJ.31.T0010, encoding for a methionine synthase, was 4.8-fold over-expressed in line S; (iii) transcript LinJ.06.T0890 (dihydrofolate reductase-thymidylate synthase, DHFR-TS) was 2.7-fold more abundant also in the S line; and (iv) transcript LinJ.36.T6670 (methylenetetrahydrofolate reductase) that showed a 4.9-fold increased level in the S line. Together, these data support the idea that resistance to antimonials may be linked to a general antioxidant response in the parasite. In this regards, there are several studies relating metacyclogenesis (differentiation process essential for infectivity) to antimony resistance. Thus, in L. donovani, it was established that the capacity of differentiation to metacyclic forms was significantly higher in antimony-resistant cells (Ouakad et al., 2011). Another work has documented that clinical isolates showing resistance to antimonials attained stationary phase at a higher parasite density, contained a higher amount of metacyclic parasites and had a greater capacity to cause in vivo infection in mice compared to the antimonial-sensitive strains (Vanaerschot et al., 2010). Furthermore, it was found that antimonial-resistant parasites cause a significant higher in vivo parasite load compared to antimonial-sensitive parasites (Vanaerschot et al., 2011).

Among the top differentially expressed transcripts, it is remarkable the presence of 5 up-regulated transcripts in the S line encoding for members of the amastin family (LinJ.29.T1450, LinJ.29.T3000, LinJ.29.T3010, LinJ.29.T3030, and LinJ.34.T1150). Transcripts encoding amastins were also found to be upregulated in the A and P lines (see below).

On the other hand, among the 30 transcripts showing highest differences in their steady-state levels between the parental and S lines (Table 2), only four were down-regulated in the resistant line. Three of them (LinJ.35.T2640, LinJ.36.T5110 and LinJ.36.T5100) encode for hypothetical proteins without putative function, and the fourth (LinJ.01.T0850) is annotated as coding for a metallopeptidase. Nevertheless, according to the GeneDB annotation, the latter may be a pseudogenic transcript. However, it is curious that the levels of this transcript (LinJ.01.T0850) remained 2.8-fold below that found in the WT strain, despite the mutant S line has two-fold more copies of chromosome 1 than the WT strain (Fig. 2B).

The transcript showing the greatest difference in expression was LinJ.13.T1430, whose steady-state level was found to be 17.3-fold higher in the resistant S line than in the parental HU3 strain. The corresponding gene is annotated in the GeneDB database as “protein associated with differentiation” (PAD), because its sequence conservation with PAD1 from T. brucei. This protein, a member of the carboxylate-transporter family, was found to be specifically expressed on the surface of the transmission-competent ‘stumpy-form’ parasites in the bloodstream as a prerequisite to initiate life-cycle development when transmitted to their tsetse fly vector (Dean et al., 2009). However, to our knowledge, the function of this protein in Leishmania was not yet experimentally explored.

In order to confirm the results obtained by RNA-Seq, the expression of the PAD gene in the S-resistant line was analyzed. For this purpose, total RNA from WT and S-resistant lines was extracted and the levels of expression of PAD gene were determined by RT-qPCR; this analysis showed that the PAD gene was expressed 2.1-fold higher in the S-resistant line than in the WT line (Fig. 3A). It was unexpected to find these differences in the magnitude of the fold change determined for this gene between the RNA-Seq and RT-qPCR methodologies. Finally, we found the explanation in the fact that the L. infantum genome contains two additional genes coding for PAD (LinJ.11.0660 and LinJ.11.0670) apart from LinJ.13.1430. All three genes are conserved in the coding regions but are divergent at their 3′UTRs. The RT-qPCR was done using oligonucleotides derived from the coding region (common to all three genes), whereas RNA-Seq reads are assigned to the individual transcripts. In the S line, only transcript LinJ.13.T1430 is differentially expressed.

Fig. 3

Gene expression analysis of Total RNA of L. donovani WT line, S, A and P resistant lines, and WT + Lexsy, A + Lexsy, WT + PAD, A + SMT, WT + D-LDH, WT + BCAT, WT + D-LDH + BCAT and WT + BCAT + D-LDH transfected lines, was extracted from promastigotes grown at log-phase as described in Materials and Methods. RT-qPCR expression values of the genes in each line were normalized with the expression of GAPDH. The relative expression of each gene was calculated as the fold-change between the resistant lines and the WT line (which was set to 1.0) (A) or transfected lines with PAD, SMT or D-LDH, BCAT, and WT + Lexsy, WT + Lexsy + pIR1SAT or A + Lexsy (set to 1.0) (B), see Material and Methods for further details. Results shown are the means ± SD from two independent experiments. In all cases, significant differences versus the controls were determined using Student's t-test (p < 0.01).

On the other hand, to directly analyze the possible contribution of this protein to the phenotype of antimony resistance, the PAD gene was cloned into pLEXSY-hyg2 expression vector and transfected into the WT line (WT + PAD line). As control, the WT line was transfected with the empty vector, creating the WT + Lexsy line. The over-expression of the corresponding transcripts in both transfected lines was determined by RT-qPCR, the results showed an increase in PAD expression levels of 18.8-fold regarding the control WT + Lexsy line (Fig. 3B). In parallel, drug susceptibility of WT, S-resistant, WT + PAD, and WT + Lexsy lines to SbIII was determined in L. donovani promastigotes (Table 3). The resistance index (RI) of the WT + PAD line was found to be only 1.5, significantly lower than the 3.7 RI value obtained for the S-resistant line (Table 3), suggesting that the PAD protein would be contributing slightly to the resistance phenotype in this resistant S line.

Table 3

Drug susceptibility in the different L. donovani lines used in this study.

Lines	SbIII EC50 μM (RI)a	AmB EC50 μM (RI)	PMM EC50 μM (RI)
WT	89.12 ± 9.67b	0.09 ± 0.01	14.14 ± 2.36
S	331.21 ± 39.08(3.72)	-c	–
A	–	0.23 ± 0.04(2.64)	–
P	–	–	132.37 ± 8.26(9.36)
WT + Lexsy	79.62 ± 12.56	–	11.98 ± 0.63
A + Lexsy	–	0.31 ± 0.02(3.56)	–
WT + PAD	117.96 ± 14.11(1.48)	–	–
A + SMT	–	0.12 ± 0.03(1.35)	–
WT + D-LDH	–	–	58.36 ± 3.07(4.87)
WT + BCAT	–	–	48.94 ± 4.52(4.08)
WT + Lexsy + pIRSAT	–	–	13.92 ± 2.36
WT + D-LDH + BCAT	–	–	56.76 ± 6.25(4.08)
WT + BCAT + D-LDH	–	–	66.24 ± 3.39(4.76)

Resistance indexes (RI) were calculated by dividing the EC50 for each resistant line, or overexpressing D-LDH, BCAT, PAD or SMT lines, by that for Leishmania control line (WT for A, P and S lines, WT + Lexsy or A + Lexsy for the transfected lines, and WT + Lexsy + pIRSAT for cotransfected lines).

Data are the averages ± SD of EC50 values from six independent experiments. Significant differences versus the controls were determined using Student's t-test (p < 0.005).

-: Not determined.

Identification of transcripts differentially expressed in the AmB resistant line (A line)

In the A line, among the differentially expressed transcripts regarding the parental line, 379 transcripts showed differences in the expression levels above 1.5-fold (supplementary file 1), and for 68 of them the FC was higher than 2 (Fig. 1). The top 24 down- and upregulated transcripts are shown in Table 4. Most of them, 20 out of the 24, correspond to transcripts whose levels were increased in the A line. Another striking finding was that 19 of them correspond to the expression of genes located at chromosome 29. Moreover, when all the differentially expressed transcripts are considered, 271 out of the 379 transcripts (71.5%) showing differential expression in the A line derive from loci located at chromosome 29. It is likely that this enrichment in transcripts derived from chromosome 29 is related to the fact that the ploidy of this chromosome in the A line was found to be two-fold higher than of the parental line (Fig. 2C). In particular, as shown in Fig. 4, the amplification observed in the A line affected to 2/3 of the chromosome 29, from position 341113 to the end of the chromosome. Remarkably, as shown also in Fig. 4, similar amplification of this large region of chromosome 29 was also observed in S line and in a lower extent in the P line. In all three cases, an abrupt change in coverage was mapped in the divergent strand switch region (dSSR) located between genes LinJ.29.0960 and LinJ.29.0970. At present, the molecular mechanism responsible for this partial amplification of this chromosome could not be depicted; the SIDER2 repeated sequences found in the vicinity (Fig. 4B) did not seem to be involved in the rearrangement. However, we found in this dSSR several sequences with potential to conform G-quadruplexes, and given the constraints that these structures impose in DNA replication (Valton and Prioleau, 2016), they may be chromosomal fragile sites. This would explain that the same chromosomal alteration is observed in the three resistant lines (A, S and P; Fig. 4B).

Table 4

Top 24 differentially expressed transcripts in the L. donovani AmB resistant line (A line).

Transcript	Log2FC (WT/A)	Product
LinJ.36.T2510	2.86	sterol 24-c- methyltransferase, putative
LinJ.36.T5110	1.26	hypothetical protein, conserved
LinJ.36.T5100	1.21	hypothetical protein, conserved
LinJ.23.T0700	1.13	hypothetical protein
LinJ.29.T1260	−1.12	hypothetical protein, conserved
LinJ.29.T1550	−1.13	phosphatidylinositol 4- kinase alpha, putative
LinJ.29.T1090	−1.14	hypothetical protein, conserved
LinJ.29.T2820	−1.14	hypothetical protein, conserved
LinJ.29.T1460	−1.15	RNA binding protein, putative
LinJ.29.T1880	−1.15	paraflagellar rod protein 1D, putative
LinJ.29.T2950	−1.15	hypothetical protein, conserved
LinJ.29.T2400	−1.16	hypothetical protein
LinJ.26.T2600	−1.17	protein kinase, putative
LinJ.29.T1600	−1.17	hypothetical protein, conserved
LinJ.29.T1020	−1.19	A-1 protein, putative
LinJ.29.T1630	−1.19	hypothetical protein, conserved
LinJ.29.T1890	−1.21	paraflagellar rod protein 1D, putative
LinJ.29.T3030	−1.27	amastin, putative
LinJ.29.T2010	−1.30	hypothetical protein, conserved
LinJ.29.T1210	−1.39	hypothetical protein, conserved
LinJ.29.T1450	−1.41	amastin-like protein
LinJ.29.T1610	−1.50	hypothetical protein, conserved
LinJ.29.T3000	−1.77	amastin, putative
LinJ.29.T3010	−2.07	amastin, putative

a Hypothetical protein: predicted bioinformatically; conserved: predicted protein presents also in Trypanosoma brucei and/or T. cruzi (GeneDB).

Fig. 4

Read coverage along the chromosome 29 in WT and A, S and P resistant lines. The coverage values obtained with genomeCoverageBed (BEDTools) were smoothed by an in-house Perl script that calculates the mean value within overlapping 5-Kb windows. Graphs were generated from the coverage values by plotly utility of R package (https://cran.r-project.org). Panel A shows the coverage along the complete chromosome 29, whereas panel B shows the region comprised between the coordinates 322,447 and 390,000 of this chromosome. Green-colored boxes denote ORF (the gene names are indicated), and red-boxes correspond to repeated sequences of the SIDER2 family (Requena et al., 2017).

The highest FC between the AmB resistant line and the parental strain corresponded to the transcript LinJ.36.T2510 (Table 4), which is about 7-fold more abundant in WT strain than in the resistant A line. This transcript codes for 24-sterol methyltranferase (SMT), an enzyme that methenylates steroids at the 24 position during ergosterol biosynthesis. Studies in L. major has shown that this enzyme is mainly localized in the endoplasmic reticulum (Jimenez-Jimenez et al., 2008). The postulated mode of action for AmB suggests that the compound forms complexes with 24-substituted sterols, such as ergosterol in the cell membrane, and then causes pores which alter ion balance and result in cell death (Croft and Coombs, 2003). Hence, taking into account the specific interaction between AmB and ergosterol, a decrease in the activity of the enzyme SMT in the A line would be a plausible resistance strategy to elude the toxicity of this polyene antimycotic. Interestingly, in a previous work dealing with the mechanism of AmB resistance in L. donovani promastigotes, it was found that the major sterol present in the membrane of resistant parasites was an ergosterol precursor, the cholesta-5, 7, 24-trien-3b-ol and not ergosterol, which is the most abundant sterol present in the membrane of the AmB-sensitive strain; therefore, it was concluded that this AmB-resistant Leishmania was defective in C-24 transmethylation of C-27 sterols (Mbongo et al., 1998).

Analysis by whole genome sequencing of the locus Ld-SMT in the WT strain and in the A line showed a sharp decrease in coverage in the A line, suggesting a possible gene deletion event in the A resistant line (Fig. 5A). In the genome of L. infantum, and other Leishmania species, there are two tandemly organized genes (LinJ.36.2510 and LinJ.36.2520) showing an absolute sequence identity in their coding regions, but different in their untranslated regions. Additional experimental evidence of a deletion event occurring in the A line was obtained after PCR amplification of the locus using two oligonucleotides designed to amplify the genomic region containing genes LinJ.36.2510 and LinJ.36.2520. As shown in Fig. 5B, an amplification product of approx. 1.6-kb long was obtained from genomic DNA of the A line, but no when using genomic DNA of the parental strain as template. In contrast, a slight amplification band around 5.5-kb was observed only in the PCR from WT DNA. The latter is the expected size for the amplification product of the complete region according to the L. infantum JPCM5 reference genome. After direct sequencing of the 1.6-kb amplification product, specific for the A line, it was confirmed that a deletion event took place in the resistant line leading to the loss of one SMT gene and the intercoding region separating the orthologues to genes LinJ.36.2510 and LinJ.36.2520. Arguably, this deletion was the result of an intrachromosomal homologous recombination event between the coding regions of genes LinJ.36.2510 and LinJ.36.2520. Hence, as a result, the A line would lose the LinJ.36.T2510 transcript, and, predictably, would have a lower expression of SMT enzyme, and the subsequent decrease of ergosterol levels at the parasite cell membranes. In agreement with this hypothesis, a study by Pourshafie et al. (2004) showed the existence in L. donovani of two transcripts encoding for SMT, named SCMT A and SCMT B; interestingly, SCMT A was found to be absent in an AmB-resistant L. donovani line, selected by stepwise drug pressure. More recently, an L. donovani clinical isolate obtained from an AmB-nonresponsive patient was investigated to understand the mechanism of AmB resistance (Purkait et al., 2012). In that study, it was found that the transcript SCMT A was absent in AmB-resistant parasites, but expressed in susceptible parasites.

Fig. 5

Deletion of a Read depth distribution of Illumina reads from WT and A in the region 948,788–962,094 of chromosome 36. An in-house Perl script was used for representing the smooth coverage with a 200 nucleotide-window and a normalization by the total number of reads from each sample. In the bottom, a graph showing the ORF of the genes located in this region is included. The ORFs of genes LinJ.36.2510 and LinJ.36.2520 are identical and code for SMT. The location and sizes of transcripts LinJ.36.T2510 and LinJ.36.T2520 are indicated by arrows with rhomboid heads. The slim arrows indicate the position of the oligonucleotides used for PCR amplification of the region. (B) PCR amplification with oligonucleotides gSMT-Fw and gSMT-Rv using as template genomic DNA from WT strain or A line. The green arrow points to the expected amplicon in the WT line, whereas the red one indicates the amplicon observed in the A line. HindIII-digested DNA of Φ29 phage was used as molecular weight marker (M). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

To validate that SMT genes present altered levels of expression in the AmB-resistant line, as deduced from the RNA-Seq data, total RNA of WT and A lines was extracted and the relative amount of SMT transcripts was analyzed by RT-qPCR. As shown in Fig. 3A, and in agreement with the RNA-Seq data, the expression levels were 3.9-fold lower in the A line than in the WT line. Next, we asked whether overexpression of SMT in the A line might restore the susceptibility of this line to AmB. For this purpose, the SMT gene was cloned into the pLexsy-hyg2 expression vector, and the construct used to transfect the AmB-resistant line, the resulting line was named A + SMT. In parallel, as control, the AmB-resistant line was transfected with the empty pLexsy-hyg2 vector to create the line A + Lexsy. The results obtained by RT-qPCR showed a 75.2-fold over-expression of SMT gene in the A + SMT line compared to the control (Fig. 3B). Finally, the susceptibility to AmB of WT, A-resistant, A + SMT, and A + Lexsy transfected lines was determined in promastigotes. As shown in Table 3, when the A line over-expressed the SMT gene (A + SMT), it lost the resistance to AmB, showing similar EC50 values for AmB that the WT line (0.12 and 0.09 μM, respectively). On the contrary, the control line (A + Lexsy) maintained the resistance to AmB and exhibited EC50 values for AmB 3.6-fold higher than the WT line. Altogether, these data suggest that the lower SMT expression in the A line respect to the WT strain is directly related to the loss of one SMT gene. Therefore, these results support the use of SMT as a drug resistance marker for AmB resistance. Thus, a PCR based on the amplification of the SMT locus (as depicted in Fig. 5B) might be useful to quickly check SMT deletion events in AmB-resistant, clinical isolates.

On the other hand, the amplification of the most part of chromosome 29, occurring in the A line, suggests that expression of some gene(s) located at this chromosome may be contributing to the phenotype of resistance to AmB. However, the identification of that/those driver gene(s) is not a simple task, taking into account that 271 transcripts derived from this chromosome show increased levels regarding the expression levels in WT L. donovani HU3 (Supplementary file 1). It is noticeable that among the most upregulated transcripts, three are annotated as coding for amastin (LinJ.29.T3000, LinJ.29.T3010, LinJ.29.T3030) and an additional overexpressed transcript (LinJ.29.T1450) is annotated as an amastin-like protein. Amastin genes were first described in Trypanosoma cruzi because of their specific expression at the amastigote stage of this parasite (Teixeira et al., 1994). Subsequently, these genes were identified in Leishmania, and genomic analysis showed that the amastin repertoire is much larger in Leishmania relative to Trypanosoma spp. (Rochette et al., 2005; Jackson, 2010). Thus, in the genomes of L. major and L. infantum up to 45 amastin genes were identified, comprising the largest gene family reported so far in Leishmania (Rochette et al., 2005). However, it was somewhat remarkable that none of the four genes, identified as up-regulated in the A line, was catalogued as coding for amastin in those previous studies. Hence, we analyzed the information regarding these four genes in the GeneDB database. Orthologues to these four genes do not exist in the L. major genome. These genes are highly conserved in sequence each other; thus, the protein encoded by gene LinJ.29.3000 shares 98, 94 and 61% of sequence identity with the proteins encoded by genes LinJ.29.1450, LinJ.29.3010 and LinJ.29.3030, respectively. In contrast, the sequence identity between the LinJ.29.3000-encoded protein and that encoded by gene LinJ.08.0820, a prototypical amastin-gene (Rochette et al., 2005), is only 49%. Although the function of amastin is unknown, its structure is typical for a membrane protein, containing four hydrophobic transmembrane domains, interspersed with serine–threonine rich, extracellular domains and probable glycosylation sites (Teixeira et al., 1994; Rochette et al., 2005). Taking into account a possible location of these molecules in the parasite membrane, and given the interaction of AmB with ergosterol, the predominant sterol in the membrane of Leishmania (Croft et al., 2006), it can be postulated that an increased expression of this class of amastins may be hampering the interaction of the drug with the membrane in the Leishmania resistant A line. Studies based on the overexpression of these genes would serve to test this hypothesis.

Among other transcripts also showing increased levels in the A line regarding the WT strain is LinJ.29.T1250 (see supplementary file 1), which codes for tryparedoxin 1 (TXN1). Interestingly, in a recent study, the protein expression levels of this antioxidant molecule was found to be 2–3 fold higher in AmB resistant isolates as compared to sensitive strains of L. donovani (Suman et al., 2016). Thus, our data reinforce the hypothesis that over-expression of TXN1 might be linked to AmB resistance.

Transcripts differentially expressed in a MIL resistant line (M line)

Eighteen transcripts were found to have above 2-fold expression differences between the M line and the parental parasites (Fig. 1), and Table 5 shows those having annotated genes. Remarkably, most of them are transcripts that result downregulated in the M parasites. MIL resistance mechanisms in Leishmania spp. are associated with a decrease in intracellular drug accumulation due to the defective inward translocation of MIL (Perez-Victoria et al., 2003b, 2006b; Sanchez-Cañete et al., 2009; Mondelaers et al., 2016) and the overexpression of a P-glycoprotein MDR1 (ABCB4), which is responsible for drug efflux (Perez-Victoria et al., 2001, 2006a). Additionally, a member of the ABC subfamily G transporter, LiABCG4, has been found to confer moderate (2-fold) resistance to MIL; this protein has been involved in the outward transport of phosphatidylcholine analogues from the cytoplasmic to the exoplasmic leaflet of the parasite membrane (Castanys-Muñoz et al., 2007). More recently, a Leishmania protein of unknown function, named P299 (LinJ.08.0630), has also been linked to MIL resistance; thus, overexpression of its gene led to a 2.6-fold increased protection against MIL (Choudhury et al., 2008).

Table 5

Protein-coding transcripts that are differentially expressed in the L. donovani miltefosine resistant line (M line).

Transcript	Log2FC (WT/A)	Producta
LinJ.17.T0690	1.84	hypothetical protein
LinJ.26.T2600	1.75	protein kinase, putative
LinJ.14.T0180	1.58	carboxypeptidase, putative,
LinJ.31.T0370	1.41	amino acid transporter AAT1.4
LinJ.17.T0990	1.29	META1
LinJ.23.T1230	1.27	small hydrophilic endoplasmic reticulum-associated protein (SHERP)
LinJ.32.T3810	1.21	3-hydroxyisobutyryl-coenzyme a hydrolase-like protein
LinJ.20.T1370	1.17	kinase-like protein
LinJ.07.T0910	1.13	flavoprotein subunit-like protein
LinJ.28.T1070	1.13	P27 protein (P27)
LinJ.10.T0380	1.00	folate/biopterin transporter, putative
LinJ.33.T1860	−1.06	OTT_1508-like deaminase, putative
LinJ.08.T0180	−1.19	hypothetical protein

Hypothetical protein: predicted bioinformatically; conserved: predicted protein presents also in Trypanosoma brucei and/or T. cruzi (GeneDB).

As LMT/LRos3 complex is pivotal for resistance towards MIL (Perez-Victoria et al., 2003b, 2006a, 2006b), we first looked for possible alterations in the steady-state levels of transcripts LinJ.13.T1590 (LMT) and LinJ.32.T0540 (LRos3). The steady-state levels of both transcripts were slightly higher in the WT strain than in the M line: 1.1-fold for LinJ.13.T1590 and 1.2-fold for LinJ.32.T0540. However, these differences were not found to be statistically significant in this MIL-resistant line.

In a recent study, the dynamics of the appearance of single nucleotide polymorphisms (SNP) in the LMT gene was analyzed during the selection process for MIL resistance in an L. infantum strain (Laffitte et al., 2016a). Interestingly, it was observed an increase in mutated alleles under high drug pressure. Additionally, in a recent report, the existence of point mutations in the LMT gene of two MIL-resistant clinical isolates has been documented (Srivastava et al., 2017). In our MIL resistant line, we examined SNPs in the coding regions of genes LMT and LRos3, comparing sequence reads derived from the WT and M DNA samples. We identified a heterozygous SNP at position 1041 in the LMT gene of the M line: 42% of the reads showed an A and the rest a G residue at that position, whereas no SNP was found in the WT. Interestingly, the G-to-A change leads to the appearance of a stop codon at triplet 347, which in the WT protein codes for Trp. In conclusion, our data suggest that the resistance to MIL in the M line may be associated with a decreased expression of the MIL transporter, leading to a mild alteration in the MIL-inward translocation system. It must be noted that the resistance index of this M line to MIL is 2.5 (data not shown), whereas L. donovani promastigotes, having a non-functional MT subunit, are 15-fold more resistant to MIL than WT parasites (Perez-Victoria et al., 2003a).

While the mechanisms of action of MIL are not yet completely understood, several studies have pointed to alterations in phospholipid metabolism, impairment of bioenergetic metabolism, and recently to the induction of apoptosis as potential modes of action (Paris et al., 2004; Rakotomanga et al., 2007; Mishra and Singh, 2013). From the differential expressed transcripts in the M line, it is not immediately intuitive to postulate a possible role of the encoded proteins in the resistance phenotype of this line. Most of the differentially expressed transcripts appeared as down-regulated in the M line (Table 5), and those that were found to be upregulated showed a modest 2-FC and, furthermore, code for proteins of unknown function. Among the downregulated transcripts, there are two transcripts, LinJ.17.T0990 and LinJ.23.T1230, which code for proteins associated with metacyclogenesis (META1 and SHERP, respectively). Thus, META1 has increased expression at the stationary phase of L. amazonensis, L. donovani and L. major promastigotes (Uliana et al., 1999). The other protein, SHERP (small hydrophilic endoplasmic reticulum (ER)- associated protein), is expressed only in metacyclics, in which it shows association with the ER and outer mitochondrial membrane (Knuepfer et al., 2001).

Among other transcripts down-regulated in the M-line, there are two coding for transporters: LinJ.31.T0370 (coding for the amino acid transporter AAT1.3) and LinJ.10.T0380 (coding for a folate/biopterin transporter). Leishmania contain a large family of folate/biopterin transporters composed by 14 different members, and it is remarkable that one of them, FT1 (encoded by gene LinJ.10.0400) is not expressed in methotrexate-resistant L. infantum (Ouameur et al., 2008). In a recent work, Vacchina and co-workers (Vacchina et al., 2016) analyzed changes in gene expression by RNA-Seq in a MIL-resistant line of L. donovani, and also they found a down-regulation of several folate/biopterin transporters in the resistant line. Additionally, coincident with our data (supplementary file 1), these authors identified, as down-regulated in their MIL-resistant line, the genes coding for an amino acid permease (AAT20.2; LinJ.10.0770), a surface antigen-like protein (LinJ.04.0170) and a hypothetical protein (LinJ.04.0160).

Transcripts differentially expressed associated with the resistance to PMM in the P-line

The mechanism of action of PMM has been well studied in Escherichia coli, where it binds to the 30S ribosomal A-site RNA, and causes misreading of the genetic code and inhibits translocation (Fourmy et al., 1996). Aminoglycosidic antibiotics were initially used against bacterial infection (currently, they are out of clinical use), but subsequently some were found to be effective against several parasites, among others, Leishmania (Davidson et al., 2009). However, the mechanism(s) of action of PMM in Leishmania spp is(are) not known yet (Croft et al., 2006), even though recent studies point to an interaction with the Leishmania ribosome A-site (Shalev et al., 2015). On the contrary, resistance to aminoglycosides in bacteria is well known and it is attained by three ways: (i) decreased uptake; (ii) alteration of the ribosomal binding site, and (iii) chemical modification of the amino and/or hydroxyl groups by bacterial enzymes (Davies and Wright, 1997).

In our study, only 13 transcripts varied above 2-fold in the PMM-resistant line (P line; Fig. 1). Table 6 shows those transcripts having annotated genes. Among them, two transcripts (LinJ.27.T1940 and LinJ.27.T1950) were strongly up-regulated in the resistant line, 56- and 32-fold, respectively, compared to their expression in the parental (WT) line. According to GeneDB annotations, the encoded products for these transcripts are D-lactate dehydrogenase-like protein (D-LDH) and an aminotransferase of branched-chain amino acids (BCAT), respectively. More strikingly, genomic DNA coverage revealed that the region encoding for both transcripts, which are collinear, was over-replicated about 150-fold in the PMM-resistant line (Fig. 6). At present, we have not investigated about the nature of this amplicon, whether intra- or extra-chromosomal, but the presence of 2 direct repeats of 424-nucleotides in length that are located in the intergenic region (Fig. 6B) is noteworthy. The genome of Leishmania contains a significant number (around 2000) of repeated DNA sequences (Ubeda et al., 2014), which provide to this parasite the capacity of generating chromosomal re-arrangements and extrachromosomal DNA elements in response to drugs pressure (Laffitte et al., 2016b).

Table 6

Protein-coding transcripts that are differentially expressed in the L. donovani paromomycin-resistant line (P line).

Transcript	Log2FC (WT/A)	Producta
LinJ.07.T0910	1.06	flavoprotein subunit-like protein
LinJ.29.T1600	−1.00	hypothetical protein, conserved
LinJ.29.T1450	−1.06	amastin-like protein
LinJ.13.T1430	−1.07	Protein Associated with Differentiation, putative
LinJ.29.T0930	−1.24	hypothetical protein, conserved
LinJ.31.T3190	−1.37	iron/zinc transporter protein-like protein (LIT1)
LinJ.29.T3000	−1.79	amastin, putative
LinJ.27.T1950	−4.98	branched-chain amino acid aminotransferase, putative
LinJ.27.T1940	−5.77	D-lactate dehydrogenase-like protein, putative

Hypothetical protein: predicted bioinformatically; conserved: predicted protein presents also in Trypanosoma brucei and/or T. cruzi (GeneDB).

Fig. 6

Amplification in the P line of the genomic region containing the genes LinJ.27.1940 ( (A) The reads coverage along a region (coordinates 853,000–872,000) of the chromosome 27 in both WT and P lines is shown. The coverage values obtained with genomeCoverageBed (BEDTools) were smoothed by an in-house Perl script that calculates the mean value within overlapping 200-nucleotide windows. Note that a logarithmic scale is used to represent the coverage values. (B) Schematic drawing of the amplified genomic region, showing the position of genes LinJ.27.1940 (D-LDH) and LinJ.27.1950 (BCAT). The red boxes correspond to two identical sequences of 424-nucleotides in length. The dotted arrow at the bottom denotes the transcriptional direction of the genes. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

The over-expression of the D-LDH and BCAT genes observed by RNA-Seq in the PMM-resistant line was confirmed by RT-qPCR (46.6- and 19.2-fold of over-expression, respectively, regarding the WT strain; Fig. 3A). In order to evaluate the possible contribution of these genes to PMM resistance in Leishmania, we created the constructs expressing either separately each one of the genes or co-expressing both D-LDH and BCAT genes. As shown in Fig. 3B, a significant increase in the expression of the genes (32.0-fold and 47.93-fold for D-LDH and BCAT, respectively) was observed after transfection of the Leishmania WT line with the D-LDH or BCAT-expressing constructs (lines WT + D-LDH and WT + BCAT, respectively, Fig. 3B). Interestingly, both transfectants acquired a significant resistance against PMM, with RI values of 4.9 and 4.1 for WT + D-LDH- and WT + BCAT-transfected lines, respectively, whereas the RI value of P line was 9.4 (Table 3). Nevertheless, co-expression of both D-LDH and BCAT genes in L. donovani WT line did not produce a synergistic effect, obtaining values of EC50 to PMM in the same range of single overexpression of either D-LDH or BCAT (Table 3). Nevertheless, we could not achieve similar expression levels to the single overexpression of BCAT gene despite using different plasmids and after 13 passages with 100 μg/ml of nourseothricin or Hyg (Fig. 3B).

To gain greater understanding on the molecular basis linking the overexpression of D-LDH and BCAT genes to PMM resistance, we analyzed the resistance of both the P line and the transfected parasites with D-LDH and BCAT to other aminoglycosidic antibiotics such as geneticin (G-418) and nourseothricin. Remarkably, three independent transfection experiments indicated a significant cross-resistance (Table 7), with average RI values to geneticin of 5.6 and 4.5 for WT + LDH- and WT + BCA-transfected lines, respectively, whereas the RI value for P line was 12.1. Regarding nourseothricin, WT + D-LDH and WT + BCAT-transfected lines exhibited average RI values of 1.9 and 2.1, respectively, whereas it was determined an RI value of 3.5 for P line (Table 7). The similar high RI values of the P line and the transfectants against PMM and geneticin might be related to the fact that their structures are very close, whereas the structure of nourseothricin is clearly different. This suggested that the resistance may arise from chemical modification of the amino and/or hydroxyl groups of these antibiotics, and that the putative over-expression of enzymes D-LDH and BCAT may be contributing to the chemical modification of these molecules. Nevertheless, other metabolic possibilities exist and they were further explored (see below).

Table 7

Aminoglycosides susceptibility in L. donovani lines.

Lines	PMM EC50 μM (RI)a	Geneticin EC50 μg/ml (RI)	Nourseothricin EC50 μg/ml (RI)
WT	12.18 ± 1.31b	0.44 ± 0.09	9.26 ± 1.52
P	122.76 ± 7.32 (10.08)	5.34 ± 0.81 (12.14)	32.81 ± 4.16 (3.54)
WT+Lexsy	13.14 ± 1.28	0.42 ± 0.09	8.13 ± 0.59
WT+D-LDH	62.15 ± 4.07 (4.73)	2.36 ± 0.23 (5.62)	15.51 ± 1.52 (1.91)
WT+BCAT	48.79 ± 4.80 (3.71)	1.88 ± 0.20 (4.48)	16.81 ± 1.77 (2.07)

Resistance indexes (RI) were calculated by dividing the EC50 for P line, or overexpressing D-LDH or BCAT lines by that for Leishmania control line (WT for P, and WT + Lexsy for WT + D-LDH and WT + BCAT).

Data are the averages ± SD of EC50 values from six independent experiments for three independent transfections. Significant differences versus the controls were determined using Student's t-test (p < 0.001).

Firstly, it was analyzed if the resistance observed in promastigote forms was also maintained in intracellular amastigotes obtained after infection of THP-1 cells with the different promastigote lines (Table 8). The results showed that the resistance to PMM in intracellular amastigotes was maintained; with RI values of 2.9 and 2.4 for WT + D-LDH and WT + BCAT transfected lines, respectively, whereas the RI for P line was 4.1 (Table 8). These values were 2.78 and 2.79 to the WT + D-LDH + BCAT and WT + BCAT + D-LDH lines, respectively, when the co-transfection with both D-LDH and BCAT genes were done. These results demonstrate that the resistance phenotype to PMM in the amastigote stage of the L. donovani P line would be due in part to the over-expression of D-LDH and BCAT.

Table 8

Paromomycin susceptibility in intracellular amastigotes of L. donovani lines.

Lines	PMM EC50 μM (RI)a
WT	9.49 ± 1.21b
P	38.51 ± 3.83 (4.05)
WT + D-LDH	29.24 ± 3.26 (2.96)
WT + BCAT	24.03 ± 2.64 (2.44)
WT + Lexsy	9.86 ± 1.30
WT + D-LDH + BCAT	36.37 ± 5.81 (2.78)
WT + BCAT + D-LDH	38.47 ± 6.32 (2.94)
WT + Lexsy + pIRSAT	13.09 ± 2.61

Resistance indexes (RI) were determined from the percentage of infected cells and the number of amastigotes per cell in drug-treated cultures versus non-treated cultures.

Data are the mean ± SD of EC50 values from four independent experiments. Significant differences versus the controls were determined using Student's t-test (p < 0.001).

Additionally, we wanted to know whether over-expression of D-LDH or BCAT genes in other Leishmania species/strains also confer resistance to PMM. For this purpose, promastigotes of L. donovani (Dd8 strain), L. major (Friedlin strain), and L. infantum (JPCM5 strain) were transfected with either D-LDH- or BCAT-expressing vectors, and the over-expression of the corresponding transcripts was confirmed by RT-qPCR (data not shown). Afterwards, drug susceptibility of WT, WT + D-LDH, WT + BCAT, and WT + pLexsy (transfection control) lines to PMM was determined in both promastigote and amastigote forms of the different Leishmania species (Table 9). In all cases, both D-LDH- and BCAT- transfected promastigotes of L. donovani Dd8, L. major and L. infantum lines showed higher values of RI than the corresponding control lines (Table 9). Also, the intracellular amastigote forms showed higher RI values than the control lines, even though only the L. donovani Dd8 amastigotes over-expressing D-LDH or BCAT genes showed statistically significant differences in RI values in regard to the control ones (Table 9). Overall, these results suggest that D-LDH and BCAT genes can contribute to the resistance phenotype to PMM in other Leishmania species, although apparently to a lesser extent than in L. donovani HU3, probably due to their different genetic background.

Table 9

Paromomycin susceptibility of different species of Leishmania transfected with D-LDH and BCAT.

Lines	L. donovani Dd8EC50 μM (RI)a		L. majorEC50 μM (RI)		L. infantumEC50 μM (RI)
	Promastigotes	Intracellular amastigotes	Promastigotes	Intracellular amastigotes	Promastigotes	Intracellular amastigotes
WT	40.48 ± 4.09b	45.97 ± 5.97	6.60 ± 0.40	15.71 ± 1.52	40.16 ± 4.69	62.19 ± 7.88
WT+D-LDH	79.99 ± 7.20***(1.85)	63.61 ± 7.36*(1.53)	16.43 ± 2.01***(2.01)	18.98 ± 2.62(1.14)	74.87 ± 9.63***(1.88)	76.71 ± 9.62(1.20)
WT+BCAT	61.53 ± 2.82***(1.43)	65.02 ± 6.11**(1.56)	10.58 ± 0.97*(1.31)	21.02 ± 2.27(1.27)	56.17 ± 5.27*(1.42)	65.19 ± 7.12(1.02)
WT+Lexsy	43.03 ± 4.11	41.52 ± 5.15	8.10 ± 0.40	16.51 ± 1.94	39.68 ± 4.19	63.68 ± 5.23

Resistance indexes were calculated as described in Table 3 (for promastigotes) and Table 8 (for amastigotes).

Data are the mean ± SD of EC50 values from three independent experiments. Significant differences versus the controls were determined using Student‘s t-test (*p < 0.05; **p < 0.01; ***p < 0.005).

In addition, we analyzed the stability of resistance to PMM after removal of the drug pressure. After 2, 3, and 4 months of culturing the P line in absence of PMM, the EC50 values were determined. As shown in Fig. 7A, the EC50 values significantly decreased from the initial EC50 value of the P line (132.4 μM) to 89.2, 70.4, and 65.1 μM, after growing the P line in absence of PMM during 2, 3 and 4 months, respectively. In parallel, the levels of expression of D-LDH and BCAT genes were determined by RT-qPCR at the same times of drug removal (Fig. 7B). These data showed that the expression of D-LDH gene decreased to 54.9, 41.5, and 25.3% of the initial value, whereas the expression of BCAT gene decreased to 65.9, 43.6, and 51.7%, in the P line grown 2, 3 and 4 months without drug pressure (Fig. 7B). Therefore, we can conclude the existence of a direct relationship between PMM resistance and expression levels of these genes in L. donovani HU3, supporting further the involvement of D-LDH and BCAT genes in the resistance to PMM in Leishmania.

Fig. 7

Stability of the resistance to PMM and levels of expression of (A) The stability of resistance with respect to the initial EC50 of P line (month 0) was determined at 2, 3, and 4 months after removal of the drug pressure, using an MTT-based assay, as described in Materials and Methods. The results shown are the mean ± SD from three independent experiments. Asterisks denotes statistically significant differences (p < 0.05; Wilcoxon–Mann–Whitney test) between the P line before and after drug removal. (B) Relative expression levels of D-LDH and BCAT transcripts was determined by RT-qPCR before and after 2, 3 and 4 months after removal of drug pressure. Total RNA from different lines was extracted as described in Materials and Methods. Also, RT-qPCR expression values of the genes were normalized with the expression level of GAPDH transcript. The expression levels of the transcripts in the P line was set at 100%.The results shown are the mean ± SD from three independent experiments. Significant differences versus the control (P line) were determined using Student's t-test (*p < 0.01; **p < 0.005; ***p < 0.001).

It is known that PMM interferes with mitochondrial activity, disrupting the ΔΨm; even though it has been proposed that this drug would be acting at a metabolic level upstream of the respiratory chain (Maarouf et al., 1997; Jhingran et al., 2009). Our results showed that the PMM-resistant line analyzed in this study presents a low accumulation of Rh123 (Fig. 8), indicating that this line also presents a diminished ΔΨm, as described for other experimental PMM resistant parasites (Jhingran et al., 2009). Considering that D-LDH and BCAT would be involved in metabolic pathways to obtaining energy (Opperdoes and Coombs, 2007), it was determined whether the over-expression of either D-LDH or BCAT might affect the ΔΨm by interfering with metabolic processes leading to mitochondrial ATP synthesis. As depicted in Fig. 8, no depolarization was observed in promastigotes over-expressing D-LDH (WT + D-LDH line). However, promastigotes that over-expressed BCAT (WT + BCAT line), similar to the PMM-resistant promastigotes (P line), showed a significant decrease in the ΔΨm with respect to the control (WT + Lexsy line) (Fig. 8). The same behavior was observed when promastigotes were exposed to PBS for 2 h (Fig. 8), as an inducer of mitochondrial stress. As described, nutritional starvation induces autophagy in Leishmania (Besteiro et al., 2006); being autophagy required for the maintenance of mitochondrial homeostasis (Williams et al., 2013). In conclusion, this finding would suggest that the alteration of ΔΨm may be a consequence of a direct interference of the BCAT activity with oxidative phosphorylation in the Leishmania mitochondrion. Enzymes for the degradation of branched-chain amino acids are predicted to localize into the mitochondrion (Afanador et al., 2014). Thus, it is plausible that the aberrant overexpression of BCAT and/or the accumulation of a metabolite produced by this enzyme in the mitochondria could determine a block of the electron transport chain. In contrast, parasites overexpressing D-LDH presented ΔΨm values similar to controls (Fig. 8), suggesting that the activation of glycolytic and methylglyoxal pathways, in which this enzyme is involved (Greig et al., 2009), does not significantly influence the mitochondrial activity.

Fig. 8

Mitochondrial membrane potential in different The different L. donovani lines were left untreated (black columns) or exposed to PBS (white columns) for 2 h. Promastigotes were labeled for 15 min with 0.5 μM of Rh123 as described in Materials and Methods. Measurements by flow-cytometry are expressed as mean of arbitrary fluorescence units ± SD from three independent experiments. Asterisks denotes statistically significant differences (p < 0.05; Wilcoxon–Mann–Whitney test) between the control line (WT + Lexsy) and both the P line and the BCAT- transfected line (WT + BCAT).

In addition, since changes on the ΔΨm have been associated with ROS production, which in turn are associated with cell death, it was analyzed whether the overexpression of D-LDH and BCAT could be related to protection of parasites from PBS-induced ROS production. For this purpose, overall intracellular ROS levels were measured using the cell-permeable probe H2DCFDA in WT + Lexsy, WT + D-LDH, WT + BCAT, and P lines (Fig. 9A). The results showed that WT + Lexsy parasites (control) incubated with PBS for 2 h exhibited a marked increase in H2DCFDA fluorescence signal (79.7 versus 9.3 of untreated parasites; Fig. 9A), supporting an increased ROS production upon exposure to nutritional stress induced by PBS. In contrast, when the P line was exposed to PBS a 2.2-fold lower accumulation of the fluorescent probe than in control parasites was observed (Fig. 9A), indicating a lower ROS production and/or a higher antioxidant response in the PMM-resistant line regarding the control line, as expected from a P resistant line (Garcia-Hernandez et al., 2012). Remarkably, levels of ROS in the transfected D-LDH and BCAT Leishmania lines (WT + D-LDH, WT + BCAT, respectively) were very low, and they did not significantly increase when the over-expressing parasites were incubated in PBS (Fig. 9A). Thus, it can be concluded that Leishmania lines over-expressing D-LDH or BCAT generated (or accumulated) lower ROS amounts than control, and also than the P line (Fig. 9A), suggesting that these enzymes, or their metabolites, are involved in metabolic processes that increase the ability of parasites to neutralize the ROS production induced by nutritional stress.

Fig. 9

PBS-induced ROS generation in The L. donovani lines were left untreated (black columns) or exposed to PBS (white columns) for 2 h. The parasites were incubated with 4 μM H2DCFDA (A) or 5 μM MitoSOX (B) for 30 min at 28 °C. The fluorescence intensity was determined by flow-cytometry analysis and expressed as arbitrary fluorescence units (a.u). The data are the mean ± SD values from three independent experiments. Significant differences (p < 0.05; Wilcoxon–Mann–Whitney test) when comparing the different untreated lines with the untreated control line (WT + Lexsy) are denoted by asterisks, whereas significant differences between untreated and PBS-treated parasites for each line are marked by an ‘a’.

Additionally, considering that the major source of oxidants in eukaryotes is the mitochondria, we repeated the assays but using the MitoSOX red dye, as a probe for measuring mitochondrial ROS levels. As shown in panel B of Fig. 9, incubation of parasites in PBS also induced a marked increase in MitoSOX red fluorescence levels in the control line versus untreated parasites, but not in Leishmania lines over-expressing D-LDH or BCAT (Fig. 9B). Also, PMM-resistant parasites accumulated 1.7-fold less probe than the control line after incubation with PBS (Fig. 9B). The results indicated the ability of D-LDH and BCAT proteins to protect Leishmania parasites against oxidative stress also inside the mitochondria.

All these results suggest that suppression of mitochondrial activity by the action of PMM (Maarouf et al., 1997) leads to an activation of glycolytic activity as the most likely alternative for the production of cellular ATP (Chawla et al., 2011). It has been reported that Leishmania produce s limited amounts of D-lactate via the methylglyoxal pathway (Darling and Blum, 1988). This pathway functions as an overflow of glycolysis, in which triose-phosphates, formed in the glycolytic pathway, are converted to inorganic phosphate and methylglyoxal. Indeed the methylglyoxal pathway enzymes glyoxalase I (LinJ.35.3060) and II (LinJ.12.0200) are present in L. infantum. The resulting D-lactate would be converted to pyruvate by a D-LDH (LinJ.27.1940). The pyruvate resulting from both the glycolytic and the methylglyoxal pathways may be converted by the enzymes of the TCA cycle into building blocks for biosynthesis or acts as a scavenger of hydrogen peroxide, helping to reduce the oxidative stress. The compromised mitochondrion is likely to be crippled in its capacity for beta-oxidation of fatty acids and therefore must resort to an increase in metabolism of the branched amino acids, leucine, isoleucine and valine. The BCAT is the first enzyme of this oxidation pathway for all three amino acids. The resulting succinyl-CoA and acetyl-CoA is then converted into necessary building blocks for biosynthesis by the TCA cycle enzymes, or excreted as free acetate with the additional advantage of the formation of an extra molecule of ATP (Van Hellemond et al., 1998). Thus, under the pressure of PMM, overexpression of D-LDH and BCAT may allow for the production of additional energy in the form of ATP to compensate for mitochondrial dysfunction. Consequently, overexpression of D-LDH and BCAT could be considered as new molecular markers of drug resistance to PMM in Leishmania that will require future studies of validation in clinical and experimental resistant lines to PMM.

Apart from transcripts for these two enzymes, D-LDH and BCAT, among the differential expression transcripts, we found the transcript encoding for a member of the ABC family transporter (ABCA10; transcript LinJ.29.T0640) was present at higher levels (∼1.8-fold) in the P line (see Supplementary file 1). It is noticeable that this protein was found to be overexpressed in an L. donovani line resistant to difluoromethylornithine (Singh et al., 2014), and it may be postulated a role of this transporter in PMM exclusion. In fact, many members of the superfamily of ABC transporters have been found to be involved in drug resistance (see above (Sauvage et al., 2009)). Taking into account that ribosome would be another target for PMM (Fernandez et al., 2011; Zhang et al., 2016), another up-regulated transcript worthy to note was LinJ.03.T0240, coding for a ribosomal protein (L38). Finally, in the P line, similar to that occurring in the S and A lines, transcripts encoding for amastins (LinJ.29.T1450 and LinJ.29.T3000) were found to be upregulated. These results may be indicating that overexpression of amastin genes in Leishmania might lead to an increased parasite fitness against the cellular stresses elicited by drugs.

Conclusions

Our study illustrates the battery of molecular adaptations exploited by Leishmania to develop drug resistance, including gene over-expression, deletion, SNPs generating stop codons or amplification of sets of genes. This also further validates our approach (selection of drug resistance followed by NGS analysis of whole genome and transcriptome) as it allows the simultaneous identification of gain-of- and loss-of-function. Our study also highlights the power of the method, as it identifies in each line several genes that could be participating to drug resistance and accompanying effects (virulence for instance). Validation of the candidate driver genes is essential and was done here for the strongest candidates, for which largest differences in gene expression were observed. This allowed to identify new drivers, like SMT, the deletion of which being associated with resistance to AmB and the tandem D-LDH-BCAT, the amplification of which being related to PMM resistance. Our results provide highly relevant results for the understanding of drug resistance in Leishmania, but also for the design of molecular tools for its monitoring in clinical conditions.