Skip to main content

Comparative analysis of the transcriptional responses of five Leishmania species to trivalent antimony



Leishmaniasis is a neglected tropical disease caused by several species of Leishmania. The resistance phenotype of these parasites depends on the characteristics of each species, which contributes to increased therapeutic failures. Understanding the mechanism used by the parasite to survive under treatment pressure in order to identify potential common and specific therapeutic targets is essential for the control of leishmaniasis. The aim of this study was to investigate the expression profiles and potential shared and specific resistance markers of the main Leishmania species of medical importance [subgenus L. (Leishmania): L. donovani, L. infantum and L. amazonensis; subgenus L. (Viannia): L. panamensis and L. braziliensis)] resistant and sensitive to trivalent stibogluconate (SbIII).


We conducted comparative analysis of the transcriptomic profiles (only coding sequences) of lines with experimentally induced resistance to SbIII from biological replicates of five Leishmania species available in the databases of four articles based on ortholog attribution. Simultaneously, we carried out functional analysis of ontology and reconstruction of metabolic pathways of the resulting differentially expressed genes (DEGs).


Resistant lines for each species had differential responses in metabolic processes, compound binding, and membrane components concerning their sensitive counterpart. One hundred and thirty-nine metabolic pathways were found, with the three main pathways comprising cysteine and methionine metabolism, glycolysis, and the ribosome. Differentially expressed orthologous genes assigned to species-specific responses predominated, with 899 self-genes. No differentially expressed genes were found in common among the five species. Two common upregulated orthologous genes were found among four species (L. donovani, L. braziliensis, L. amazonensis, and L. panamensis) related to an RNA-binding protein and the NAD(P)H cytochrome-B5-oxidoreductase complex, associated with transcriptional control and de novo synthesis of linoleic acid, critical mechanisms in resistance to antimonials.


Herein, we identified potential species-specific genes related to resistance to SbIII. Therefore, we suggest that future studies consider a treatment scheme that is species-specific. Despite the limitations of our study, this is the first approach toward unraveling the pan-genus genetic mechanisms of resistance in leishmaniasis.

Graphical Abstract


Leishmania (Trypanosomatida: Trypanosomatidae) comprises more than 50 species (grouped into five subgenera), at least 20 of which are considered pathogenic for humans [1]. It is the etiological agent of leishmaniasis, a group of tropical infectious diseases transmitted by the bite of the sand fly [2]. It is estimated that 2 million new cases occur annually worldwide, with 70,000 deaths and 350 million people at risk of acquiring the disease [3]. Leishmaniasis, in turn, presents the particularity of a broad spectrum of clinical manifestations that may depend on the Leishmania species and/or the host's immune response [4, 5].

Leishmania species are classified into subgenera according to propagation within the vector's digestive tract. In Leishmania (Leishmania), the development occurs in the mid- and foregut, and in Leishmania (Viannia) in the hindgut [6]. Similarly, both subgenera present a disjunctive geographical distribution pattern: L. (Viannia) is restricted to the neotropics and L. (Leishmania) is found in both New World and Old World regions [7]. Its associated clinical manifestations are linked to the specific species of the subgenera. In L. (Leishmania), we find New World species in the L. mexicana complex, where L. amazonensis stands out, one of the dominant pathogenic species in Latin America, responsible for diffuse forms of cutaneous leishmaniasis (CL). Here, we also find Old World species, the L. donovani complex that causes the visceral form of the disease, with L. infantum and L. donovani, species with particular endemicity in southern Asia (India, Bangladesh, and Nepal) and north-east Africa (Sudan and Ethiopia) [8, 9]. In L. (Viannia), L. braziliensis and L. panamensis stand out, representing the most frequent causes of CL and mucosal leishmaniasis (MCL) in Latin America, especially in Colombia, Ecuador, and Brazil [10].

The number of chromosomes of each species also differs in L. infantum and L. donovani (36), L. amazonensis (34), and L. braziliensis and L. panamensis (35). This occurs due to the fusion of chromosomes 8 and 39 and 20 and 36 in the L. mexicana complex, compared to the L. major genome. In the L. braziliensis complex, the fusion occurs with chromosomes 20 and 34 [11]. Leishmania chromosomes are made up of sizable polycistronic transcription units (PTUs) of functionally unrelated genes. Although not tightly regulated, RNA synthesis has been associated with an accumulation of acetylated histone H3 at the end of the chromosome and in the strand-change regions, marking the start and end of PTU transcription. The presence of a novel modified nucleotide base (β-d-glucopyranosyloxymethyluracil), known as Base J, plays a critical role in ensuring that transcription termination occurs only at the end of each PTU, rather than at the polyadenylation sites of each gene [12,13,14]. Polycistronic pre-mRNAs are processed into mature trans-spliced mRNAs coupled by a leader RNA and polyadenylation as means of genetic regulation [15].

Around 200 genes are reported with differential distribution between L. major and L. infantum compared with L. braziliensis, where the majority encode for the parasite's survival in the macrophage in addition to a series of unique genes for each of the species [16]. A marker of the subgenus Leishmania (Viannia) is the RNA interference machinery, which functions as an antiviral defense mechanism against the Leishmania RNA virus-1 (LRV1), which seems to impact the virulence of the parasite [17, 18]. Another signature of this subgenus is high expression of molecules such as glycoprotein (GP) GP63 and GP46, known surface components that act on the tropism of CL, development inside the vector, and neutralization of the host macrophage defense system [19,20,21,22,23]. The subgenus Leishmania (Leishmania) has targets associated with visceral leishmaniasis (VL), such as the genes encoding proteins A2 and 6-phosphogluconate dehydrogenase (6PGDH), besides differential expression of genes such as the catalytic subunit of DNA polymerase A (POLA), heat-shock protein 20 (HSP20), or N-acetylglucosamine-phosphotransferase (NAGT) [24,25,26].

There is still no vaccine approved for human use; thus chemotherapy and various drugs are the only alternatives for treating the different clinical manifestations of leishmaniasis [27]. These include antibiotics with leishmanicidal action (pentamidine or paromomycin), antitumor treatments (miltefosine), antifungal drugs (imidazoles and triazoles), and essential leishmanicidal drugs (amphotericin B) [28]. Over the past 60 years, antimonials (Sb) (e.g., sodium stibogluconate [SSG] and meglumine antimoniate) have been used as the standard treatment in many countries, mainly due to poor socioeconomic conditions [29,30,31,32]. However, with antimonials, an extensive series of therapeutic failures have been reported, associated with multiple factors, including in the patients (adherence to treatment, geographical location, and immune status), the practitioners (failures in the antimonial administration protocols), and the characteristics of the parasite (drug resistance) [32,33,34,35,36,37,38,39].

Multiple studies have focused on determining the mechanisms used by this pathogen to survive under the pressure of antimonials using next-generation sequencing (NGS). Genomic studies have elucidated genetic variations in the number of copies of genes and specific mutations in resistant strains [40,41,42]. The transcriptomic profile (by RNA-seq) of species exposed to trivalent antimony (SbIII) has suggested that Leishmania uses differential pathways and global changes in gene expression to reach an adaptive phenotype, as reported in L. infantum [43], L. donovani [44], L. amazonensis [45], L. panamensis, and L. braziliensis [46]. These and other studies have demonstrated that this behavior generally goes hand in hand with an alteration in the number of copies of particular genes locally or throughout the chromosome [47,48,49]. Here, the transcripts associated with resistant lines involved a series of metabolic pathways directly associated with drugs and antibiotics or involved in glycolysis and fatty acid processing, in contrast to reducing the expression of factors associated with motility and replicative processes, related to the ribosomal machinery. In addition, some specific genes or processes have been reported as potential resistance targets, including transporters of molecules and metal compounds, heat-shock proteins such as HSP70, and genes encoding for amastins or transmembrane channels such as aquaporin 1 (AQP1) [50, 51]. Other interesting results include a higher number of transcripts associated with resistance to trivalent antimony than other treatments for L. donovani. Autophagy induction by L. amazonensis represents a survival or cell death strategy or a different expression of genes that encode for some ribosomal proteins whose function in the SSG-resistant Leishmania species is still unknown [44, 46, 52, 53].

Factors that may contribute to virulence and pathogenicity have been investigated in different species responsible for the presentations of the disease [54, 55]. Some of these studies also include lines resistant to SbIII and other drugs or focus on factors closely related to treatment, thus establishing the influence of conserved and species-specific traits [56,57,58]. Despite this, to date, no interspecific comparisons of transcriptomic profiles have been conducted among Leishmania species with induced resistance to SbIII. Therefore, we studied the expression profiles and potential shared and specific resistance markers of the main Leishmania species of medical importance [subgenus L. (Leishmania): L. donovani, L. infantum, and L. amazonensis; subgenus L. (Viannia): L. panamensis and L. braziliensis] resistant to SbIII.


Sequence and transcriptome data download

To homogenize the data (RNA-seq reads, stranded, paired) for each of the five Leishmania species (L. donovani, L. infantum, L. amazonensis, L. panamensis, and L. braziliensis), we used the limiting number of replicates of the articles with fewer samples (L. infantum and L. donovani), and the others were adapted. Then, two biological replicates corresponding to populations with experimentally induced resistance to SbIII and another two to their wild-type counterparts (sensitive to SbIII) were used per species (each one composed of two paired RNA-seq reads), available in the European Nucleotide Archive (ENA) databases and National Center for Biotechnology Information (NCBI) in fastq.gz format ( (Additional file 1: Table S1). These were downloaded using GNU Wget software that is part of the application console. Finally, the data were selected, taking as criteria the half maximal effective concentration (EC50), and using promastigotes in the log phase. The quality control of the crude sequences was determined using the FastQC tool, using all parameters except those involving GC content [59]. Subsequently, we removed the adapters previously found with the FastQC Adapter Content report using the CutAdapt tool (, and the low-quality reads were removed using Trimmomatic software (TOPHRED33; SLIDINGWINDOW:4:28; MINLEN:36) [60, 61]. Simultaneously, the reference genomes (.fasta) and the annotation files (.gff) of each species were downloaded from the kinetoplastid parasite database, EuPathDB ( ( [62].

Differential expression analysis

Initially, the annotated genomes were converted to Gene transfer format (.gtf) using the gffread v0.12.1 package for different gene-specific conventions such as chromosome location [63]. Then STAR v.2.7 software was used to index the reference genomes using the annotation files (.gtf). Finally, the previously mentioned fastq.gz files were aligned, generating a BAM file, again with STAR v.2.7 software (SortedByCoordinate; runThreadN:10; Per-sample 2-pass mapping) [64]. We then used the binary file and the indexed reference genome as input for HTSeq in the generation of a text file with a table containing the count for each gene. Using the intersection parameter of all non-empty sets and “order” options for the paired-end data, four were generated per species, two from wild-type reads and the other two from SbIII-resistant reads, counting only the coding sequences (CDS) [65].

The text file served as input for the R package, Bioconductor DESeq2 v.1.30.0, where the differential expression of the genes was evaluated from a binomial distribution model. Here, the wild-type samples served as a point of comparison for samples resistant to SbIII, generating a .csv file with the statistical information for the differentially expressed genes (DEGs) [66]. Afterwards, we filtered the sequences, leaving only those genes that encode known proteins, removing hypothetical proteins. This facilitates analyses such as determining orthologous groups and ontology, comparing profiles between species, and constructing tables. To visualize the DEGs, a volcano diagram was constructed showing the average normalized count under the following criteria of statistical significance: a cut-off point of the level of change > 1 or < 1, for the upregulated and downregulated genes, respectively, including a Benjamini–Hochberg test with an adjusted P-value < 0.01 [67].

Determining orthologous groups and ontology

We performed ontology enrichment analysis using the TriTrypDB tools ( with Fisher’s exact test with P-value < 0.05, using the DEG identification numbers (IDs) as an input. The gene ontology terms were classified into three categories (biological processes, cellular components, and molecular functions). These terms were refined in REVIGO (, where redundant terms were cut through the filter of essential genes. Subsequently, a bar diagram of the gene ontology terms with the highest gene count was made using the ggplot2 visualization package [68].

We assigned an orthologous group to each DEG to homogenize the data for later comparison among species using the additional parameters in the search for gene ID in TriTrypDB ( A bar diagram was constructed using the ggplot2 visualization package for the orthologous genes to illustrate the relationships with the Leishmania species. Those that were upregulated and downregulated were included based on the level of change. A heat map of the DEGs for each chromosome distributed among the five species was generated with the R ggplot2 package [68].

Establishing metabolic pathways

Using a list of DEG IDs (two per state, upregulated and downregulated for each species), FASTA files were generated through the EupathDB TriTryp tool. Each FASTA file (10 total) was uploaded to the Kyoto Encyclopedia of Genes and Genomes (KEGG) Automatic Annotation Server (KAAS) tool, where BLAST was used as the search algorithm and the kinetoplastid data set available in the manually curated KEGG database. Subsequent signaling pathway reconstruction was analyzed under the option "Pathway map" allowing a graphical visualization based on data grouping in the different metabolic pathways. We used Adobe Photoshop CC software to condense the most relevant data obtained in each graph, including the genes of the different species. Then, the pathways where they were involved were chosen (> 11 DEGs) [64]. Additionally, a table was generated that included the metabolic pathways where DEGs were involved and their location and description (Additional file 5: Table S5).


DEGs among Leishmania species resistant to SbIII

No significant differences between the biological replicates of each species were found with their limited sampling. A total of 1222 differentially expressed genes were identified, counting all species (P-value < 0.01), of which 711 were upregulated and 511 downregulated. Leishmania amazonensis presented the highest number of DEGs in both categories (248 and 225, respectively), while L. panamensis and L. infantum presented the lowest number of upregulated (33) and downregulated (16) DEGs. In addition, L. donovani had 172 and 73 and L. braziliensis 125 and 80 DEGs. Genes can be visualized in volcanoes, showing the relationship between mean expression and level of change for each one (Fig. 1; Additional file 2: Table S2; Additional file 3: Table S3). The distribution of the DEGs on the chromosomes of the species belonging to the subgenus L. (Leishmania) shows a homogeneous pattern for L. infantum, a significant number of DEGs on chromosome 29 of L. donovani (70), and a considerable number of genes on chromosome 34 for L. amazonensis (64). The New World species (L. amazonensis, L. braziliensis, and L. panamensis) also coincide with a high number of DEGs on chromosome 23 (Fig. 2; Additional file 2: Table S2, Additional file 3: Table S3).

Fig. 1

Differential expression between SbIII-resistant and SbIII-sensitive lines of different Leishmania species. Each point represents one gene. Gray dots indicate the genes that were not differentially expressed, and the blue dots, located above and below the black lines (cut-off for the fold change [> 1 and < − 1]), represent differentially expressed genes, adjusted P-value < 0.01

Fig. 2

DEG distribution on the chromosomes of Leishmania species. Darker colors represent a larger number of genes

Gene ontology (GO) of DEGs among Leishmania species resistant to SbIII

Different numbers of DEGs associated with ontological processes were found, as follows, with the first value indicating downregulated processes and the second upregulated: L. donovani (137, 128), L. infantum (46, 64), L. amazonensis (160, 69), L. braziliensis (62, 97), and L. panamensis (87, 35) (Additional file 4: Table S4). In addition, the prominent roles of the DEGs of the SbIII -resistant populations in each category of GO functional groups were determined (Fig. 3).

Fig. 3

Gene ontology of DEGs among Leishmania species resistant to SbIII. Green bars represent those in upregulated functional groups and red bars represent those in downregulated groups

A significant number of upregulated genes were observed in the L. (Leishmania) subgenus related to biological processes. In contrast, the pattern was reversed with the L. (Viannia) subgenus species, where the downregulated genes were related to different metabolic processes. Antimonial induction affects the metabolism of Leishmania (organic substances, drugs, and antibiotics) and cellular processes associated with redox, protein folding, and replication. The group of cellular components was characterized by an increase in the expression of the genes that encode membrane proteins and their integral components and their movement structure. On the other hand, the downregulated DEGs were related to internal components like cytoplasm and ribosome. The most representative upregulated genes with molecular function were associated with catalytic and transport processes and binding to small molecules, metal compounds, and drugs.

The downregulated genes followed a structural pattern, and we found a reduction in protein and nucleic acid binding, hydrolase, and oxidoreductase activity (Additional file 4: Table S4). At the individual level, the gene ontology with the highest enrichment was L. donovani (regulation of DNA replication [6899]), L. infantum (protein-binding complex [33,193]), L. amazonensis (metabolism of isocitrate [2288]), L. braziliensis (metabolism of cysteine and glycerol [5267]), and L. panamensis (metabolism of sphingolipids [8566]). This also highlights the few upregulated processes associated with the cellular component in these species (Additional file 4: Table S4).

Comparison of the transcriptomic profiles of Leishmania species

The orthologous genes of each species were predominant in proportion to the shared ones. This occurred for both downregulated (red bars) (361) and upregulated (green bars) (491) DEGs. They were followed by the shared genes between any two species (downregulated: 84; upregulated: 141), where each combination between two species had at least one shared ortholog. In L. donovani and L. infantum, involved in VL, three orthologs were found (OG6_100045, OG6_101028, OG6_101741), related to the transport of substances and catalytic events (Fig. 4, Additional file 2: Table S2).

Fig. 4

Number of orthologous genes shared and exclusive to each Leishmania species (full purple points). Green bars represent those that are upregulated and red bars represent those that are downregulated

Shared orthologs were found in three combinations of species. The first, among the subgenus L. (Viannia) species together with L. amazonensis, are distributed in the New World and are generally associated with CL (OG6_156932, OG6_100025, OG6_142640), related to microtubule motor proteins and molecular binding enzymes. Another combination of similarities in orthologs was between L. amazonensis, L. donovani, and L. panamensis (OG6_110995, OG6_127587, OG6_101746, and OG6_157042). Only one common ortholog was found in L. braziliensis, L. donovani, and L. infantum (OG6_119258). Finally, no common orthologs of the subgenus L. (Leishmania) were found (Fig. 4, Additional file 3: Table S3). No orthologous genes were shared among all species. Only two (OG6_152462 and OG6_200283) were shared among four species (L. donovani, L. braziliensis, L. amazonensis, and L. panamensis), both upregulated, associated with an RNA-binding protein and cytochrome b5, respectively. A bar plot illustrates the number of shared and self-DEGs of lines resistant to SbIII of Leishmania species (Fig. 4). A list of the specific and shared orthologs and the statistical information is also provided (Additional file 2: Table S2; Additional file 3: Table S3).

Metabolic pathways associated with DEGs of Leishmania species resistant to SbIII

A total of 139 metabolic pathways were found; 14 had more than 11 DEGs involved, and some were associated with the same component of the pathway. Of these, the three most representative were glycolysis/gluconeogenesis (ko00010) (26), cysteine and methionine (ko00270) (21), and ribosome (ko03010) (54) (Fig. 5a–c respectively). In these cases, differential responses were found between species, even in the same gene (e.g., acetyl-CoA synthetase in the glycolysis/gluconeogenesis pathway was upregulated in L. braziliensis but downregulated in L. panamensis). On the other hand, referring directly to resistance to treatments and pathology, we found signaling pathways of the metabolism of drugs and other enzymes (ko00983), resistance to antifolates (ko01523), resistance to platinum drugs (ko01524), and leishmaniasis (ko05140) (Additional file 5: Table S5).

Fig. 5

Most representative metabolic pathways for the differentially expressed genes. Glycolysis/gluconeogenesis (ko00010) (a), cysteine and methionine metabolism (ko00270) (b), and ribosome (ko03010) (c)


Here, we identified a large number (1222) of DEGs between sensitive and resistant lines in the studied Leishmania species (Fig. 1). This is possibly a consequence of the basic properties of the parasite, starting with the repertoire of unique genes and the expression patterns that impact its virulence, which do not necessarily have to be exclusively linked to resistance to treatments. Furthermore, Leishmania is exposed to other stressful conditions (e.g., oxidative stress in macrophages or contamination by arsenic in the environment), so the resistance mechanisms could be repeated in both situations [69,70,71]. Both repertoires of DEGs could be discriminated in future studies using cell lines involved during the infection process of the parasite without the treatment intervention. It is worth mentioning that herein we employed sequence data from different studies, and these could be subject to heterogeneous experimental conditions with regard to their conservation and manipulation in methodological processes not defined in the articles, as well as subsequent filtering of the raw sequences, a concerning limitation that we want to highlight. Additionally, we suggest the use of amastigotes in the future when treating the parasite's infective stage and underexposure from the treatment in a real case [72].

In species lacking transcriptional regulation such as Leishmania, one of the adaptive mechanisms to modulate the gene dose is the change in the number of chromosome copies, whose base number is a chromosomal mosaic where monosomies (like chromosome 2 in L. major) to supernumerary chromosomes (like chromosome 5 in L. major) can occur [73, 74]. Nevertheless, we cannot be sure that the change in chromosomal somy is directly responsible for the resistance phenotype [75]. In addition, there may be changes associated with more specific processes, such as the variation in the number of copies of regions or even variants generated at the nucleotide level [46, 76]. Leishmania can amplify and remove some loci regions through rearrangement and recombination events as a strategy against environmental changes [77, 78]. These could be some of the reasons behind the overrepresentation of DEGs at chromosome 23 in New World species and the relationship to CL (L. amazonensis, L. braziliensis, and L. panamensis) (Fig. 2).

A large proportion of the genes correspond to ABC transporters involved in the entry and distribution of the drug in the parasite cell [79, 80], which agrees with previous findings in L. infantum and L. major. This points to a potential biomarker of resistance beyond the pathogen's clinical manifestation and geographical distribution [81, 82], but future mechanistic studies are needed to confirm this hypothesis. Chromosome 34 of L. amazonensis and L. braziliensis has genes related to arginine metabolism. This pathway contributes to the maintenance of the reduction–oxidation balance and the expression of virulence factors such as GP63 and other genes that encode for amastins, the latter of which is vital in the differentiation of Leishmania [83], suggesting the need to understand the effect of these genes in the resistance phenotype.

We found GO terms associated mainly with metabolic processes, probably due to the inhibitory action on glycolysis or oxidation of fatty acids by trivalent antimony. This is consistent with previous studies in L. donovani and L. major (Fig. 3; Additional file 4: Table S4) [84,85,86]. Regarding molecular functions, we observed the gene ontology associated with binding to different compounds, especially those with ionic properties, probably related to metal-thiol conjugate sequestration by ATP binding cassettes (ABC transporters), a resistance mechanism widely reported in L. infantum and L. major [87,88,89]. For the structural components of the pathogen, the most significant was the membrane, likely due to properties such as the number of ergosterols required for protection against reactive oxygen species (ROS) generated in the biological activation of antimony and/or the presence of transporters [90, 91]. In future studies, each of these specific processes could be explored, seeking to inhibit or reduce the impact of treatment actions that could directly and/or indirectly interfere, setting up their role in the parasite's virulence and survival rate, including studies in vivo if possible.

As some product descriptions were redundant even after REVIGO scrubbing, we performed a signaling pathway reconstruction as a supplemental functional analysis. Three pathways stood out for the number of DEGs involved. The first was glycolysis/gluconeogenesis, the main basis of Leishmania promastigote metabolism. The latter is also considered a virulence factor, as it facilitates replication in the phagolysosome of macrophages and the generation of typical lesions in mice [92, 93] (Fig. 5a; Additional file 5: Table S5). Glucose transporters and fructose-1,6-bisphosphate, as gluconeogenic enzymes, have been validated as pharmacological flanks; in turn, this pathway indirectly helps to synthesize other flanks, such as mannan from hexoses [94]. The cysteine and methionine pathway has been found to be modulated in L. donovani when exposed to amphotericin B and oxidative stress in general. It also contributes to apoptosis inhibition, one of the mechanisms of action of miltefosine [95, 96] (Fig. 5b; Additional file 5: Table S5). Finally, studies related to the ribosome have found differential expression of the proteins of the large subunit (e.g., L23a) involved in translation and with a role in kinase activity in Leishmania resistant to pentavalent antimonials, miltefosine, and amphotericin [97, 98]. Likewise, ribosomal proteins play a role in cell growth and apoptosis control [99] (Fig. 5c; Additional file 5: Table S5). Also, in trypanosomatids, regulation of proteins in a specific subcellular compartment, such as ribosomes, is usually exponential in stages of parasite development such as early metacyclogenesis (eg., L5, L11, Rpf2, and Rrs1), showing its importance in translation repression, interaction with specific essential proteins, and its susceptibility to conditions such as oxidative stress [100,101,102,103].

Signaling pathway reconstruction showed a more detailed picture of the processes associated with GO terms, where agreement with reports of resistance and participation from biological processes, molecular functions, and cellular components are maintained. Some molecular targets involved in previously studied resistance mechanisms such as multidrug resistance-related protein A (MRPA) or key enzymes in trypanothione synthesis (e.g., gamma-glutamylcysteine synthase and spermidine) were differentially regulated (Additional file 3: Table S3; Additional file 5: Table S5). MRPA contributes to drug transport in the sequestration of SbIII–trypanothione complexes in compartments close to the flagellar pocket of the parasite [88]. On the other hand, trypanothione acts as a nonenzymatic antiparasitic reducing agent, counteracting the oxidative effects of antimonial components, and even preventing their action when conjugated with them [104]. However, we did not find differential expression of aquaglyceroporin (AQP1), a fundamental element in resistance to antimony mediated by drug entry [105]. Further proteomic studies might be conducted to investigate the impact between RNA-seq data and protein expression.

When comparing orthologous groups, we found two shared orthologs between species of both subgenera (L. amazonensis, L. braziliensis, L. donovani, and L. panamensis) (Fig. 5, Additional file 3: Table S3). OG6_152462, encoding for RNA-binding proteins (RBPs), was present in the post-transcriptional control processes of gene expression, something fundamental in trypanosomatids. It has been described as a therapeutic target of interest, as it affects the parasite's virulence because of altered expression [106, 107]. The other ortholog was OG6_200283, encoding for a component of the complex NAD(P)H-cytochrome-B5 oxidoreductase, whose upregulation can favor the survival of the parasite in a medium with lipid deficiency, something that it experiences not only when exposed to trivalent antimony but also during its life cycle. It is common in several species, as it participates in the de novo synthesis of linolenic acid, the most common unsaturated fatty acid in Leishmania, related to a decrease in oxidative stress ROS [108].

Our results show a more significant proportion of orthologs unique to each species than those shared between two or more species. At the subgenus level, there is no evidence of several significant orthologous relationships: within L. (Leishmania), there are no proper orthologs, and in L. (Viannia), only seven good orthologous genes are found (Fig. 4). This is probably because the mechanism of antimony resistance is complex and multifactorial, where the characteristics of each species can contribute to the final response, which reflects the taxonomic classification, the geographical location, and the clinical manifestation. However, a population-level analysis would provide further information [109,110,111]. On the other hand, some of these orthologs may be associated with a general stress response, and others with the resistance mechanisms of Leishmania to antimony, where resistant strains could modify the expected behavior of the metabolic pathways to elaborate their phenotype [112]. This goes together with the different mechanisms of action that the drug must have, including the fragmentation of genetic material, oxidative stress induction, interference in protein interactions, or an increase in intracellular Ca+ [113,114,115]. Mechanistic studies are needed to respond to these questions and fully understand the concept of resistance that seems species-specific.

Possible limitations of our methodology include the exclusion of non-coding regions (UTR) in the count of gene reads, which provide heterogeneity in transcription by regulating the RNA stability and/or the efficiency of translation [116, 117]. Limitations in functional annotation are inherent in one or several species. Also, there is remarkable specificity of orthologous groups, where two or more genes belong to the same metabolic pathway and/or molecular/structural process [118]. For example, the orthologs OG6_101492 and OG6_101499 both putatively encode for the structural maintenance of the chromosome in L. amazonensis, but one of them belongs to a protein family (Additional file 2: Table S2). Using only the sequence of two biological replicates for each condition can also have repercussions in generating potential deficiencies in sensitivity and statistical specificity or the absence of data on other differentially expressed genes; thus, additional validations are required. Another scenario that may contribute to a higher proportion of self-genes per species is that the ortholog is shared, but its differential expression is different, as occurs, for example, with OG6_119258, where this nucleobase transporter is downregulated in L. braziliensis, L. donovani, and L. infantum, but upregulated in L. panamensis (Additional file 2: Table S2; Additional file 3: Table S3). These changes in the transcriptome could also be correlated to alternative trans-splicing, a mechanism of heterogeneity described in trypanosomatids such as L. major and T. brucei [119, 120]. Also, changes may occur at the transcript level due to differences in mRNA maturation and stability mediated mainly by polyadenylation or deadenylation events, mainly by RNA-binding proteins [121].

The use of promastigotes instead of amastigotes must also be considered a limitation of the current study. Transcriptomic differences have been demonstrated in both stages, and amastigotes, known as the infective and intracellular stage of the parasite, represent the best alternative in in vitro drug exposure studies [72]. Experimental conditions that are likely to affect the results include the stability of the sensitive and resistant lines determined by their respective media. It is worth mentioning that the resistance selection rounds and their duration also varied depending on the species, but this went hand in hand with the necessary doses until reaching a common cut-off point (EC50). It should also be noted that resistance and other virulence factors can vary even intraspecifically. Thus, future studies should consider more than just coding sequences and should consider a more significant number of strains of each of the species evaluated in this work, preferably collected from common geographical areas [122, 123]. Establishing relationships between our results and those of studies carried out at the proteome and chromosomal structure level (general and local) could contribute to a better understanding of resistance mechanisms across species.

In addition, future studies should consider the role of hypothetical proteins and genes with unknown functions, as they represent a significant percentage of potential individual markers. Despite the limitations highlighted herein, this study is the first to provide a comprehensive comparative analysis between several species of different subgenera, providing a functional perspective of general (ontology) and specific (reconstruction of signaling pathways) forms, contributing to the basic knowledge necessary for the development of potential treatments based on epidemiology and the behavior of this parasite resistance linked to a public health problem.


Our study contributes to a better understanding of the trivalent antimony resistance patterns of various Leishmania species. We found that there was a species-specific response predominance, mainly associated with biological processes and cellular components. Likewise, some probable potential DEGs were evidenced that could be evaluated as therapeutic targets directed to a specific subgenus or clinical manifestation, which based on their function could be regulated in stressful conditions beyond exposure to SbIII, such as lipid deficiency and post-transcriptional control (RNA and cytochrome B5 binding protein). Furthermore, we highlight some limitations of the NGS approach in determining patterns of resistance mechanisms among different Leishmania species, and therefore the proposal of treatment strategies (species-specific or a holistic response), as well as other variables that should be considered for future research in this field.

Availability of data and materials

The data sets used and/or analyzed during the current study are available in Additional file 1: Table S1.



Differentially expressed genes


Trivalent antimony


Cutaneous, visceral, mucosal leishmaniasis


Polycistronic transcription units


Leishmania RNA virus-1

GP63, GP46:

Glycoprotein 63, 46

HSP70, HSP20:

Heat-shock protein 70, 20


6-Phosphogluconate dehydrogenase


DNA polymerase A


N-acetylglucosamine phosphotransferase


Aquaporin 1

EC50 :

Half maximal effective mean concentration


Reactive oxygen species


Coding sequence


  1. 1.

    Akhoundi M, Kuhls K, Cannet A, Votýpka J, Marty P, Delaunay P, et al. A historical overview of the classification, evolution, and dispersion of Leishmania parasites and sandflies. PLoS Negl Trop Dis. 2016;10: e0004349.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  2. 2.

    Alemayehu B, Alemayehu M. Leishmaniasis: a review on parasite, vector and reservoir host. Heal Sci J. 2017;11:519.

    Article  Google Scholar 

  3. 3.

    Oryan A, Akbari M. Worldwide risk factors in leishmaniasis. Asian Pac J Trop Med. 2016;9:925–32.

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Hashiguchi Y, Gomez EA. Importance of Leishmania species and vector sand fly (Diptera: Psychodidae) identification. J Med Entomol. 2018;55:773–4.

    Article  PubMed  Google Scholar 

  5. 5.

    Galluzzi L, Ceccarelli M, Diotallevi A, Menotta M, Magnani M. Real-time PCR applications for diagnosis of leishmaniasis. Parasit Vectors. 2018;11:1–13.

    CAS  Article  Google Scholar 

  6. 6.

    Dostálová A, Volf P. Leishmania development in sand flies: Parasite-vector interactions overview. Parasit Vectors. 2012;5:276.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  7. 7.

    Steverding D. The history of leishmaniasis. Parasit Vectors. 2017;10:82.

    Article  PubMed  PubMed Central  Google Scholar 

  8. 8.

    Burza S, Croft SL, Boelaert M. Leishmaniasis. Lancet. 2018;392:951–70.

    Article  PubMed  Google Scholar 

  9. 9.

    Arenas R, Torres-Guerrero E, Quintanilla-Cedillo MR, Ruiz-Esmenjaud J. Leishmaniasis: a review. F1000Res. 2017;6:750.

  10. 10.

    Valero NNH, Uriarte M. Environmental and socioeconomic risk factors associated with visceral and cutaneous leishmaniasis: a systematic review. Parasitol Res. 2020;119:365–84.

    Article  PubMed  Google Scholar 

  11. 11.

    Britto C, Ravel C, Bastien P, Blaenau C, Pagès M, Dedet JP, et al. Conserved linkage groups associated with large-scale chromosomal rearrangements between Old World and New World Leishmania genomes. Gene. 1998;222:107.

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Thomas S, Green A, Sturm NR, Campbell DA, Myler PJ. Histone acetylations mark origins of polycistronic transcription in Leishmania major. BMC Genomics. 2009;10:152.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Grünebast J, Clos J. Leishmania: responding to environmental signals and challenges without regulated transcription. Computat Struct Biotechnol J. 2020;18:4016–23.

    CAS  Article  Google Scholar 

  14. 14.

    Jensen BC, Phan IQ, McDonald JR, Sur A, Gillespie MA, Ranish JA, et al. Chromatin-associated protein complexes link DNA base J and transcription termination in leishmania. mSphere. 2021;6:e01204-e1220.

    CAS  Article  PubMed  Google Scholar 

  15. 15.

    Clayton C. Regulation of gene expression in trypanosomatids: living with polycistronic transcription. Open Biol. 2019;9:1900072.

    CAS  Article  Google Scholar 

  16. 16.

    Peacock CS, Seeger K, Harris D, Murphy L, Ruiz JC, Quail MA, et al. Comparative genomic analysis of three Leishmania species that cause diverse human disease. Nat Genet. 2007;39:839–47.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  17. 17.

    Brettmann EA, Shaik JS, Zangger H, Lye LF, Kuhlmann FM, Akopyants NS, et al. Tilting the balance between RNA interference and replication eradicates Leishmania RNA virus 1 and mitigates the inflammatory response. Proc Natl Acad Sci USA. 2016;113:11998–2005.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  18. 18.

    Kariyawasam R, Mukkala AN, Lau R, Valencia BM, Llanos-Cuentas A, Boggild AK. Virulence factor RNA transcript expression in the Leishmania Viannia subgenus: influence of species, isolate source, and Leishmania RNA virus-1. Trop Med Health. 2019;47:25.

    Article  PubMed  PubMed Central  Google Scholar 

  19. 19.

    Llanes A, Restrepo CM, Del VG, Anguizola FJ, Lleonart R. The genome of Leishmania panamensis: insights into genomics of the L. (Viannia) subgenus. Sci Rep. 2015;5:8550.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  20. 20.

    Diotallevi A, Buffi G, Ceccarelli M, Neitzke-Abreu HC, Gnutzmann LV, da Costa Lima MS, et al. Real-time PCR to differentiate among Leishmania (Viannia) subgenus, Leishmania (Leishmania) infantum and Leishmania (Leishmania) amazonensis: Application on Brazilian clinical samples. Acta Trop. 2020;201: 105178.

    CAS  Article  PubMed  Google Scholar 

  21. 21.

    Soulat D, Bogdan C. Function of macrophage and parasite phosphatases in leishmaniasis. Front Immunol. 2017;8:1838.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  22. 22.

    Ovalle-Bracho C, Camargo C, Díaz-Toro Y, Parra-Muñoz M. Molecular typing of Leishmania (Leishmania) amazonensis and species of the subgenus Viannia associated with cutaneous and mucosal leishmaniasis in Colombia: a concordance study. Biomedica. 2018;38:86–95.

    Article  PubMed  Google Scholar 

  23. 23.

    Marín M, Aguilar YA, Ramírez JR, Triana O, Muskus CE. Molecular and immunological analyses suggest the absence of hydrophilic surface proteins in Leishmania (Viannia) panamensis. Biomedica. 2008;28:423–32.

    Article  PubMed  Google Scholar 

  24. 24.

    Akhoundi M, Downing T, Votýpka J, Kuhls K, Lukeš J, Cannet A, et al. Leishmania infections: molecular targets and diagnosis. Mol Aspects Med. 2017;57:1–29.

    CAS  Article  PubMed  Google Scholar 

  25. 25.

    Fraga J, Montalvo AM, Van der Auwera G, Maes I, Dujardin JC, Requena JM. Evolution and species discrimination according to the Leishmania heat-shock protein 20 gene. Infect Genet Evol. 2013;18:229–37.

    CAS  Article  PubMed  Google Scholar 

  26. 26.

    Fernandes AP, Canavaci AMC, McCall LI, Matlashewski G. A2 and other visceralizing proteins of Leishmania: Role in pathogenesis and application for vaccine development. Subcell Biochem. 2014;74:77–101.

    CAS  Article  PubMed  Google Scholar 

  27. 27.

    Jain K, Jain NK. Vaccines for visceral leishmaniasis: A review. J Immunol Methods. 2015;422:1–12.

    CAS  Article  PubMed  Google Scholar 

  28. 28.

    Lindoso JAL, Costa JML, Queiroz IT, Goto H. Review of the current treatments for leishmaniases. Res Rep Trop Med. 2012;3:69–77.

    Article  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Eddaikra N, Ait-Oudhia K, Kherrachi I, Oury B, Multi-Mati F, Benikhlef R, et al. Antimony susceptibility of Leishmania isolates collected over a 30-year period in Algeria. PLoS Negl Trop Dis. 2018;12: e0006310.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  30. 30.

    Balasegaram M, Ritmeijer K, Lima MA, Burza S, Ortiz Genovese G, Milani B, et al. Liposomal amphotericin B as a treatment for human leishmaniasis. Expert Opin Emerg Drugs. 2012;17:493–510.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Burza S, Sinha PK, Mahajan R, Lima MA, Mitra G, Verma N, et al. Five-year field results and long-term effectiveness of 20 mg/kg liposomal amphotericin B (Ambisome) for visceral leishmaniasis in Bihar. India PLoS Negl Trop Dis. 2014;8: e2603.

    CAS  Article  PubMed  Google Scholar 

  32. 32.

    Vanlerberghe V, Diap G, Guerin PJ, Meheus F, Gerstl S, Van Der SP, et al. Drug policy for visceral leishmaniasis: a cost-effectiveness analysis. Trop Med Int Heal. 2007;12:274–83.

    CAS  Article  Google Scholar 

  33. 33.

    Sundar S, Chakravarty J. An update on pharmacotherapy for leishmaniasis. Expert Opin Pharmacother. 2015;16:237–52.

    CAS  Article  PubMed  Google Scholar 

  34. 34.

    Sundar S, Chakravarty J, Meena LP. Leishmaniasis: treatment, drug resistance and emerging therapies. Expert Opin Orphan Drugs. 2019;7:1–10.

    CAS  Article  Google Scholar 

  35. 35.

    Rojas R, Valderrama L, Valderrama M, Varona MX, Ouellette M, Saravia NG. Resistance to antimony and treatment failure in human Leishmania (Viannia) infection. J Infect Dis. 2006;193:1375–83.

    CAS  Article  PubMed  Google Scholar 

  36. 36.

    Denis S, Carla M, Khatima AO. Antimony resistance and environment: elusive links to explore during Leishmania life cycle. Int J Parasitol Drugs Drug Resist. 2012;2:200–3.

    CAS  Article  Google Scholar 

  37. 37.

    Romero GA, Vinitius De Farias Guerra M, Gomes Paes M, de Oliveira Macêdo V. Comparison of cutaneous leishmaniasis due to Leishmania (Viannia) braziliensis and L. (V.) guyanensis in Brazil: clinical findings and diagnostic approach. Clin Infect Dis. 2001;32:1304–12.

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Arevalo J, Ramirez L, Adaui V, Zimic M, Tulliano G, Miranda-Verástegui C, et al. Influence of Leishmania (Viannia) species on the response to antimonial treatment in patients with American tegumentary leishmaniasis. J Infect Dis. 2007;195:1846–51.

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    de Vries HJC, Reedijk SH, Schallig HDFH. Cutaneous leishmaniasis: recent developments in diagnosis and management. Am J Clin Dermatol. 2015;16:99–109.

    Article  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Borsari C, Jiménez-Antón MD, Eick J, Bifold E, Torrado JJ, Olías-Molero AI, et al. Discovery of a benzothiophene-flavonol halting miltefosine and antimonial drug resistance in Leishmania parasites through the application of medicinal chemistry, screening and genomics. Eur J Med Chem. 2019;183: 111676.

    CAS  Article  PubMed  Google Scholar 

  41. 41.

    Zheng Z, Chen J, Ma G, Satoskar AR, Li J. Integrative genomic, proteomic and phenotypic studies of Leishmania donovani strains revealed genetic features associated with virulence and antimony-resistance. Parasit Vectors. 2020;13:510.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  42. 42.

    Dumetz F, Cuypers B, Imamura H, Zander D, D’Haenens E, Maes I, et al. Molecular preadaptation to antimony resistance in Leishmania donovani on the Indian Subcontinent. mSphere. 2018;3:e000548-17.

    Article  Google Scholar 

  43. 43.

    Andrade JM, Gonçalves LO, Liarte DB, Lima DA, Guimarães FG, de Melo RD, et al. Comparative transcriptomic analysis of antimony resistant and susceptible Leishmania infantum lines. Parasit Vectors. 2020;13:600.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  44. 44.

    Rastrojo A, García-Hernández R, Vargas P, Camacho E, Corvo L, Imamura H, et al. Genomic and transcriptomic alterations in Leishmania donovani lines experimentally resistant to antileishmanial drugs. Int J Parasitol Drugs Drug Resist. 2018;8:246–64.

    Article  PubMed  PubMed Central  Google Scholar 

  45. 45.

    Patino LH, Muskus C, Ramírez JD. Transcriptional responses of Leishmania (Leishmania) amazonensis in the presence of trivalent sodium stibogluconate. Parasit Vectors. 2019;12:348.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  46. 46.

    Patino LH, Imamura H, Cruz-Saavedra L, Pavia P, Muskus C, Méndez C, et al. Major changes in chromosomal somy, gene expression and gene dosage driven by SbIII in Leishmania braziliensis and Leishmania panamensis. Sci Rep. 2019;9:9485.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  47. 47.

    Iantorno SA, Durrant C, Khan A, Sanders MJ, Beverley SM, Warren WC, et al. Gene expression in Leishmania is regulated predominantly by gene dosage. MBio. 2017;8:e01393-e1417.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  48. 48.

    Restrepo CM, Llanes A, Cedeño EM, Chang JH, Álvarez J, Ríos M, et al. Environmental conditions may shape the patterns of genomic variations in Leishmania panamensis. Genes (Basel). 2019;10:838.

    CAS  Article  Google Scholar 

  49. 49.

    Leprohon P, Légaré D, Raymond F, Madore É, Hardiman G, Corbeil J, et al. Gene expression modulation is associated with gene amplification, supernumerary chromosomes and chromosome loss in antimony-resistant Leishmania infantum. Nucleic Acids Res. 2009;37:1387–99.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  50. 50.

    Monte-Neto R, Laffitte MCN, Leprohon P, Reis P, Frézard F, Ouellette M. Intrachromosomal amplification, locus deletion and point mutation in the Aquaglyceroporin AQP1 gene in antimony resistant Leishmania (Viannia) guyanensis. PLoS Negl Trop Dis. 2015;9: e0003476.

    Article  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Maharjan M, Madhubala R. Heat shock protein 70 (HSP70) expression in antimony susceptible/resistant clinical isolates of Leishmania donovani. Nepal J Biotechnol. 2015;3:22–8.

    Article  Google Scholar 

  52. 52.

    Ponte-Sucre A, Gamarro F, Dujardin JC, Barrett MP, López-Vélez R, García-Hernández R, et al. Drug resistance and treatment failure in leishmaniasis: a 21st century challenge. PLoS Negl Trop Dis. 2017;11: e0006052.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  53. 53.

    Verma A, Bhandari V, Deep DK, Sundar S, Dujardin JC, Singh R, et al. Transcriptome profiling identifies genes/pathways associated with experimental resistance to paromomycin in Leishmania donovani. Int J Parasitol Drugs Drug Resist. 2017;3:370–7.

    Article  Google Scholar 

  54. 54.

    Depledge DP, Evans KJ, Ivens AC, Aziz N, Maroof A, Kaye PM, et al. Comparative expression profiling of Leishmania: modulation in gene expression between species and in different host genetic backgrounds. PLoS Negl Trop Dis. 2009;3: e476.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Rochette A, Raymond F, Ubeda JM, Smith M, Messier N, Boisvert S, et al. Genome-wide gene expression profiling analysis of Leishmania major and Leishmania infantum developmental stages reveals substantial differences between the two species. BMC Genomics. 2008;9:255.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  56. 56.

    Hefnawy A, Berg M, Dujardin JC, De Mulder G. Exploiting knowledge on Leishmania drug resistance to support the quest for new drugs. Trends Parasitol. 2017;33:162–74.

    CAS  Article  PubMed  Google Scholar 

  57. 57.

    Andrade JM, Murta SMF. Functional analysis of cytosolic tryparedoxin peroxidase in antimony-resistant and -susceptible Leishmania braziliensis and Leishmania infantum lines. Parasit Vectors. 2014;7:406.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  58. 58.

    Matrangolo FSV, Liarte DB, Andrade LC, De Melo MF, Andrade JM, Ferreira RF, et al. Comparative proteomic analysis of antimony-resistant and-susceptible Leishmania braziliensis and Leishmania infantum chagasi lines. Mol Biochem Parasitol. 2013;190:63–75.

    CAS  Article  PubMed  Google Scholar 

  59. 59.

    Andrews S. FASTQC a quality control tool for high throughput sequence data. Babraham Inst. 2015.

  60. 60.

    Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17:10–2.

  61. 61.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  62. 62.

    Aurrecoechea C, Barreto A, Basenko EY, Brestelli J, Brunk BP, Cade S, et al. EuPathDB: the eukaryotic pathogen genomics database resource. Nucleic Acids Res. 2017;45:D581–91.

    CAS  Article  PubMed  Google Scholar 

  63. 63.

    Pertea M, Pertea G. GFF Utilities: GffRead and GffCompare. F1000Res. 2020;9:ISCB Comm J-304.

  64. 64.

    Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.

    CAS  Article  PubMed  Google Scholar 

  65. 65.

    Anders S, Pyl PT, Huber W. HTSeq-A Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.

    CAS  Article  PubMed  Google Scholar 

  66. 66.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  67. 67.

    Batut B, Hiltemann S, Bagnacani A, Baker D, Bhardwaj V, Blank C, et al. Community-driven data analysis training for biology. Cell Syst. 2018;6:752–8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  68. 68.

    Wickham H. ggplot2: elegant graphics for data analysis. Springer-Verlag New York. 2016. ISBN 978-3-319-24277-4.

  69. 69.

    Uliana SRB, Trinconi CT, Coelho AC. Chemotherapy of leishmaniasis: present challenges. Parasitology. 2018;145:464–80.

    CAS  Article  PubMed  Google Scholar 

  70. 70.

    Decuypere S, Vanaerschot M, Brunker K, Imamura H, Müller S, Khanal B, et al. Molecular mechanisms of drug resistance in natural Leishmania populations vary with genetic background. PLoS Negl Trop Dis. 2012;6: e1514.

    Article  PubMed  PubMed Central  Google Scholar 

  71. 71.

    Yardley V, Ortuño N, Llanos-Cuentas A, Chappuis F, De Doncker S, Ramirez L, et al. American tegumentary leishmaniasis: is antimonial treatment outcome related to parasite drug susceptibility? J Infect Dis. 2006;194:1168–75.

    Article  PubMed  Google Scholar 

  72. 72.

    Vermeersch M, da Luz RI, Toté K, Timmermans JP, Cos P, Maes L. In vitro susceptibilities of Leishmania donovani promastigote and amastigote stages to antileishmanial reference drugs: practical relevance of stage-specific differences. Antimicrob Agents Chemother. 2009;53:3855–9.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  73. 73.

    Rogers MB, Hilley JD, Dickens NJ, Wilkes J, Bates PA, Depledge DP, et al. chromosome and gene copy number variation allow major structural change between species and strains of Leishmania. Genome Res. 2011;21:2129–42.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  74. 74.

    Sterkers Y, Lachaud L, Crobu L, Bastien P, Pagès M. FISH analysis reveals aneuploidy and continual generation of chromosomal mosaicism in Leishmania major. Cell Microbiol. 2011;13:274–83.

    CAS  Article  PubMed  Google Scholar 

  75. 75.

    Downing T, Stark O, Vanaerschot M, Imamura H, Sanders M, Decuypere S, et al. Genome-wide SNP and microsatellite variation illuminate population-level epidemiology in the Leishmania donovani species complex. Infect Genet Evol. 2012;12:149–59.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Papadopoulou B, Ouellette M, Laffitte MCN, Leprohon P. Plasticity of the Leishmania genome leading to gene copy number variations and drug resistance. F1000Res 2016;5:2350.

  77. 77.

    Ubeda JM, Raymond F, Mukherjee A, Plourde M, Gingras H, Roy G, et al. Genome-wide stochastic adaptive DNA amplification at direct and inverted DNA repeats in the parasite Leishmania. PLoS Biol. 2014;12: e1001868.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  78. 78.

    Papadopoulou B, Ouellette M, Laffitte MCN, Leprohon P. Plasticity of the Leishmania genome leading to gene copy number variations and drug resistance. F1000Res. 2016;5:2350.

  79. 79.

    Douanne N, Wagner V, Roy G, Leprohon P, Ouellette M, Fernandez-Prada C. MRPA-independent mechanisms of antimony resistance in Leishmania infantum. Int J Parasitol Drugs Drug Resist. 2020;13:28–37.

    Article  PubMed  PubMed Central  Google Scholar 

  80. 80.

    Barrera MC, Rojas LJ, Weiss A, Fernandez O, McMahon-Pratt D, Saravia NG, et al. Profiling gene expression of antimony response genes in Leishmania (Viannia) panamensis and infected macrophages and its relationship with drug susceptibility. Acta Trop. 2017;176:355–63.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  81. 81.

    Brotherton MC, Bourassa S, Leprohon P, Légaré D, Poirier GG, Droit A, et al. Proteomic and genomic analyses of antimony resistant Leishmania infantum mutant. PLoS ONE. 2013;8: e81899.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  82. 82.

    Mukherjee A, Boisvert S, Monte-Neto RL do, Coelho AC, Raymond F, Mukhopadhyay R, et al. Telomeric gene deletion and intrachromosomal amplification in antimony-resistant Leishmania. Mol Microbiol. 2013;88:189–202.

  83. 83.

    Acuña SM, Aoki JI, Laranjeira-Silva MF, Zampieri RA, Fernandes JCR, Muxel SM, et al. Arginase expression modulates nitric oxide production in Leishmania (Leishmania) amazonensis. PLoS One. 2017;12:e0187186. https://

  84. 84.

    Singh N, Sundar S. Integrating genomics and proteomics permits identification of immunodominant antigens associated with drug resistance in human visceral leishmaniasis in India. Exp Parasitol. 2017;176:30–45.

    CAS  Article  PubMed  Google Scholar 

  85. 85.

    Rabhi I, Rabhi S, Ben-Othman R, Rasche A, Daskalaki A, Trentin B, et al. Transcriptomic signature of Leishmania infected mice macrophages: a metabolic point of view. PLoS Negl Trop Dis. 2012;6: e1763.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  86. 86.

    Ghosh AK, Sardar AH, Mandal A, Saini S, Abhishek K, Kumar A, et al. Metabolic reconfiguration of the central glucose metabolism: a crucial strategy of Leishmania donovani for its survival during oxidative stress. FASEB J. 2015;29:2081–98.

    CAS  Article  PubMed  Google Scholar 

  87. 87.

    Manzano JI, García-Hernández R, Castanys S, Gamarro F. A new ABC half-transporter in Leishmania major is involved in resistance to antimony. Antimicrob Agents Chemother. 2013;57:3719–30.

  88. 88.

    El Fadili K, Messier N, Leprohon P, Roy G, Guimond C, Trudel N, et al. Role of the ABC transporter MRPA (PGPA) in antimony resistance in Leishmania infantum axenic and intracellular amastigotes. Antimicrob Agents Chemother. 2005;49:1988–93.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  89. 89.

    Légaré D, Richard D, Mukhopadhyay R, Stierhof YD, Rosen BP, Hammer A, et al. The Leishmania ATP-binding Cassette protein PGPA is an intracellular metal-thiol transporter ATPase. J Biol Chem. 2001;276:26301–7.

    Article  PubMed  Google Scholar 

  90. 90.

    Mathur R, Das RP, Ranjan A, Shaha C. Elevated ergosterol protects Leishmania parasites against antimony-generated stress. FASEB J. 2015;29:4201–13.

    CAS  Article  PubMed  Google Scholar 

  91. 91.

    Frézard F, Monte-Neto R, Reis PG. Antimony transport mechanisms in resistant Leishmania parasites. Biophys Rev. 2014;6:119–32.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  92. 92.

    Naderer T, Ellis MA, Serene MF, De Souza DP, Curtis J, Handman E, et al. Virulence of Leishmania major in macrophages and mice requires the gluconeogenic enzyme fructose-1,6-bisphosphatase. Proc Natl Acad Sci U S A. 2006;103:5502–7.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  93. 93.

    Biyani N, Singh AK, Mandal S, Chawla B, Madhubala R. Differential expression of proteins in antimony-susceptible and -resistant isolates of Leishmania donovani. Mol Biochem Parasitol. 2011;179:91–9.

    CAS  Article  PubMed  Google Scholar 

  94. 94.

    Opperdoes FR, Michels P a M. The metabolic repertoire of Leishmania and implications for drug discovery. Leishmania. Caister Academic Press. 2008.

  95. 95.

    Singh K, Ali V, Pratap Singh K, Gupta P, Suman SS, Ghosh AK, et al. Deciphering the interplay between cysteine synthase and thiol cascade proteins in modulating Amphotericin B resistance and survival of Leishmania donovani under oxidative stress. Redox Biol. 2017;12:350–66.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  96. 96.

    Kumar R, Tiwari K, Dubey VK. Methionine aminopeptidase 2 is a key regulator of apoptotic like cell death in Leishmania donovani. Sci Rep. 2017;7:95.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  97. 97.

    Shalev-Benami M, Zhang Y, Rozenberg H, Nobe Y, Taoka M, Matzov D, et al. Atomic resolution snapshot of Leishmania ribosome inhibition by the aminoglycoside paromomycin. Nat Commun. 2017;8:1589.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  98. 98.

    Das S, Shah P, Baharia RK, Tandon R, Khare P, Sundar S, et al. Over-expression of 60s ribosomal L23a is associated with cellular proliferation in SAG resistant clinical isolates of Leishmania donovani. PLoS Negl Trop Dis. 2013;7: e2527.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  99. 99.

    Ranjan R, Das P, Vijayakumar S. Differentially modulated proteins associated with leishmaniasis—a systematic review of in-vivo and in-vitro studies. Mol Biol Rep. 2020;47:9159–78.

    CAS  Article  PubMed  Google Scholar 

  100. 100.

    Avila CC, Mule SN, Rosa-Fernandes L, Viner R, Barisón MJ, Costa-Martins AG, et al. Proteome-wide analysis of Trypanosoma cruzi exponential and stationary growth phases reveals a subcellular compartment-specific regulation. Genes (Basel). 2018;9:413.

    CAS  Article  Google Scholar 

  101. 101.

    Jaremko D, Ciganda M, Christen L, Williams N. Trypanosoma brucei L11 is essential to ribosome biogenesis and interacts with the kinetoplastid-specific proteins P34 and P37. mSphere. 2019;4:e00475-e519.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  102. 102.

    Jaremko D, Ciganda M, Williams N. Trypanosoma brucei homologue of regulator of ribosome synthesis 1 (Rrs1) has direct interactions with essential trypanosome-specific proteins. mSphere. 2019;4:e00453-e519.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  103. 103.

    Romaniuk MA, Frasch AC, Cassola A. Translational repression by an RNA-binding protein promotes differentiation to infective forms in Trypanosoma cruzi. PLoS Pathog. 2018;14: e1007059.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  104. 104.

    Mohebali M, Kazemirad E, Hajjaran H, Kazemirad E, Oshaghi MA, Raoofian R, et al. Gene expression analysis of antimony resistance in Leishmania tropica using quantitative real-time PCR focused on genes involved in trypanothione metabolism and drug transport. Arch Dermatol Res. 2019;311:9–17.

    CAS  Article  PubMed  Google Scholar 

  105. 105.

    Potvin J-E, Leprohon P, Queffeulou M, Sundar S, Ouellette M. Mutations in an aquaglyceroporin as a proven marker of antimony clinical resistance in the parasite Leishmania donovani. Clin Infect Dis. 2021;72:e526–32.

    CAS  Article  PubMed  Google Scholar 

  106. 106.

    Rashidi S, Kalantar K, Fernandez-Rubio C, Anvari E, Nguewa P, Hatam G. Chitin binding protein as a possible RNA binding protein in Leishmania parasites. Pathog Dis. 2020;78:ftaa007.

  107. 107.

    Mukherjee S, Santara S Sen, Das S, Bose M, Roy J, Adak S. NAD(P)H Cytochrome b5 oxidoreductase deficiency in Leishmania major results in impaired linoleate synthesis followed by increased oxidative stress and cell death. J Biol Chem. 2012;287:34992–35003.

  108. 108.

    Mukherjee A, Adhikari A, Das P, Biswas S, Mukherjee S, Adak S. Loss of virulence in NAD(P)H cytochrome b5 oxidoreductase deficient Leishmania major. Biochem Biophys Res Commun. 2018;503:371–7.

    CAS  Article  PubMed  Google Scholar 

  109. 109.

    Fernández OL, Diaz-Toro Y, Ovalle C, Valderrama L, Muvdi S, Rodríguez I, et al. Miltefosine and antimonial drug susceptibility of Leishmania Viannia Species and populations in regions of high transmission in colombia. PLoS Negl Trop Dis. 2014;8: e2871.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  110. 110.

    Jeddi F, Mary C, Aoun K, Harrat Z, Bouratbine A, Faraut F, et al. heterogeneity of molecular resistance patterns in antimony-resistant field isolates of Leishmania species from the western Mediterranean area. Antimicrob Agents Chemother. 2014;58:4866–74.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  111. 111.

    Bañuls AL, Hide M, Prugnolle F. Leishmania and the leishmaniases: a parasite genetic update and advances in taxonomy, epidemiology and pathogenicity in humans. Adv Parasitol. 2007;64:1–109.

    Article  PubMed  Google Scholar 

  112. 112.

    Chakravarty J, Sundar S. Drug resistance in leishmaniasis. J Glob Infect Dis. 2010;2:167–76.

    Article  PubMed  PubMed Central  Google Scholar 

  113. 113.

    Haldar AK, Sen P, Roy S. Use of antimony in the treatment of leishmaniasis: current status and future directions. Mol Biol Int. 2011;2011: 571242.

    Article  PubMed  PubMed Central  Google Scholar 

  114. 114.

    Croft SL, Sundar S, Fairlamb AH. Drug resistance in leishmaniasis. Clin Microbiol Rev. 2006;19:111–26.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  115. 115.

    Aslett M, Aurrecoechea C, Berriman M, Brestelli J, Brunk BP, Carrington M, et al. TriTrypDB: a functional genomic resource for the Trypanosomatidae. Nucleic Acids Res. 2009;38:D457–62.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  116. 116.

    Quijada L, Soto M, Alonso C, Requena JM. Identification of a putative regulatory element in the 3’-untranslated region that controls expression of HSP70 in Leishmania infantum. Mol Biochem Parasitol. 2000;110:79–91.

    CAS  Article  PubMed  Google Scholar 

  117. 117.

    Rastrojo A, Carrasco-Ramiro F, Martín D, et al. The transcriptome of Leishmania major in the axenic promastigote stage: transcript annotation and relative expression levels by RNA-seq. BMC Genomics. 2013;14:223.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  118. 118.

    Dillon LAL, Okrah K, Hughitt KV, Suresh R, Li Y, Fernandes MC, et al. Transcriptomic profiling of gene expression and RNA processing during Leishmania major differentiation. Nucleic Acids Res. 2015;43:6799–813.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  119. 119.

    Michaeli S. Trans-splicing in trypanosomes: Machinery and its impact on the parasite transcriptome. Future Microbiol. 2011;6:459–74.

    CAS  Article  PubMed  Google Scholar 

  120. 120.

    Clayton CE. Gene expression in Kinetoplastids. Curr Opin Microbiol. 2016;32:46–51.

    CAS  Article  PubMed  Google Scholar 

  121. 121.

    Urrea DA, Duitama J, Imamura H, Álzate JF, Gil J, Muñoz N, et al. Genomic analysis of Colombian Leishmania panamensis strains with different level of virulence. Sci Rep. 2018;8:17336.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  122. 122.

    Rugani JN, Quaresma PF, Gontijo CF, Soares RP, Monte-Neto RL. Intraspecies susceptibility of Leishmania (Viannia) braziliensis to antileishmanial drugs: antimony resistance in human isolates from atypical lesions. Biomed Pharmacother. 2018;108:1170–80.

    CAS  Article  PubMed  Google Scholar 

  123. 123.

    Downing T, Imamura H, Decuypere S, Clark TG, Coombs GH, Cotton JA, et al. Whole-genome sequencing of multiple Leishmania donovani clinical isolates provides insights into population structure and mechanisms of drug resistance. Genome Res. 2011;21:2143–56.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

Download references


We thank the authors who generated the data and sequences used in this article.


This work was funded by Dirección de Investigación e Innovación from Universidad del Rosario.

Author information




JM, LC, LHP, MM analyzed the data and performed bioinformatics analyses. JM wrote the manuscript. JDR led the project and wrote the manuscript. All authors have read and approved the final version of the manuscript.

Corresponding author

Correspondence to Juan David Ramírez .

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

General information about the samples and reference genomes used in this study. Samples were obtained from the European Nucleotide Archive (ENA) and reference genomes from TriTrypDB.

Additional file 2: Table S2.

Shared orthologous shared between Leishmania species.

Additional file 3: Table S3.

Exclusive orthologs of Leishmania species.

Additional file 4: Table S4.

Gene ontology (three functional groups: cellular component, molecular function, and biological process) of Leishmania species.

Additional file 5: Table S5.

Metabolic pathways and differentially expressed genes involved among Leishmania species.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Medina, J., Cruz-Saavedra, L., Patiño, L.H. et al. Comparative analysis of the transcriptional responses of five Leishmania species to trivalent antimony. Parasites Vectors 14, 419 (2021).

Download citation


  • Leishmania
  • Transcriptomic profile
  • Trivalent antimony
  • Orthologous groups
  • Pathway reconstruction