# Comparative transcriptomic analysis of antimony resistant and susceptible Leishmania infantum lines

Juvana Moreira Andrade, Leilane Oliveira Gonçalves, Daniel Barbosa Liarte, Davi Alvarenga Lima, Frederico Gonçalves Guimarães, Daniela de Melo Resende, Ana Maria Murta Santi, Luciana Marcia de Oliveira, João Paulo Linhares Velloso, Renato Guimarães Delfino, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Juvana Moreira Andrade, Leilane Oliveira Gonçalves, Daniel Barbosa Liarte, Davi Alvarenga Lima, Frederico Gonçalves Guimarães, et al. Comparative transcriptomic analysis of antimony resistant and susceptible Leishmania infantum lines. Parasites & Vectors, 2020, 13 (1), pp.600. 10.1186/s13071-020-04486-4. pasteur-03109036

HAL Id: pasteur-03109036

https://pasteur.hal.science/pasteur-03109036

Submitted on 13 Jan 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
RESEARCH

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

# Comparative transcriptomic analysis of antimony resistant and susceptible *Leishmania infantum* lines

Juvana Moreira Andrade<fnref n="1" />†, Leilane Oliveira Gonçalves<fnref n="2" />,³†, Daniel Barbosa Liarte⁴, Davi Alvarenga Lima<fnref n="1" />,<fnref n="2" />, Frederico Gonçalves Guimarães<fnref n="2" />, Daniela de Melo Resende<fnref n="1" />,<fnref n="2" />, Ana Maria Murta Santi<fnref n="1" />, Luciana Marcia de Oliveira⁵, João Paulo Linhares Velloso<fnref n="2" />, Renato Guimarães Delfino<fnref n="2" />, Pascale Pescher⁶, Gerald F. Späth⁶, Jeronimo Conceição Ruiz<fnref n="2" />,³*† and Silvane Maria Fonseca Murta<fnref n="1" />*† ID

**Abstract**

**Background:** One of the major challenges to leishmaniasis treatment is the emergence of parasites resistant to antimony. To study differentially expressed genes associated with drug resistance, we performed a comparative transcriptomic analysis between wild-type and potassium antimony tartrate (Sb³⁺)-resistant *Leishmania infantum* lines using high-throughput RNA sequencing.

**Methods:** All the cDNA libraries were constructed from promastigote forms of each line, sequenced and analyzed using STAR for mapping the reads against the reference genome (*L. infantum* JPCM5) and DESeq2 for differential expression statistical analyses. All the genes were functionally annotated using sequence similarity search.

**Results:** The analytical pipeline considering an adjusted *p*-value < 0.05 and fold change > 2.0 identified 933 transcripts differentially expressed (DE) between wild-type and Sb³⁺-resistant *L. infantum* lines. Out of 933 DE transcripts, 504 presented functional annotation and 429 were assigned as hypothetical proteins. A total of 837 transcripts were upregulated and 96 were downregulated in the Sb³⁺-resistant *L. infantum* line. Using this DE dataset, the proteins were further grouped in functional classes according to the gene ontology database. The functional enrichment analysis for biological processes showed that the upregulated transcripts in the Sb³⁺-resistant line are associated with protein phosphorylation, microtubule-based movement, ubiquitination, host–parasite interaction, cellular process and other categories. The downregulated transcripts in the Sb³⁺-resistant line are assigned in the GO categories: ribonucleoprotein complex, ribosome biogenesis, rRNA processing, nucleosome assembly and translation.

**Conclusions:** The transcriptomic profile of *L. infantum* showed a robust set of genes from different metabolic pathways associated with the antimony resistance phenotype in this parasite. Our results address the complex and multifactorial antimony resistance mechanisms in *Leishmania*, identifying several candidate genes that may be further evaluated as molecular targets for chemotherapy of leishmaniasis.

<footnote>*Correspondence: jeronimo.ruiz@fiocruz.br; silvane.murta@fiocruz.br</footnote>

<footnote>†Juvana Moreira Andrade, Leilane Oliveira Gonçalves, Jeronimo Conceição Ruiz and Silvane Maria Fonseca Murta contributed equally to this work</footnote>

<footnote n="1">¹Genômica Funcional de Parasitos, Instituto René Rachou, Fiocruz Minas, Belo Horizonte, MG, Brazil</footnote>

<footnote n="2">²Grupo Informática de Biossistemas, Instituto René Rachou, Fiocruz Minas, Belo Horizonte, MG, Brazil</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
**Keywords:** *Leishmania infantum*, Trivalent antimony, Resistance, RNA sequencing, Transcriptome, Differentially expressed genes

**Background**

Leishmaniasis is a complex of diseases caused by different species of the protozoan parasite *Leishmania* (Kinetoplastida, Trypanosomatidae) that represents one of the major public health problems in developing countries according to the World Health Organization [1]. Human leishmaniasis has a prevalence of 12 million cases and an incidence of 0.7–1.0 million new cases annually from nearly 100 endemic countries, with an estimated population of more than 1 billion people at risk of infection [1, 2]. Depending on genetic and environmental factors, the host immune response and mainly on the *Leishmania* species involved, the disease can comprise two main clinical forms: visceral or cutaneous leishmaniasis [3]. Visceral leishmaniasis (VL)—caused by *Leishmania (Leishmania) donovani* in Asia and Africa and *Leishmania (Leishmania) infantum* (syn. *L. (L. chagasi*) in the Mediterranean Basin, the Middle East, Central Asia, South America and Central America—is the most severe, systemic form and is lethal if not treated [3, 4]. The estimated incidence of VL is approximately 50,000 to 90,000 cases per year, and this disease remains endemic in more than 60 countries [1]. More than 95% of global VL cases occur in ten countries: Brazil, China, Ethiopia, India, Iraq, Kenya, Nepal, Somalia, South Sudan and Sudan [1].

There is no human vaccine available against *Leishmania* infections, and the control is based mainly on chemotherapy. Pentavalent antimony-containing compounds such as sodium stibogluconate (Pentostam®) and N-methyl-glucamine (Glucantime®) have been used as the first-line therapies against all forms of leishmaniasis. Although antimony's action has not been fully elucidated, studies suggest that pentavalent antimony (SbV) is reduced in vivo to the trivalent active form (SbIII) [5]. Literature data have indicated that antimony inhibits macromolecule biosynthesis in amastigotes, possibly via inhibition of glycolysis and fatty acid oxidation [6]. An earlier report indicated that antimonials cause perturbations in the thiol redox potential, driving to parasite death by oxidative stress [7]. Other studies have shown that antimony causes DNA fragmentation and can kill the parasite by an apoptotic process [8, 9]. In addition, zinc finger proteins have also been recognized as potential targets of SbIII [10].

One major challenge for leishmaniasis treatment is the emergence of parasites resistant to SbV [11, 12]. Treatment failure with SbV has been reported in Bihar (India), where more than 60% of patients with VL are

unresponsive to SbV [4, 13]. Different mechanisms of drug resistance have been reported [11], including decreased SbIII entry into the cell due to reduced expression of aquaglyceroporin (AQP1) [14–16] or sequestration of the metal–thiol conjugate into vesicular membranes of *Leishmania* by the ATP-binding cassette transporter [17, 18] and greater SbIII efflux due to amplification of ABC transporters [19].

Decuypere et al. [20] showed that the molecular changes associated with antimonial resistance in *Leishmania* isolates depend on their genetic background. To understand the mechanisms responsible for drug resistance in *Leishmania*, different approaches have been used. These include gene expression analyses of antimony-resistant *L. amazonensis* by DNA microarrays [21], proteomic analyses of SbIII-resistant *L. braziliensis* and *L. infantum* [22] and phosphoproteomic analysis of SbIII-resistant and -susceptible *L. braziliensis* [23].

Several studies showed the use of next-generation sequencing technologies to contribute to a better understanding of *Leishmania* biology. These have been widely used for comparing the gene expression profiles of primary cutaneous lesions from patients infected with *L. braziliensis* [24], to analyze the global changes in gene expression during *L. major* differentiation from procyclic to metacyclic forms [25] and for an overview of the global transcriptome of the *L. major* promastigote stage [26]. Recently, transcriptomic changes in an in vitro-adapted *L. amazonensis* in response to SbIII and comparative genomic and transcriptome analysis of SbIII-resistant and -susceptible *L. braziliensis* and *L. panamensis* were performed using DNA and RNA sequencing [27, 28]. To the best of our knowledge, the transcriptomic analysis associated with antimony resistance in *L. infantum* has not yet been addressed. Thus, this study attempts to perform a comparative transcriptomic analysis (RNA-seq) between SbIII-resistant and wild-type *L. infantum* lines.

**Methods**

*Leishmania* samples

This study used lines of *L. infantum* (MHOM/BR/74/PP75) wild type (LiWTS) and resistant to potassium antimonyl tartrate (SbIII) (LiSbR). The resistant line (LiSbR) was previously selected in vitro by a step-wise increase of SbIII drug pressure [29]. These parasites were further maintained in culture under SbIII pressure,
and the effective concentration required to decrease growth by 50% ($EC_{50}$) was determined using a model Z1 Coulter Counter (Beckman Coulter, Fullerton, CA, USA). $EC_{50}$ values for LiWTS and LiSbR obtained in this study were 0.12 mg/ml and 1 mg/ml, respectively, with an eight-fold resistance index (data not shown). Then, promastigote forms of these lines were grown at 26 °C in M199 medium supplemented with 40 mM HEPES pH 7.4, 1 µg/ml biotin, 5 µg/ml hemin, 2 µg/ml biopterin, 2 mM L-glutamine, 500 U penicillin, 50 µg/ml streptomycin and 10% (v/v) heat-inactivated fetal bovine serum [29]. Three independent biological replicates of each line were cultured. Based on previous studies of our group [29], wild-type *L. infantum* parasites were incubated for 24 h in the absence of drug (LiWTS 0), and resistant parasites were treated with 0.06 mg/ml $Sb^{III}$ ($LiSbR$ 0.06), which corresponds to half of the $Sb^{III}$ $EC_{50}$ for the LiWTS line. Cells were washed in RPMI medium, pelleted by centrifugation and frozen at -70 °C.

### RNA-Seq library preparation and sequencing

Promastigote forms were harvested, lysed and homogenized in the presence of guanidine-thiocyanate-containing buffer, and total RNA was extracted using the RNA extraction kit (RNeasy-QIAGEN, Valencia, CA, USA), according to the manufacturer's instructions. After extraction, total RNA was analyzed on the Agilent Bioanalyzer (Santa Clara, CA, USA) for quality and integrity assessment and, after this, submitted to cDNA synthesis. All samples presented an RNA integrity number ($RIN$) ≥ 6.8.

The construction of six non-directional libraries was prepared using the TruSeq RNA Sample preparation v2 protocol (Illumina, Inc., San Diego, CA), using 5 µg of total RNA for each library. The Illumina HiSeq2000 technology (Illumina, Inc., San Diego, CA) of Sequencing Platform of the Institut Pasteur was used to sequence the samples, based on directional sequencing of 100-bp-long reads of retro-transcribed mRNAs.

### Genome data used

*Leishmania infantum* JPCM5 genome data were downloaded from *European Nucleotide Archive* (ENA; http://www.ebi.ac.uk/ena/) under accession number ENSTUDY-CBMSO-04-04-2017 [30]. This genome version refers to the resequencing of the *L. infantum* JPCM5 genome at the end of 2017. According to Fuente et al. [30], 495 new genes have been annotated, 100 have been corrected, and 75 previous annotated genes have been discontinued. The TriTrypDB version contains a chromosome LinJ.00 formed by 34 genomic regions of uncertain chromosomal location, and in the new genome

version all regions were identified in the correct chromosomal location.

### Data quality control

Raw sequence reads in FASTQ format were evaluated in terms of read quality (per base sequence quality, per base G+C content, sequence length distribution, sequence duplication levels, Kmer content and low complexity sequences) with PRINSEQ [31]. Data filtering and trimming were performed with Trimmomatic [32]. Sequence artifacts such as sequencing adapters were removed using data available in the Trimmomatic software package.

One cDNA library from the $LiSbR$ line sequenced was removed from RNA seq analysis, since it showed a low throughput and small coverage (< 60×) compared to the other two libraries from the same *L. infantum*-resistant line (approximately 200×).

A curated General Feature Format (GFF) file was generated from the updated genome annotation and used to guide the alignment process. Reads were aligned in the reference genome with STAR [33], allowing up to three mismatches per read.

Mapped reads were converted to SAM format with SAMtools [34] and visualized with the Integrative Genome Viewer (IGV) [35].

### Differential expression analysis

To perform the differential gene expression analysis, HTSeq-count [36] was used to count the total number of mapped reads for each annotated gene in the GFF file. For differential gene expression discovery, DESeq2 [37] was used. To identify differentially expressed (DE) genes, an adjusted *p*-value < 0.05 and fold-change (FC) > 2.0 were set as thresholds to define the significance.

### Functional analysis

Blast2GO [38, 39] version 5.1.13 was used to map the DE genes in the gene ontology (GO) database [40]. The functional enrichment analysis (Fisher's exact test) was performed in the Blast2GO software using as test set the lists of DE genes and as reference the *L. infantum* JPCM5 predicted proteome. A false discovery rate (FDR) and an adjusted *p*-value < 0.05 were set as thresholds to define the functional enrichment significance.

## Results

### Overview of samples sequencing

The purpose of this study was to compare the transcriptome of $Sb^{III}$-resistant ($LiSbR$) and wild-type (LiWTS) *L. infantum* lines. For this, cDNA libraries from these samples were constructed, sequenced and analyzed, allowing the identification of differential gene expression associated with $Sb^{III}$ resistance mechanisms.
Three independent biological replicates of the wild-
type parasites and two of the Sb<sup>III</sup>-resistant parasites were
sequenced, producing ~500 million 100 base pair reads.
After quality trimming (adaptor removal and Phred qual-
ity cutoff ≥ 25), approximately 2% of the reads were lost.
After mapping, approximately 388 million reads were
aligned to a reference genome (*L. infantum* JPCM5).

Differential expression analysis

In the dataset comprising 8591 protein coding tran-
scripts (obtained from ENA database), 933 (933/8591,
10.86%) DE transcripts were identified between wild-
type and Sb<sup>III</sup>-resistant *L. infantum* lines considering the
applied cutoffs of adjusted *p*-value < 0.05 and FC > 2.0.
Out of 933 DE transcripts, 504 (504/933, 54.01%) pre-
sented functional annotation and 429 (429/933, 45.99%)
were assigned as hypothetical proteins without predicted
function (Table 1). A total of 837 (837/933, 89.71%)
transcripts were upregulated and 96 (96/933, 10.29%)
were downregulated in the Sb<sup>III</sup>-resistant *L. infantum*
(Table 1).

The functional enrichment analysis of 933 transcripts obtained in this study was performed on Blast2GO software [38, 39]. Out of 933 transcripts, 644 (644/933, 69.02%) were associated with some GO ontology related to the biological process (BP), molecular function (MF) or cellular component (CC) and 289 (298/933, 30.98%) did not present any associated GO term (Table 1). A total of 207 (207/644, 32.14%) DE transcripts presented GO ontology on biological processes, 181 (181/207, 87.44%) being functionally enriched, and 26 (26/207, 12.56%) did not show enrichment (Table 1).

The distribution of 644 differentially expressed tran-
scripts in the three different GO categories, BP, MF and
CC, is represented in the Venn diagram (Fig. 1).

Many transcripts are associated with more than one
GO ontology. The majority of transcripts correspond to
cellular components (282/644, 47.79%), followed by bio-
logical processes (147/644, 22.83%) and molecular func-
tion (127/644, 19.72%). Fourteen transcripts present all
three categories.

Gene ontology enrichment analysis (Fisher's exact
test) was performed using as "test set" the list of upreg-
ulated (Fig. 2) and downregulated (Fig. 3) differentially
expressed transcripts (DE) and as "reference set" (back-
ground) the *L. infantum* JPCM5 predicted proteome.
FDR < 0.05 was set as the threshold to define the func-
tional enrichment significance.

Gene ontology enrichment analysis of 837 transcripts
upregulated in the Sb<sup>III</sup>-resistant *L. infantum* showed
that the overrepresented terms were related to quorum

sensing (GO:0009372, GO:0052097 and GO:0052106),
host–parasite interaction (GO:0044764, GO:0044114,
GO:0044115 and GO:0044145), protein modifications
(GO:1903320 and GO:0032446), post-translational modi-
fications, protein phosphorylation (GO:0006468) and
protein ubiquitination (GO:0016567), microtubule-based
movement (GO:0007018) and regulation of membrane
lipid distribution (GO:0097035) (Figs. 2, 4).

In contrast, GO enrichment analysis of 96 tran-
scripts downregulated in the Sb<sup>III</sup>-resistant *L. infan-
tum* showed overrepresentation of all GO terms linked
with rRNA processing (GO:0006364), nucleosome
assembly (GO:0034622, GO:0022607, GO:0065003,
GO:0022618, GO:0070925, GO:0042255, GO:0000028,
GO:0006334, GO:0031497, GO:0006333, GO:0065004
and GO:0034063) and maintenance of translational fidel-
ity (GO:1990145) (Figs. 3, 4).

RNA profiling of *L. infantum* (MHOM/BR/74/PP75)

Genes upregulated in the LiSbR line

Out of 431 (431/837, 51.49%) upregulated enriched genes
in the Sb<sup>III</sup>-resistant *L. infantum* line, 124 (124/431,
28.77%) presented GO ontology on the biological pro-
cess (111 enriched and 13 without enrichment), 289
(289/431, 67.05%) genes did not present GO ontology on
the biological process, and 18 (18/431, 4.18%) genes had
not GO ontology associated (Additional file 1: Table S1,
Additional file 2: Table S2 and Additional file 3: Table S3,
respectively).

According to GO enrichment analysis, some terms
related to biological processes were under- or overrepre-
sented in the differentially expressed genes (Fig. 2). The
most representative GO terms were protein phospho-
rylation, microtubule-based movement, protein ubiqu-
uitination, cellular process, quorum sensing involved
in interaction with hosts and others (Fig. 4). Data from
other DE genes up- and downregulated in the LiSbR line
that presented GO enrichment are described on Tables 2
and 3. These data are representative of the complete
results given in Additional file 1: Table S1.

Thirty-seven transcripts belonging to the protein phos-
phorylation category were upregulated in the LiSbR line
(Table 2, Additional file 1: Table S1 and Additional file 2:
Table S2). This group includes five transcripts encod-
ing phosphatidylinositol kinase (PIK) (2.52 to 3.97-fold
upregulated); RAB GTPases (2-fold upregulated); dual
specificity protein phosphatase (DUSP) (8.93-fold upreg-
ulated); protein phosphatase and protein phosphatase
2C (2.28 to 17.52-fold upregulated, respectively); cyclins
10, 11 and CYC2-like (2.27 to 5.56-fold upregulated) and
elongation factor 2 (EF2) (3.23-fold upregulated).
<table><caption>**Table 1** Transcripts differentially expressed between wild-type and Sb<sup>III</sup>-resistant *L. infantum* lines and Gene Ontology (GO) functional enrichment analysis</caption><thead><tr><th rowspan="2"></th><th colspan="2">GO category: biological process</th><th rowspan="2">GO category: cellular component or<br/>molecular function</th><th rowspan="2">Without GO</th><th rowspan="2">Total</th></tr><tr><th>Enriched</th><th>Not enriched</th></tr></thead><tbody><tr><td>Functional annotation<sup>a</sup></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Upregulated</td><td>111 (11.90%)</td><td>13 (1.39%)</td><td>289 (31.94%)</td><td>18 (1.93%)</td><td>431 (46.19%)</td></tr><tr><td>Downregulated</td><td>37 (3.96%)</td><td>0</td><td>6 (0.64%)</td><td>30 (3.21%)</td><td>73 (7.82%)</td></tr><tr><td>Hypothetical proteins<sup>b</sup></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Upregulated</td><td>32 (3.43%)</td><td>13 (1.39%)</td><td>142 (15.22%)</td><td>219 (23.47%)</td><td>406 (43.51%)</td></tr><tr><td>Downregulated</td><td>1 (0.11%)</td><td>0</td><td>0</td><td>22 (2.36%)</td><td>23 (2.79%)</td></tr><tr><td>Total number of DE tran-<br/>scripts identified</td><td>181</td><td>26</td><td>437</td><td>289</td><td>933</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Venn diagram of shared and specific Gene Ontology terms for the differentially expressed transcripts. The 644 differentially expressed genes (FC ≥ 2) of antimony-resistant *L. infantum* were compared and grouped together using the Gene Ontology categories (BP biological process, MF molecular function, CC cellular component). The total amount of shared and specific sequences in each ontology group is depicted in the figure. In addition, 289 differentially expressed genes (FC ≥ 2) were not assigned to any gene ontology category (WGO)</figcaption></figure>

In the microtubule-based movement category, 20 transcripts were upregulated in the LiSbR line (Table 2 and Additional file 1: Table S1), including dyneins (2.04 to 5.9-fold upregulated); ten transcripts encoding kinesins (2.02 to 5.87-fold upregulated); tryptophan-aspartic acid (WD) protein (2-fold upregulated) and tetratricopeptide repeat domain (TPR) (two-fold upregulated).

GO enrichment analysis also showed that transcripts related to protein ubiquitination were upregulated in the LiSbR line. Twenty transcripts (2.03 to 9.14-fold

upregulated in the LiSbR line) were assigned for this category, such as ubiquitin, ubiquitin-transferase, cullin protein and zinc finger-containing proteins. Four transcripts encoding different zinc finger proteins (C3HC4 type—RING finger and FYVE) were 2.15 to 3.77-fold upregulated in the LiSbR line (Table 2 and Additional file 1: Table S1). Glycosomal transporter (GAT3) was 3.39-fold upregulated in the LiSbR line.

Other categories were also present among the upregulated transcripts. Serine palmitoyltransferase, included in the biosynthetic process category, was 2.99-fold upregulated in the LiSbR line (Table 2 and Additional file 1: Table S1). Transcripts encoding ATP-dependent RNA helicase were 2.38 to 3.34-fold upregulated in the LiSbR line (Table 2, Additional file 1: Table S1 and Additional file 2: Table S2) and were included in the ribonucleoprotein complex assembly category. In the stress granule assembly category, one transcript assigned as pumilio protein was 5.59-fold upregulated in the LiSbR line. Four transcripts related to phospholipid-transporting ATPase/P-type were upregulated in the LiSbR line. In the cellular metabolic process category, a transcript encoding a 100 kDa heat shock protein was 2.86-fold upregulated in the LiSbR line. In the categories primary metabolic process, cellular macromolecule biosynthetic process and cellular nitrogen compound metabolic process, DNAJ was found to be 2.13-fold upregulated in the LiSbR line. Many transcripts belonging to ATP-binding cassette (ABC) transporters were 2.2 to 4.6-fold upregulated in the LiSbR line and were included in the category regulation of membrane lipid distribution and phospholipid translocation. Among the transcripts implicated in quorum sensing involved in interaction with host and multi-organism cellular processes, four transcripts of the RNA recognition motif (RRM) were 2.87 to 14.4-fold upregulated in the LiSbR line.

<footnote><sup>a</sup> Transcripts with functional annotation</footnote>

<footnote><sup>b</sup> Transcripts with no assigned function</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Gene Ontology enrichment analysis for the upregulated transcripts. The GO enrichment analysis (Fisher's exact test) was performed using as test set the list of upregulated transcripts and as reference set the *L. infantum* JPCM5 predicted proteome. FDR < 0.05 was set as a threshold to define the functional enrichment significance. The percentage of sequences in each GO category is described in the Y axis. Red bars represent the percentage of sequences classified in each GO term for the "reference set" group, and the blue bars represent the percentage of sequences classified in each GO term for the "test set" group (DE genes)</figcaption></figure>

## Genes downregulated in the LiSbR line

Out of 73 (73/96, 76.01%) enriched transcripts downregulated in the Sb<sup>III</sup>-resistant *L. infantum* line, 37 (37/73, 50.68%) presented GO enrichment on biological process, six (6/73, 8.22%) genes did not present GO ontology on this category, and 30 (30/73, 41.09%) genes had no GO term associated (Additional file 1: Table S1, Additional file 2: Table S2 and Additional file 3: Table S3, respectively). According to GO enrichment analysis, some terms related to biological processes were underrepresented in the DE genes (Fig. 3, Table 3 and Additional file 1: Table S1). The most representative GO terms were ribonucleoprotein complex subunit organization, rRNA processing, ribosome biogenesis, translation and nucleosome assembly.

According to GO enrichment analysis, a group of terms related to ribosomes were overrepresented in the DE dataset. The transcripts encoding ribosomal proteins such as ribosomal proteins 40S and 60S and nucleolar and nuclear proteins are downregulated in the LiSbR line.

Fourteen transcripts encoding ribosomal proteins of the small ribosomal 40S subunit, 40S ribosomal S3a, S4, S15, S16, S17, S18, S19, S21, S23 and S33 were 2.01 to 3.12-fold downregulated in the LiSbR line (Table 3).

In addition, 11 transcripts encoding ribosomal protein components of the 60S subunit of the large ribosomal subunit—60S ribosomal L5, L7a, L11, L13, L18a, L21, L22, L31, L35 and L37—were 2.02 to 3.0-fold downregulated in the LiSbR line (Table 3).

Transcripts encoding the histones H2A, H2B, H3 and H4 were found 2.0 to 2.56-fold downregulated in anti-mony-resistant *L. infantum* line (Table 3, Additional file 1: Table S1 and Additional file 3: Table S3) and were included in the nucleosome assembly category.

## Proteins without GO enrichment for biological process

Some differentially expressed transcripts were not related to any category in the GO enrichment analysis (Additional file 2: Table S2), including, for example: 60S ribosomal L23a (2.07-fold upregulated in the LiSbR line); cytochrome b5 and cytochrome P450 reductase (6.37 and 4.33-fold upregulated in the LiSbR line, respectively); gamma-glutamylcysteine synthetase (2.6-fold upregulated in the LiSbR line); mannosyltransferase (2.54-fold upregulated in the LiSbR line) and two transcripts encoding protein classified as amastin (3.33- and 2.97-fold upregulated in the LiSbR).
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

## Hypothetical proteins

Data from comparative transcriptomic analysis of susceptible and Sb<sup>III</sup>-resistant *L. infantum* lines showed that 429 differentially expressed transcripts were assigned as hypothetical protein without predicted function. From these, 406 transcripts were upregulated

and 23 were downregulated in the LiSbR line (Table 1). Out of 429 DE transcripts, 46 presented GO ontology on biological process (32 functionally enriched and 13 without enrichment), 142 transcripts did not present GO ontology on biological process, and 241 genes had no GO term associated (Additional file 4: Table S4).
<table><caption>**Table 2** Upregulated enriched genes associated to the Gene Ontology biological process category</caption><thead><tr><th>ID</th><th>Description</th><th>Fold change</th><th>P-value</th></tr></thead><tbody><tr><td colspan="4"><strong>Protein phosphorylation</strong></td></tr><tr><td>LINF_320013700</td><td>CYC2-like cyclin—putative<sup>a</sup></td><td>5.56</td><td>3.34E-42</td></tr><tr><td>LINF_340027100</td><td>Dual specificity protein phosphatase—putative</td><td>8.93</td><td>4.34E-107</td></tr><tr><td>LINF_360054700</td><td>Elongation factor-2 kinase-like protein</td><td>3.23</td><td>2.10E-16</td></tr><tr><td>LINF_350046100</td><td>Kinetoplastid kinetochore protein 3—putative</td><td>2.02</td><td>9.69E-13</td></tr><tr><td>LINF_250012200</td><td>Myosin heavy chain kinase c-like protein</td><td>3.37</td><td>3.77E-20</td></tr><tr><td>LINF_300023300</td><td>Phosphatidylinositol kinase—putative<sup>a</sup></td><td>3.97</td><td>4.10E-42</td></tr><tr><td>LINF_360034100</td><td>Protein kinase—putative<sup>a</sup></td><td>3.88</td><td>5.44E-42</td></tr><tr><td>LINF_010011100</td><td>Rab-GTPase-TBC domain-containing protein</td><td>2.13</td><td>1.17E-16</td></tr><tr><td>LINF_340047100</td><td>Target of rapamycin kinase 3<sup>a</sup></td><td>2.74</td><td>2.33E-21</td></tr><tr><td colspan="4"><strong>Microtubule-based movement</strong></td></tr><tr><td>LINF_070010300</td><td>Dynein heavy chain<sup>a</sup></td><td>5.90</td><td>6.51E-54</td></tr><tr><td>LINF_130011900</td><td>Kinesin—putative<sup>a</sup></td><td>5.87</td><td>1.11E-70</td></tr><tr><td>LINF_330038800</td><td>Present in the outer mitochondrial membrane</td><td>4.06</td><td>1.74E-37</td></tr><tr><td>LINF_260024800</td><td>Microtubule-associated protein—putative</td><td>3.73</td><td>1.67E-59</td></tr><tr><td>LINF_270017700</td><td>Intraflagellar transport protein 88—putative</td><td>2.16</td><td>2.53E-19</td></tr><tr><td>LINF_040010500</td><td>Tetratricopeptide repeat</td><td>2.14</td><td>2.06E-19</td></tr><tr><td>LINF_210017400</td><td>WD40 repeat-containing protein</td><td>2.06</td><td>4.35E-21</td></tr><tr><td colspan="4"><strong>Biosynthetic process</strong></td></tr><tr><td>LINF_350008300</td><td>Serine palmitoyltransferase—putative</td><td>2.99</td><td>2.23E-35</td></tr><tr><td>LINF_170011100</td><td>Phenazine biosynthesis-like protein</td><td>2.15</td><td>3.76E-08</td></tr><tr><td colspan="4"><strong>Cellular process</strong></td></tr><tr><td>LINF_360020600</td><td>N-terminal region of chorein</td><td>4.95</td><td>3.27E-59</td></tr><tr><td>LINF_190015000</td><td>FYVE zinc finger-containing protein</td><td>3.66</td><td>2.95E-27</td></tr><tr><td>LINF_270009800</td><td>Glycosomal transporter (GAT3)—putative</td><td>3.39</td><td>1.48E-39</td></tr><tr><td>LINF_240020050</td><td>Multi-drug resistance protein-like</td><td>3.27</td><td>1.11E-33</td></tr><tr><td>LINF_350025400</td><td>Ankyrin repeat protein—putative</td><td>2.67</td><td>1.11E-21</td></tr><tr><td colspan="4"><strong>Cellular component organization</strong></td></tr><tr><td>LINF_130020500</td><td>Phosphoprotein phosphatase—putative</td><td>2.28</td><td>1.71E-11</td></tr><tr><td colspan="4"><strong>Protein ubiquitination</strong></td></tr><tr><td>LINF_350059500</td><td>SPRY domain; HECT-domain (ubiquitin-transferase)<sup>a</sup></td><td>9.14</td><td>4.18E-81</td></tr><tr><td>LINF_360053000</td><td>Ubiquitin protein ligase—putative<sup>a</sup></td><td>7.30</td><td>3.59E-111</td></tr><tr><td>LINF_220017000</td><td>Zinc finger—C3HC4 type (RING finger)</td><td>3.77</td><td>2.20E-40</td></tr><tr><td>LINF_110017600</td><td>Protein transport protein SEC31—putative</td><td>3.57</td><td>6.80E-54</td></tr><tr><td>LINF_160017100</td><td>WW domain-containing protein</td><td>2.63</td><td>5.06E-15</td></tr><tr><td>LINF_160018000</td><td>Cullin family</td><td>2.23</td><td>1.06E-15</td></tr><tr><td colspan="4"><strong>Ribonucleoprotein complex assembly</strong></td></tr><tr><td>LINF_170016300</td><td>ATP-dependent RNA helicase—putative</td><td>2.38</td><td>1.49E-26</td></tr><tr><td colspan="4"><strong>Stress granule assembly</strong></td></tr><tr><td>LINF_180019600</td><td>Pumilio protein 2—putative</td><td>5.59</td><td>6.98E-64</td></tr><tr><td colspan="4"><strong>Cellular metabolic process</strong></td></tr><tr><td>LINF_270018400</td><td>RING-H2 zinc finger</td><td>3.33</td><td>9.20E-24</td></tr><tr><td>LINF_290018600</td><td>Heat shock protein 100 kDa</td><td>2.86</td><td>5.67E-32</td></tr><tr><td colspan="4"><strong>Primary metabolic process; cellular macromolecule biosynthetic process; cellular nitrogen compound metabolic process</strong></td></tr><tr><td>LINF_300022800</td><td>DNAJ domain protein—putative</td><td>2.13</td><td>9.30E-14</td></tr><tr><td colspan="4"><strong>Primary metabolic process; macromolecule metabolic process; nitrogen compound metabolic process</strong></td></tr><tr><td>LINF_120008300</td><td>Myotubularin-related protein—putative</td><td>4.00</td><td>4.38E-52</td></tr><tr><td colspan="4"><strong>Multi-organism cellular process</strong></td></tr><tr><td>LINF_310038500</td><td>Acetyl-CoA carboxylase—putative</td><td>2.74</td><td>5.28E-33</td></tr></tbody></table>
<table><caption>**Table 2 (continued)**</caption><thead><tr><th>ID</th><th>Description</th><th>Fold change</th><th>P-value</th></tr></thead><tbody><tr><td colspan="4">Regulation of membrane lipid distribution; phospholipid translocation</td></tr><tr><td>LINF_340032300</td><td>Phospholipid-transporting ATPase-like protein<sup>a</sup></td><td>3.69</td><td>2.87E-39</td></tr><tr><td colspan="4">Regulation of membrane lipid distribution</td></tr><tr><td>LINF_260032600</td><td>ATP-binding cassette protein subfamily B—member 2—putative</td><td>2.21</td><td>7.28E-13</td></tr><tr><td colspan="4">Quorum sensing involved in interaction with host; multi-organism cellular process; modulation of symbiont involved in interaction with host</td></tr><tr><td>LINF_330022900</td><td>RNA recognition motif. (a.k.a. RRM-RBD or RNP domain)<sup>a</sup></td><td>14.4</td><td>1.07E-153</td></tr><tr><td>LINF_330023000</td><td>RNA-binding protein—putative</td><td>4.01</td><td>1.25E-12</td></tr></tbody></table>

<table><caption>**Table 3** Downregulated enriched genes associated to the Gene Ontology biological process category</caption><thead><tr><th>ID</th><th>Description</th><th>Fold change</th><th>P-value</th></tr></thead><tbody><tr><td colspan="4">Nucleosome assembly</td></tr><tr><td>LINF_100016800</td><td>Histone H3—putative</td><td>2.54</td><td>5.91E-28</td></tr><tr><td>LINF_170019500</td><td>Histone H2B<sup>a</sup></td><td>2.34</td><td>4.29E-09</td></tr><tr><td>LINF_310040900</td><td>Histone H4</td><td>2.15</td><td>6.51E-18</td></tr><tr><td colspan="4">Ribonucleoprotein complex subunit organization; ribosome biogenesis</td></tr><tr><td>LINF_110013300</td><td>40S ribosomal protein S21—putative<sup>a</sup></td><td>2.77</td><td>2.71E-31</td></tr><tr><td>LINF_070010600</td><td>60S ribosomal protein L7a—putative<sup>a</sup></td><td>2.41</td><td>3.03E-17</td></tr><tr><td>LINF_010009200</td><td>Ribosomal protein S7—putative</td><td>2.36</td><td>9.19E-23</td></tr><tr><td>LINF_360008600</td><td>Nuclear protein family a—putative</td><td>2.11</td><td>4.18E-08</td></tr><tr><td colspan="4">Translation</td></tr><tr><td>LINF_290032300</td><td>60S ribosomal protein L13—putative<sup>a</sup></td><td>2.19</td><td>2.28E-16</td></tr><tr><td>LINF_130017200</td><td>40S ribosomal protein S4—putative</td><td>2.02</td><td>8.88E-14</td></tr><tr><td colspan="4">Ribosome biogenesis</td></tr><tr><td>LINF_060009400</td><td>Ribosomal protein L19e—putative<sup>a</sup></td><td>3.14</td><td>1.54E-35</td></tr><tr><td>LINF_350023700</td><td>60S ribosomal protein L5—putative</td><td>3.00</td><td>9.31E-33</td></tr><tr><td>LINF_360015100</td><td>40S ribosomal protein S18—putative</td><td>2.35</td><td>9.75E-17</td></tr><tr><td colspan="4">Ribonucleoprotein complex subunit organization</td></tr><tr><td>LINF_210018200</td><td>40S ribosomal protein S23—putative<sup>a</sup></td><td>2.62</td><td>1.59E-24</td></tr><tr><td>LINF_030007300</td><td>Ribosomal protein L38—putative</td><td>2.24</td><td>1.49E-22</td></tr><tr><td colspan="4">Ribonucleoprotein complex subunit organization; rRNA processing; ribosome biogenesis</td></tr><tr><td>LINF_260021300</td><td>40S ribosomal protein S33—putative</td><td>2.66</td><td>1.95E-25</td></tr><tr><td>LINF_340050700</td><td>Nucleolar protein family a—putative</td><td>2.10</td><td>2.15E-09</td></tr><tr><td colspan="4">Cytoplasmic translation</td></tr><tr><td>LINF_330028300</td><td>60S ribosomal protein L37<sup>a</sup></td><td>2.27</td><td>7.94E-20</td></tr><tr><td>LINF_350009300</td><td>40S ribosomal protein S3A—putative</td><td>2.00</td><td>3.20E-12</td></tr></tbody></table>

According to GO enrichment analysis, some terms related to biological processes were under- or over-represented in the differentially expressed transcripts (Figs. 2, 3, 4). Similar to DE genes, the main terms enriched were microtubule-based movement, protein phosphorylation, protein ubiquitination, quorum sensing involved in interaction with host, ribonucleoprotein complex and others.

## Discussion

Chemotherapy against leishmaniasis remains the main strategy to manage disease control, but several implications regarding the treatment should be considered [11–17]. Pentavalent antimonials are considered one of the main options of treatment; however, these drugs have several toxic side effects and high resistance rates [9, 11–13, 16, 17]. Thus, the comprehension of resistance

<footnote>a Genes with more than one copy in the corresponding GO category. See Additional file 1: Table S1 for complete data</footnote>

<footnote>a Genes with more than one copy in the corresponding GO category. See Additional file 1: Table S1 for complete data</footnote>
molecular mechanisms in *Leishmania* spp. is crucial to identify potential drug targets to prevent or reverse such mechanisms. In this study, RNA Seq has been used successfully to quantify transcript levels of Sb<sup>III</sup>-resistant and wild-type *L. infantum* lines. Our results showed that many pathways upregulated in the antimony-resistant *L. infantum* line are associated to signaling networks, such as kinases and phosphatases; microtubule-based movement, such as dyneins and kinesins; protein ubiquitination; stress response, such as HSP-100 and DNAJ; regulation of membrane lipid distribution, such as ATP-binding cassette; proteins associated to RNA metabolism, such as RNA-binding proteins, pumilio and other proteins involved in important metabolic pathways. Interestingly, our data revealed that the transcripts encoding ribosomal proteins such as 40S and 60S ribosomal proteins, nucleolar and nuclear proteins and histones are downregulated in the antimony-resistant *L. infantum* line. These results show downregulation of genes involved in translation and ribosome biogenesis then modulating important pathways associated with antimony resistance phenotype in *L. infantum*.

Some of the differentially expressed transcripts identified in this study corroborate previous proteomic analysis of antimony-resistant and -susceptible *L. infantum* lines [22]. In addition, some upregulated transcripts identified in this study were also previously investigated by our group, such as gamma-glutamylcysteine synthetase [41] elongation factor 2 [42] and mannosyltransferase [43], and these previous results confirmed that they are associated with antimony resistance phenotype in *Leishmania* spp.

An interesting category demonstrated to have differentially expressed transcripts was the “protein phosphorylation” category. Protein phosphorylation is one of the most important post-translational modifications regulating various signaling processes. GO enrichment analysis showed that 37 transcripts belonging to the protein phosphorylation category are 2.04 to 8.93-fold upregulated in the LiSbR line (Table 2 and Additional file 1: Table S1). Protein kinases are important regulators of many different cellular processes, such as transcriptional control, cell cycle progression, differentiation and response to stress [44, 45]. They represent promising drug targets for trypanosomes and *Leishmania*, since some of them are essential for viability of parasites and have significant sequence differences from mammalian homologues [44].

A comparative proteomic and phosphoproteomic analyses of Sb<sup>III</sup>-resistant and -susceptible lines of *L. braziliensis* identified several potential candidates for biochemical or signaling networks associated with antimony-resistance phenotype in this parasite [22, 23]. In the antimony-resistant *L. infantum* line, we also observed

that different kinases and phosphatases are differentially expressed in this parasite (Additional file 1: Table S1 and Additional file 2: Table S2). RAB GTPases (whose transcripts were shown to be two-fold upregulated in the LiSbR) play a key role in regulation of exocytic and endocytic pathways in eukaryotic cells. This protein was also more abundant in the LiSbR line [22, 46]. It has been shown that RAB GTPases of *L. major* are highly immunogenic in individuals immune to cutaneous and visceral leishmaniasis [47]. *L. donovani* overexpressing RAB6 showed a resistant phenotype by allowing trans-dibenzalacetone-treated parasites to both increase internal thiol levels and enhance MRP pump activity [47].

Elongation factor 2 (EF2), a relevant factor for production of proteins, can be regulated through inhibitory phosphorylation at threonine 56 by EF2 kinase [48]. The transcripts of this enzyme were 3.23-fold upregulated in the LiSbR line. Our group showed that the EF2-overexpressing *L. braziliensis* clone was slightly more resistant to EF2K inhibitor than the WTS line. Surprisingly, this inhibitor increased the antileishmanial effect of Sb<sup>III</sup>, suggesting that this association might be a valuable strategy for leishmaniasis chemotherapy [22, 42].

Other transcripts associated with dephosphorylation, such as protein phosphatase and protein phosphatase 2C, were 17.52 to 2.28-fold upregulated in the LiSbR line, respectively (Additional file 1: Table S1 and Additional file 2: Table S2). Proteomic analysis using these same *L. infantum* lines showed that both enzymes were also more abundant in the Sb<sup>III</sup>-resistant line [22].

The category of “protein ubiquitination” was also a category whose transcripts were differentially expressed in the resistant parasites. Ubiquitination is a crucial process in all eukaryotic organisms. It is involved in several essential functions, such as degradation of denatured proteins, DNA repair, endocytosis, regulation of protein levels, transcription, and apoptosis [49]. Twenty transcripts that are 2.03 to 9.14-fold upregulated in the LiSbR line were assigned to this category, described as: ubiquitin, ubiquitin-transferase (HECT domain—homologous to the E6-AP carboxyl terminus and SPRY domain—SPla and the RYanodine receptor), cullin protein (involved in ubiquitination through participation in multisubunit ubiquitin ligase complexes), zinc finger-containing proteins and others (Table 2 and Additional file 1: Table S1). Similar to our results, antimony-resistant *L. tropica* isolate also showed overexpression of ubiquitin [50]. These data suggest that increased levels of protein ubiquitination may contribute to degradation of oxidized proteins, protecting the parasite against oxidative stress from antimony.

The zinc finger proteins, serine palmitoyltransferase and ATP-dependent RNA helicase, grouped respectively in the categories of “cellular process,” “biosynthetic
process” and “ribonucleoprotein complex assembly,” also had their transcripts differentially expressed in the present work. Zinc finger proteins are RNA-binding proteins involved in many biological processes by binding nucleic acids or participating in transcriptional or translational processes by mediating protein-protein interactions and membrane association [51]. Zinc finger domains in proteins were recently proposed as potential targets for Sb<sup>III</sup> because of the ability of Sb<sup>III</sup> to compete with Zn<sup>II</sup>. A previous study suggested that the interaction of Sb<sup>III</sup> with zinc finger proteins may modulate the pharmacological action of antimonal drugs [10, 14, 52]. In our study, four transcripts encoding different zinc finger proteins (C3HC4 type—RING finger and FYVE) were 2.15 to 3.77-fold upregulated in the LiSbR line (Table 2 and Additional file 1: Table S1). The FYVE domain is a small zinc-binding module that recognizes phosphatidylinositol 3-phosphate, and the majority of these proteins are implicated in membrane trafficking, protein sorting and signaling transduction [53].

Serine palmitoyltransferase catalyzes the first rate-limiting step in the synthetic pathway of *de novo* sphingolipid biosynthesis [54]. This enzyme was 2.99-fold upregulated in the LiSbR line (Table 2 and Additional file 1: Table S1). Metabolomic analysis revealed differences in the phospholipid and sphingolipid contents between antimony-susceptible and -resistant *L. donovani* isolates [55]. According to Zhang and Beverley [56], these two lipid classes are both abundant and critical to virulence and viability in *Leishmania*.

Transcripts encoding ATP-dependent RNA helicase were 2.38 to 3.34-fold upregulated in the LiSbR line (Table 2, Additional file 1: Table S1 and Additional file 2: Table S2). It plays an essential function in RNA metabolism, including RNA degradation, translation, regulation and RNA editing [57, 58]. A member of RNA helicases “DDX3 DEAD-Box” of *Leishmania* have been shown to play a central role in preventing reactive oxygen species-mediated damage and in maintaining mitochondrial protein quality control [58].

Interestingly, four transcripts related to phospholipid-transporting ATPase/P-type ATPase were upregulated in the LiSbR line. These transcripts were grouped into the category “regulation of membrane lipid distribution; phospholipid-translocating.” ATPases are membrane proteins that perform active ion transport across biological membranes, which are found in bacteria and all eukaryotic cells, including *Leishmania* [59]. Fernandez-Prada et al. [60] demonstrated that different point mutations in a P-type ATPase transporter in *L. infantum* are implicated with cross-resistance to miltefosine and amphotericin B.

The categories of “cellular metabolic process” and “primary metabolic process; cellular macromolecule biosynthetic process; cellular nitrogen compound metabolic process” also showed transcripts differentially expressed in parasites resistant to Sb<sup>III</sup> such as HSPs and DNAJ proteins. The HSPs have important functions in folding, secretion, assembly, intracellular localization, regulation and degradation of other proteins [61]. In general, the heat shock response is a homeostatic mechanism that protects cells from the deleterious effects of environmental stress, such as heat and drug exposure [62]. Several authors reported the overexpression of HSPs in antimony-resistant isolates of *L. donovani* [63–65], *L. braziliensis* and *L. infantum* lines [22, 66]. Here, a transcript encoding a 100 KDa heat shock protein was 2.86-fold upregulated in the LiSbR line. HSP100 has the unique capability of recognizing misfolded proteins within an aggregate and actively unfolding them, ultimately disassembling the insoluble structure and delivering substrates into refolding pathways [67].

DNAJ proteins, also known as HSP40s, are crucial partners for HSP70 chaperones, and much of the functional diversity of the HSP70s is driven by this diverse class of cofactors [67]. Here, DNAJ was 2.13-fold upregulated in the LiSbR line. This protein plays a relevant role in the differentiation process from the promastigote to amastigote stage in *L. infantum*, since it suffers a dramatic increase in phosphorylation [68]. Interestingly, HSP40 was also found increased in artemisinin-resistant *L. donovani* [69].

Interestingly, in our study many transcripts belonging to ATP-binding cassette (ABC) transporters were upregulated in the LiSbR line. These transcripts were grouped in the category of “regulation of membrane lipid distribution; phospholipid translocation.” ABC transporters comprise a superfamily of integral membrane proteins involved in the ATP-dependent transport of a variety of molecules across biological membranes, including amino acids, sugars, peptides, lipids, ions and chemotherapeutic drugs [70]. They have been associated with drug resistance in various diseases. In *Leishmania*, the first ABC protein identified was MRPA (multidrug resistance protein, PgpA) [71] which is a member of the ABCC subfamily, able to confer resistance to antimonials by sequestering thiol-metal conjugates into an intracellular vesicle [71, 72]. Our previous data showed an association of chromosomal amplification of *MRPA* gene with the drug resistance phenotype in Sb<sup>III</sup>-resistant *Leishmania* spp. lines [22, 72]. Similarly, it has been shown that *L. infantum* knockout for the *MRPA* gene is more sensitive to Sb [73]. As ABC transporters are important regulators of drug susceptibility, they are excellent candidates for inhibitor design [74].
Since the regulation of gene expression in trypanosomatids occurs largely at post-transcriptional levels, the main control points in gene expression are mRNA degradation and translation [75]. The RNA-binding proteins (RBP) play essential roles in regulating RNA processing, transport, localization, translation and degradation. RBPs contain various structural motifs, such as RNA recognition motif (RRM), dsRNA-binding domain, zinc finger and others [76]. Four transcripts of RRM were 2.87 to 14.4-fold upregulated in the LiSbR line.

Other transcripts differentially expressed in parasites resistant to Sb<sup>III</sup> could not be classified in any GO enrichment for biological processes, such as ribosomal proteins, cytochrome b5, cytochrome P450 reductase, gamma-glutamylcysteine synthetase and mannosyltransferase. Ribosomal proteins play an important role in the translation, and they also regulate cell growth and apoptosis. In our study, the 60S ribosomal L23a, a component of the 60S subunit of the ribosome large subunit, was found 2.07-fold upregulated in the LiSbR line (Additional file 2: Table S2). In agreement with our results, proteomic analysis showed that this protein also was overexpressed in antimony-resistant *L. donovani* isolates [77]. Interestingly, 60S ribosomal L23a-overexpressing *L. donovani* line was more resistant to sodium antimony gluconate (Sb<sup>V</sup>), miltefosine and paromomycin [78].

Cytochrome b5 and cytochrome P450 reductase, which are involved in oxidoreductase activity, were 6.37- and 4.33-fold upregulated in the LiSbR line, respectively (Additional file 2: Table S2). Cytochrome b5 is a flavohemoprotein associated with oxidative reactions such as catabolism of xenobiotics and compounds of endogenous metabolism [79]. Mukherjee et al. [80] observed that *L. major* deficient in cytochrome b5 oxidoreductase domain presents decreasing of linoleate synthesis followed by increased oxidative stress and cell death by apoptosis. Cytochrome P450 reductase is located on the endoplasmic reticulum in many types of cells and is also related to drug metabolism [80].

Gamma-glutamylcysteine synthetase (γ-GCS) is the first enzyme of the glutathione pathway that produces γ-glutamylcysteine, a direct precursor of glutathione [81]. Our results showed that one transcript encoding this enzyme was found 2.6-fold upregulated in the LiSbR line (Additional file 2: Table S2). γ-GCS has been shown to be essential for *L. infantum*, where it confers protection against oxidative stress and Sb<sup>V</sup> [81]. An increase of GSH1 mRNA levels also has been reported in some *L. tarentolae* samples with in vitro-induced resistance to antimony [82] and some Sb<sup>V</sup>-resistant *L. donovani* field isolates [83]. Overexpression of γ-GCS is associated with antimony resistance phenotype in *L. guyanensis* [41].

Glycosylphosphatidylinositol is a surface molecule important for host–parasite interactions. Mannosyl-transferase (GPI-14) is an essential enzyme for adding mannose on the glycosylphosphatidyl group. Our data showed that one transcript encoding this enzyme is 2.54-fold upregulated in the LiSbR line (Additional file 2: Table S2). Interestingly, our group overexpressed the GPI-14 gene in *L. braziliensis* and observed the involvement of the GPI-14 enzyme in the Sb<sup>III</sup> resistance phenotype of *L. braziliensis* [43].

## Conclusions

Transcriptomic profiling represents an important technology for analyzing the global changes in gene expression and regulation of the main metabolic pathways. This study allowed us to compare the transcriptome data from Sb<sup>III</sup>-resistant and wild-type *L. infantum* lines and identify a robust set of transcripts from several biochemical pathways that are associated with the antimony resistance phenotype in this parasite. Overall, our results support the idea that the antimony resistance mechanism in *Leishmania*, similar to other organisms, is complex and multifactorial. The proteins encoded by DE genes may be further evaluated as molecular targets for new drugs against leishmaniasis. In addition, functional studies will be performed to determine the role of some hypothetical proteins and genes with unknown function in the antimony resistance phenotype in *Leishmania*.

## Supplementary information

**Supplementary information** accompanies this paper at https://doi.org/10.1186/s13071-020-04486-4.

**Additional file 1: Table S1.** Enriched genes for biological process category with Gene Ontology assigned terms.

**Additional file 2: Table S2.** Not enriched genes for the biological process category with Gene Ontology assigned terms.

**Additional file 3: Figure S1.** Not enriched genes for the biological process category without Gene Ontology assigned terms.

**Additional file 4: Table S3.** Hypothetical proteins: enriched for the biological process (BP) category, not enriched for BP with and without Gene Ontology assigned terms.

## Abbreviations

ABC transporter: ATP-binding cassette transporter; DE: Differentially expressed; EC<sub>50</sub>: Effective concentration required to decrease growth by 50%; FC: Fold change; GFF: General feature format; GO: Gene ontology; RIN: RNA integrity number; SAM: Sequence alignment map; Sb<sup>V</sup>: Pentavalent antimony; Sb<sup>III</sup>: Trivalent antimony; VL: Visceral leishmaniasis.

## Acknowledgements

The authors thank the FIOCRUZ - Institut Pasteur Program. We thank also the Sequencing Platform of the Institut Pasteur-Paris and the Program for Technological Development in Tools for Health-PDTIS-FIOCRUZ for use of their facilities.
**Authors' contributions**

SMFM, JCR, PP, GFS and DMR conceived, designed and supervised the project. SMFM, JMA, LOG, PP and DBL carried out the experiments. JCR, DMR, LOG, DAL, DBL, FGG, LMO, JPLV and RGD carried out the bioinformatics analyses. SMFM, JCR, DMR, DBL, LOG, AMMS, DAL, PP and GFS wrote the manuscript and analyzed the results. All authors read and approved the final manuscript.

**Funding**

This investigation received financial support from the following agencies: FIOCRUZ-Institut Pasteur Program, Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG—CBB-PPM 00610/15, PPM-00710-15 and APQ-00931-18), Convênio UGA/FAPEMIG, Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq 310104/2018-1 and 304158/2019-4), P3D-Programa de Descoberta e Desenvolvimento de Drogas (PROEP/CNPq/FIOCRUZ 401988/2012-0), INOVA FIOCRUZ/Fundação Oswaldo Cruz and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES)—Finance Code 001. S. M. F. Murta, J. C. Ruiz and D. A. Lima are supported by CNPq (National Council for the Development of Research of Brazil). D. M. Resende is supported by Programa de Pós-Graduação em Ciências da Saúde—IRR and PNPD/CAPES. A. M. M. Santi, L. O. Gonçalves, J. P. L. Velloso, L. M. Oliveira and J. M. Andrade are supported by CAPES and F. G. Guimarães by VPEIC Fiocruz.

**Availability of data and materials**

The datasets supporting the conclusions of this article are included within the article and its additional files. Sequences generated during the present study were deposited at NCBI database under SRA accession numbers SRX2833233, SRX2833324, SRX2833326, SRX2833322, SRX2833327. BioSample accession numbers SAMN07137473, SAMN07137475 and BioProject accession number PRJNA348689.

**Ethics approval and consent to participate**

Not applicable.

**Consent for publication**

Not applicable.

**Competing interests**

The authors declare that they have no competing interests.

**Author details**

¹ Genómina Funcional de Parasitos, Instituto René Rachou, Fiocruz Minas, Belo Horizonte, MG, Brazil. ² Grupo Informática de Biossistemas, Instituto René Rachou, Fiocruz Minas, Belo Horizonte, MG, Brazil. ³ Programa de Pós-graduação em Biologia Computacional e Sistemas, Instituto Oswaldo Cruz, Fiocruz, Rio de Janeiro, RJ, Brazil. ⁴ Universidade Federal do Piauí, Teresina, PI, Brazil. ⁵ Unité Biologie des ARN des Pathogènes Fongiques, Département de Mycologie, Institut Pasteur, Paris, France. ⁶ Unité de Parasitologie moléculaire et Signalisation, Département de Parasitologie et Mycologie, Institut Pasteur, Paris, France.

Received: 7 August 2020 Accepted: 17 November 2020

Published online: 30 November 2020

**References**

1. World Health Organization. Leishmaniasis overview. https://www.who.int/health-topics/leishmaniasis#tab=tab. Leishmaniasis Key facts. https://www.who.int/news-room/factsheets/detail/leishmaniasis. Accessed 25 June 2020.

2. Alvar J, Vélez ID, Bern C, Herrero M, Desjeux P, Cano J, et al. Leishmaniasis worldwide and global estimates of its incidence. PLoS ONE. 2012;7:356–71.

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

4. Guerin PJ, Olliaro P, Sundar S, Boelaert M, Croft SL, Desjeux P, et al. Visceral leishmaniasis: current status of control, diagnosis, and treatment, and a proposed research and development agenda. Lancet Infect Dis. 2002;2:494–501.

5. Shaked-Mishan P, Ulrich N, Ephros M, Zilberstein D. Novel intracellular SbV reducing activity correlates with antimony susceptibility in *Leishmania donovani*. J Biol Chem. 2001;276:3971–6.

6. Berman JD, Gallalee JV, Best JM. Sodium stibogluconate (Pentostam) inhibition of glucose catabolism via the glycolytic pathway, and fatty acid beta-oxidation in *Leishmania mexicana* amastigotes. Biochem Pharmacol. 1987;36:197–201.

7. Wyllie S, Cunningham ML, Fairlamb AH. Dual action of antimonial drugs on thiol redox metabolism in the human pathogen *Leishmania donovani*. J Biol Chem. 2004;279:39925–32.

8. Sereno D, Holzmuller P, Mangot I, Cuny G, Ouaissi A, Lemesre J-L. Antimicrobial-mediated DNA fragmentation in *Leishmania infantum* amastigotes. Antimicrob Agents Chemother. 2001;45:2064–9.

9. Sudhandiran G, Shaha C. Antimicrobial-induced increase in intracellular Ca²⁺ through non-selective cation channels in the host and the parasite is responsible for apoptosis of intracellular *Leishmania donovani* amastigotes. J Biol Chem. 2003;278:25120–32.

10. Frézard F, Silva H, de Castro Pimenta AM, Farrell N, Demicheli C. Greater binding affinity of trivalent antimony to a CCCH zinc finger domain compared to a CCHC domain of kinetoplastid proteins. Met Integr Biometal Sci. 2012;4:433–40.

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

12. Croft SL, Olliaro P. Leishmaniasis chemotherapy-challenges and opportunities. Clin Microbiol Infect. 2011;17:1478–83.

13. Sundar S. Drug resistance in Indian visceral leishmaniasis. Trop Med Int Health. 2001;6:849–54.

14. Gourbal B, Sonuc N, Bhattacharjee H, Legare D, Sundar S, Ouellette M, et al. Drug uptake and modulation of drug resistance in *Leishmania* by an aquaglyceroporin. J Biol Chem. 2004;279:31010–7.

15. Marquis N, Gourbal B, Rosen BP, Mukhopadhyay R, Ouellette M. Modulation in aquaglyceroporin AQP1 gene transcript levels in drug-resistant *Leishmania*. Mol Microbiol. 2005;57:1690–9.

16. Andrade JM, Baba EH, Machado-de-Avila RA, Chavez-Olortegui C, Demicheli CP, Frézard F, et al. Silver and nitrate oppositely modulate antimony susceptibility through aquaglyceroporin 1 in *Leishmania* (Viannia) species. Antimicrob Agents Chemother. 2016;60:4482–9.

17. Haimeur A, Brochu C, Genest P-A, Papadopoulou B, Ouellette M. Amplification of the ABC transporter gene PGPA and increased trypanothione levels in potassium antimony tartrate (SbIII) resistant *Leishmania tarentolae*. Mol Biochem Parasitol. 2000;108:131–5.

18. Légaré D, Richard D, Mukhopadhyay R, Stierhof Y-D, Rosen BP, Haimeur A, et al. The *Leishmania* ABC protein PGPA is an intracellular metal-thiol transporter ATPase. J Biol Chem. 2001;276:26301–7.

19. Perea A, Manzano JI, Castanys S, Gamarro F. The LABCG2 transporter from the protozoan parasite *Leishmania* is involved in antimony resistance. Antimicrob Agents Chemother. 2016;60:3489–96.

20. 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:1514.

21. do Monte-Neto RL, Coelho AC, Raymond F, Légaré D, Corbeil J, Melo MN, et al. Gene expression profiling and molecular characterization of antimony resistance in *Leishmania amazonensis*. PLoS Negl Trop Dis. 2011;5:1167.

22. 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.

23. Moreira DDS, Pescher P, Laurent C, Lenormand P, Späth GF, Murta SMF. Phosphoproteomic analysis of wild-type and antimony-resistant *Leishmania braziliensis* lines by 2D-DIGE technology. Proteomics. 2015;15:2999–3019.

24. Maretti-Mira AC, Bittner J, Oliveira-Neto MP, Liu M, Kang D, Li H, et al. Transcriptome patterns from primary cutaneous *Leishmania braziliensis* infections associated with eventual development of mucosal disease in humans. PLoS Negl Trop Dis. 2012;6:e1816.

25. Dillon LAL, Okrah K, Hughitt VK, 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.
26. Rastrojo A, Carrasco-Ramiro F, Martín D, Crespillo A, Reguera RM, Aguado B, et al. The transcriptome of *Leishmania major* in the axenic promastigote stage: transcript annotation and relative expression levels by RNA-seq. BMC Genom. 2013;14:223.

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

28. 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 Sb III in *Leishmania braziliensis* and *Leishmania panamensis*. Sci Rep. 2019;9:9485.

29. Liarte DB, Murta SMF. Selection and phenotype characterization of potassium antimony tartrate-resistant populations of four New World *Leishmania* species. Parasitol Res. 2010;107:205–12.

30. la Fuente SG, Peiró-Pastor R, Rastrojo A, Moreno J, Carrasco-Ramiro F, Requena JM, et al. Resequencing of the *Leishmania infantum* (strain JPCM5) genome and de novo assembly into 36 contigs. Sci Rep. 2017;7:18050.

31. Schmieder R, Edwards R. Quality control and preprocessing of metagenomic datasets. Bioinform Oxf Engl. 2011;27:863–4.

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

33. Dobin A, Gingeras TR. Mapping RNA-seq reads with STAR. Curr Protoc Bioinform. 2015;15:11–4.

34. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

35. Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29:24–6.

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

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

38. Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21:3674–6.

39. Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, et al. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36:3420–35.

40. The Gene Ontology Consortium. Gene Ontology Consortium: going forward. Nucleic Acids Res. 2015;43:1049–56.

41. Fonseca MS, Comini MA, Resende BV, Santi AMM, Zoboli AP, Moreira DS, et al. Ornithine decarboxylase or gamma-glutamylcysteine synthetase overexpression protects *Leishmania* (*Vianna*) *guyanensis* against antimony. Exp Parasitol. 2017;175:36–43.

42. Moreira DS, Murta SM. Involvement of nucleoside diphosphate kinase b and elongation factor 2 in *Leishmania braziliensis* antimony resistance phenotype. Parasites Vectors. 2016;9:641.

43. Ribeiro CV, Rocha BFB, Moreira DDS, Peruhype-Magalhães V, Murta SMF. Mannosyltransferase (GPI-14) overexpression protects promastigote and amastigote forms of *Leishmania braziliensis* against trivalent antimony. Parasites Vectors. 2019;12:60.

44. Naula C, Parsons M, Mottram JC. Protein kinases as drug targets in trypanosomes and *Leishmania*. Biochim Biophys Acta BBA Proteins Proteom. 2005;1754:151–9.

45. Ardito F, Giuliani M, Perrone D, Troiano G, Muzio LL. The crucial role of protein phosphorylation in cell signaling and its use as targeted therapy. Int J Mol Med. 2017;40:271–80.

46. Chamakh-Ayari R, Chenik M, Chakroun AS, Bahi-Jaber N, Aoun K, Meddeb-Garnaoui A. *Leishmania major* large RAB GTPase is highly immunogenic in individuals immune to cutaneous and visceral leishmaniasis. *Parasites Vectors. 2017;10:185.

47. Chauhan IS, Rao GS, Singh N. Enhancing the copy number of Ldrab6 gene in *Leishmania donovani* parasites mediates drug resistance through drug thiol conjugate dependent multidrug resistance protein A (MRPA). Acta Trop. 2019;199:105158.

48. Ryazanov AG, Davydova EK. Mechanism of elongation factor 2 (EF-2) inactivation upon phosphorylation phosphorylated EF-2 is unable to catalyze translocation. FEBS Lett. 1989;251:187–90.

49. Marín I. RBR ubiquitin ligases: diversification and streamlining in animal lineages. J Mol Evol. 2009;69:54–64.

50. Kazemi-Rad E, Mohebali M, Khadem-Erfan MB, Hajjaran H, Hadighi R, Khamesipour A, et al. Overexpression of ubiquitin and amino acid permease genes in association with antimony resistance in *Leishmania tropica* field isolates. Korean J Parasitol. 2013;51:413–9.

51. Laity JH, Lee BM, Wright PE. Zinc finger proteins: new insights into structural and functional diversity. Curr Opin Struct Biol. 2001;11:39–46.

52. Demicheli C, Frézard F, Mangrum JB, Farrell NP. Interaction of trivalent antimony with a CCHC zinc finger domain: potential relevance to the mechanism of action of antimonial drugs. Chem Commun. 2008;39:4828–30.

53. Kutateladze T. Phosphatidylinositol 3-phosphate recognition and membrane docking by the FYVE domain. Biochim Biophys Acta BBA Mol Cell Biol Lipids. 2006;1761:868–77.

54. Denny PW, Goulding D, Ferguson MAJ, Smith DF. Sphingolipid-free *Leishmania* are defective in membrane trafficking, differentiation and infectivity. Mol Microbiol. 2004;52:313–27.

55. t'Kindt R, Scheltema RA, Jankevics A, Brunker K, Rijal S, Dujardin J-C, et al. Metabolomics to unveil and understand phenotypic diversity between pathogen populations. PLoS Negl Trop Dis. 2010;4:904.

56. Zhang K, Beverley SM. Phospholipid and sphingolipid metabolism in *Leishmania*. Mol Biochem Parasitol. 2010;170:55–64.

57. Moller JV, Juul B, le Maire M. Structural organization, ion transport, and energy transduction of P-type ATPases. Biochim Biophys Acta. 1996;1286:1–51.

58. Padmanabhan PK, Zghidi-Abouzid O, Samant M, Dumas C, Aguiar BG, Estaquier J, et al. DDX3 DEAD-box RNA helicase plays a central role in mitochondrial protein quality control in *Leishmania*. Cell Death Dis. 2016;7:2406.

59. Ouameur A, Girard I, Legare D, Ouellette M. Functional analysis and complex gene rearrangements of the folate/biopterin transporter (FBT) gene family in the protozoan parasite *Leishmania*. Mol Biochem Parasitol. 2008;162:155–64.

60. Fernandez-Prada C, Vincent IM, Brotherton M-C, Roberts M, Roy G, Rivas L, et al. Different mutations in a P-type ATPase transporter in *Leishmania* parasites are associated with cross-resistance to two leading drugs by distinct mechanisms. PLoS Negl Trop Dis. 2016;10:0005171.

61. Young JC, Agashe VR, Siegers K, Hartl FU. Pathways of chaperone-mediated protein folding in the cytosol. Nat Rev Mol Cell Biol. 2004;5:781–91.

62. Folgueira C, Requena JM. A postgenomic view of the heat shock proteins in kinetoplastids. FEMS Microbiol Rev. 2007;31:359–77.

63. Vergnes B, Gourbal B, Girard J, Sundar S, Drummelsmith J, Ouellette M. A Proteomics screen implicates HSP83 and a small kinetoplastid calpain-related protein in drug resistance in *Leishmania donovani* clinical field isolates by modulating drug-induced programmed cell death. Mol Cell Proteom. 2007;6:88–101.

64. 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.

65. Kumar D, Singh R, Bhandari V, Kulshrestha A, Negi NS, Salotra P. Biomarkers of antimony resistance: need for expression analysis of multiple genes to distinguish resistance phenotype in clinical isolates of *Leishmania donovani*. Parasitol Res. 2012;111:223–30.

66. Brotherton M-C, Bourassa S, Légaré D, Poirier GG, Droit A, Ouellette M. Quantitative proteomic analysis of amphotericin B resistance in *Leishmania infantum*. Int J Parasitol Drugs Drug Resist. 2014;4:126–32.

67. Requena JM, Montalvo AM, Fraga J. Molecular chaperones of *Leishmania*: central players in many stress-related and -unrelated physiological processes. BioMed Res Int. 2015;2015:1–21.

68. Tsigankov P, Gherardini PF, Helmer-Citterich M, Späth GF, Myler PJ, Zilberstein D. Regulation dynamics of *Leishmania* differentiation: deconvoluting signals and identifying phosphorylation trends. Mol Cell Proteom. 2014;13:1787–99.

69. Verma A, Ghosh S, Salotra P, Singh R. Artemisinin-resistant *Leishmania* parasite modulates host cell defense mechanism and exhibits altered expression of unfolded protein response genes. Parasitol Res. 2019;118:2705–13.

70. Higgins CF. ABC transporters: from microorganisms to man. Annu Rev Cell Biol. 1992;8:67–113.
71. Ouellette M, Légaré D, Haimeur A, Grondin K, Roy G, Brochu C, et al. ABC transporters in *Leishmania* and their role in drug resistance. Drug Resist Updates. 1998;1:43–8.

72. Moreira DS, Monte Neto RL, Andrade JM, Santi AMM, Reis PG, Frézard F, et al. Molecular characterization of the MRPA transporter and antimony uptake in four New World *Leishmania* spp. susceptible and resistant to antimony. Int J Parasitol Drugs Drug Resist. 2013;3:143–53.

73. 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.

74. Singh S, Mandlik V. Structure based investigation on the binding interaction of transport proteins in leishmaniasis: insights from molecular simulation. Mol Biosyst. 2015;11:1251–9.

75. Clayton C, Shapira M. Post-transcriptional regulation of gene expression in trypanosomes and leishmanias. Mol Biochem Parasitol. 2007;156:93–101.

76. Lunde BM, Moore C, Varani G. RNA-binding proteins: modular design for efficient function. Nat Rev Mol Cell Biol. 2007;8:479–90.

77. Kumar A, Sisodia B, Misra P, Sundar S, Shasany AK, Dube A. Proteome mapping of overexpressed membrane-enriched and cytosolic proteins in sodium antimony gluconate (SAG) resistant clinical isolate of *Leishmania donovani*: overexpressed proteins in SAG resistant *L. donovani*. Br J Clin Pharmacol. 2010;70:609–17.

78. 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:2527.

79. Schenkman JB, Jansson I. The many roles of cytochrome b5. Pharmacol Ther. 2003;97:139–52.

80. Mukherjee S, Santara SS, 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–5003.

81. Mukherjee A, Roy G, Guimond C, Ouellette M. The γ-glutamylcysteine synthetase gene of *Leishmania* is essential and involved in response to oxidants. Mol Microbiol. 2009;74:914–27.

82. Guimond C, Trudel N, Brochu C, Marquis N, El Fadili A, Peytavi R, et al. Modulation of gene expression in *Leishmania* drug resistant mutants as determined by targeted DNA microarrays. Nucleic Acids Res. 2003;31:5886–96.

83. Mukherjee A, Padmanabhan PK, Singh S, Roy G, Girard I, Chatterjee M, et al. Role of ABC transporter MRPA, γ-glutamylcysteine synthetase and ornithine decarboxylase in natural antimony-resistant isolates of *Leishmania donovani*. J Antimicrob Chemother. 2007;59:204–11.

## Publisher's Note

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

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_8.png" /></figure>