Research ArticleLeishmaniasis

Variable gene expression and parasite load predict treatment outcome in cutaneous leishmaniasis

See allHide authors and affiliations

Science Translational Medicine  20 Nov 2019:
Vol. 11, Issue 519, eaax4204
DOI: 10.1126/scitranslmed.aax4204

Calling healing lesions

Cutaneous leishmaniasis results in skin sores that can be difficult to treat, resulting in further complications. Amorim et al. performed dual RNA sequencing of host and parasite gene expression in pretreatment lesional skin biopsies from two cohorts of patients infected with Leishmania braziliensis. A prognostic signature comprising expression of three cytolytic genes plus pathogen load predicted treatment response in both cohorts and could potentially be used to triage patients who are unlikely to respond to conventional treatment as candidates for alternate therapies. This study also provides evidence that inflammasome-mediated IL-1β production influences treatment outcome in patients, strengthening the argument that therapies targeting these components could treat cutaneous leishmaniasis.


Patients infected with Leishmania braziliensis develop chronic lesions that often fail to respond to treatment with antiparasite drugs. To determine whether genes whose expression is highly variable in lesions between patients might influence disease outcome, we obtained biopsies of lesions from patients before treatment with pentavalent antimony and performed transcriptomic profiling on these clinical samples. We identified genes that were highly variably expressed between patients, and the variable expression of these genes correlated with treatment outcome. Among the most variable genes in all the patients were components of the cytolytic pathway, and the expression of these genes correlated with parasite load in the skin. We demonstrated that treatment failure was linked to the cytolytic pathway activated during infection. Using a host-pathogen marker profile of as few as three genes, we showed that eventual treatment outcome could be predicted before the start of treatment in two separate cohorts of patients with cutaneous leishmaniasis (n = 21 and n = 25). These findings raise the possibility of point-of-care diagnostic screening to identify patients at high risk of treatment failure and provide a rationale for a precision medicine approach to drug selection in cutaneous leishmaniasis. This work more broadly demonstrates the value of identifying genes of high variability in other diseases to better understand and predict diverse clinical outcomes.


Cutaneous leishmaniasis (CL) is characterized by the development of cutaneous lesions that may resolve in several months or, in some cases, lead to the development of metastatic lesions (1, 2). The response to drug treatment is highly variable, and antimony, the current first-line treatment for CL in Brazil, is not only toxic but also has a high failure rate (16). Because the current drug therapy for leishmaniasis, pentavalent antimony, is often unsuccessful, defining the pathways leading to disease could potentially identify targets for host-directed therapies. Studies with patients combined with parallel studies in experimental murine models have revealed that tissue destruction in CL is initiated by cytolytic CD8+ T cells that lead to inflammasome activation and consequent interleukin-1β (IL-1β) production (713). IL-1β then promotes increased inflammation, enhancing disease severity. Blocking either Nlrp3 or IL-1β protects mice against CD8+ T cell–mediated immunopathology (8), indicating that these proteins might be important targets in host-directed therapies for the disease. However, despite overwhelming direct evidence that this cytolytic pathway promotes increased disease in experimental models (7, 8, 1421) and indirect evidence in patients (8, 1012, 22), a direct demonstration that this pathway influences disease outcome in patients is lacking.

In the past decade, advanced sequencing platforms have accelerated the identification of both host and pathogen gene expression signatures associated with immunopathogenesis, clinical progression, and response to treatment of infectious diseases (2327). We predicted that an assessment of genes that are differentially expressed between infected patients might identify genes associated with disease outcome that could be useful as targets for therapeutic design or as predictors of treatment failure. Given that cytolytic genes are induced early in CL lesions (9), we specifically hypothesized that variations in the magnitude of expression of such genes might influence disease outcome and provide the potential markers to identify patients who may fail conventional therapy.

To address this hypothesis, we performed RNA sequencing (RNA-seq) of lesion biopsies taken from patients with CL at the initiation of treatment with pentavalent antimony and identified highly variable genes up-regulated relative to healthy skin. After RNA-seq of biopsies of two distinct patient cohorts and statistical filter strategies, we focused on a set of genes that were up-regulated in lesions from patients who did not respond to treatment, which included genes involved in cytolysis. Using dual RNA-seq and quantitative polymerase chain reaction (qPCR) to identify both parasite and host transcripts, we discovered that small differences in the parasite load correlated with cytolytic gene expression and treatment outcome, suggesting that parasites were driving this cytolytic pathway. Together, these studies demonstrate that differences in the expression of the cytolytic pathway leading to IL-1β production influence disease outcome and that variability in the expression of a small number of genes, as well as parasite load, can predict responsiveness to treatment.


Transcriptional profiling indicates that cytotoxicity-related pathways are enriched in Leishmania braziliensis lesions

To evaluate the transcriptional signatures associated with disease, we carried out RNA-seq on skin biopsies collected from 21 patients infected with Leishmania braziliensis before pentavalent antimony treatment and 7 uninfected endemic controls. Transcriptional profiles from lesions were distinguishable from healthy skin by principle components analysis (Fig. 1A), with 4255 genes identified as differentially expressed between the two groups (Fig. 1B). Confirming previous findings (9), we identified two cytotoxicity-related pathways included in the top 10 enriched pathways in lesions by gene set enrichment analysis (GSEA) (Fig. 1C) and found that many genes associated with cytolysis, including GZMB, PRF1, GNLY, and IL1B, exhibited significant enrichment over healthy controls [false discovery rate (FDR) < 0.01] (Fig. 2A). Genes associated with B cell activation, phagocytosis, complement activation, and Fc gamma receptor signaling were among the most highly expressed genes in lesions (fig. S1) (28).

Fig. 1 Cytotoxicity-related pathways are enriched in CL lesions relative to healthy skin.

(A) Principal components analysis showing principal component 1 (PC1) and PC2 for RNA-seq data from healthy participants (n = 7) and CL lesions (n = 21) from this study. A PERMANOVA statistical test was used to compare the two clinical groups. (B) A volcano plot showing up-regulated (red; n = 2037) and down-regulated genes (blue; n = 2218) in biopsies from patients relative to healthy participants (FDR < 0.01 and logFC > 1). (C) Gene set enrichment analysis (GSEA) showing two cytotoxicity-related pathways included in the top 10 enriched pathways in lesions relative to healthy skin (HS). NES, normalized enrichment score; adj.P.value, adjusted P value; TCR, T cell receptor; KEGG, Kyoto Encyclopedia of Genes and Genomes; pathway interaction database (PID).

Fig. 2 Identification of highly variable gene expression signatures in CL lesions.

(A) Box plots showing expression of GZMB, PRF1, GNLY, and IL1B in healthy skin (n = 7; gray) and CL lesions (n = 21; red). FC comparing the sample with the highest versus lowest expression for each of the four genes in CL lesions. (B) Coefficient of variation (cV) for the 2037 up-regulated genes in CL lesions relative to healthy participants (healthy skin). Histograms at the top and right show scatter plot distributions. The dashed line indicates an arbitrary cutoff (cV = 0.65) for the most variable genes defined here as ViTALs, with 250 ViTAL genes above the line. Colored points indicate statistically significant pathways (FDR < 0.01) within ViTALs. GZMB, PRF1, GNLY, and IL1B genes are labeled. Legends are indicated inside the plot identifying enriched pathways or clusters of ViTALs. CPM, counts per million; ViTALs, variable transcripts associated with lesions.

Highly variably expressed genes are present in patient L. braziliensis lesions

We previously reported that patient-to-patient variability in gene expression during CL could not be explained by the stage of disease (9). Because cytolysis promotes pathology through the release of IL-1β (8), we hypothesized that variability in the expression of genes central to these pathways may affect clinical outcome in patients. We found that genes associated with cytotoxicity (GZMB, PRF1, and GNLY) and IL1B, all of which promote pathology in CL, varied substantially between patients (Fig. 2A). To explore the variability of gene expression more broadly in patients with CL, we calculated the coefficient of variation (cV) for all 2037 genes up-regulated in lesions (Fig. 2B). We defined the top 250 most variable genes (cV ≥ 0.65) as “variable transcripts associated with lesions” (ViTALs). Gene Ontology (GO) analysis of ViTALs identified an enrichment of genes involved in “B cell responses” (cluster 1, yellow) and “cell chemotaxis signaling-related” genes (cluster 2, blue) (FDR < 0.05; data file S1). GZMB, PRF1, GNLY, and IL1B (cluster 3, red) also exhibited highly variable expression between patients (Fig. 2B). In contrast, genes that were induced in lesions but showed low variability among patients (cV < 0.35) were enriched for type I and II interferon signaling, antigen processing and presentation, and histone-related gene expression, suggesting that induction of these pathways is common to all lesions but is not associated with clinical outcome.

Variability in cytotoxicity-related gene expression is associated with treatment outcome

Because biopsies were collected before treatment, we could only stratify the samples after monitoring the patient’s treatment outcome. Fourteen individuals had lesions that resolved within 90 days of treatment, whereas the lesions of seven individuals failed to resolve in this time period and required additional rounds of chemotherapy. These patient groups are hereafter referred to as “cured” and those for whom treatment “failed,” respectively (Fig. 3A). None of the 250 ViTALs were associated with patient age, delayed-type hypersensitivity response, or lesion size (Spearman correlation test, FDR > 0.05) nor did these clinical covariates show any correlation with sex or treatment outcome (Mann-Whitney U test, P > 0.05) (fig. S2). In contrast, 31 of 250 ViTALs (12.4%) were differentially expressed between patients for whom treatment failed and those that were cured after the first round of pentavalent antimony treatment (Fig 3B). To determine whether any of these genes also correlated with treatment failure in a second dataset, we performed the same analysis on a cohort that we previously studied (hereafter referred to as “2016 dataset”) (28). Eight of the 31 ViTALs significantly correlated with treatment failure in both the current dataset and the 2016 dataset [log2 fold change (logFC) ≥1 and P ≤ 0.05). These included key components of the cytolytic machinery (GZMB, PRF1, and GNLY); genes associated with cytolytic T cells (UNC13A, APOBEC3A, and KIR2DL4) (2932); C-C motif chemokine ligand 4 (CCL4), which binds to CCR5 and is expressed on CD8+ T cells; and interferon-γ. Consistent with these results, GO enrichment analysis carried out on all 250 ViTALs in our current dataset showed significant enrichment of several pathways, including “graft-versus-host disease” and “natural killer (NK) cell–mediated cytotoxicity” (FDR < 0.0001). Similarly, GO analysis of 595 ViTALs identified in the 2016 dataset (cV < 0.19) confirmed enrichment of the same pathways (FDR < 0.02). Together, these data identify a key set of transcripts and pathways whose expression is associated with treatment failure across multiple studies. We used CombiROC analysis (33) to determine the number of markers needed to predict pentavalent antimony treatment outcome and found in the current cohort of patients that a combination of three genes—GNLY, IL1B, and IFNG—produced an area under the curve of 0.86 to 0.96, indicating that these genes predict treatment outcome in this group with 86 to 96% accuracy (fig. S3).

Fig. 3 Expression of cytotoxicity-related genes is associated with treatment outcome.

(A) The chart defines the overall steps of study design from day 0 to day 90. A representative photograph of a lesion shows in a white circle where the 4-mm punch biopsies were collected. At day 90, patients with complete re-epithelialization of lesions were considered cured (n = 14), and those with active lesions were determined to have failed therapy and subsequently underwent an additional round of antimony treatment (n = 7). (B) Volcano plot showing 250 ViTALs according to treatment outcome. Dashed lines indicate P < 0.05 and P < 0.01. Each point indicates a gene, and point colors reflect the cluster assignment from Fig. 2. Purple text identifies genes found to be differentially overexpressed between antimony treatment outcomes (P < 0.05) in an independent study (28).

Variable expression of cytotoxic genes is associated with increased abundance of CD8+ effector T cells and NK cells in lesions of patients who do not respond to treatment

Variability in expression of cytotoxic genes could be explained by either a difference in the per cell expression of these genes between patients or differences in the number of cytotoxic cells present in lesions from different patients. To distinguish between these two possibilities, we used ImmQuant (34) to infer the cell type composition of lesions using cell type–specific marker genes present in our RNA-seq data (Fig. 4). ImmQuant analysis identified high variability in the predicted relative abundance of CD8+ T cells and NK cells between patients with CL (Fig. 4A). Furthermore, patients who did not respond to treatment exhibited elevated proportions of CD8+ T cells and NK cells but not CD4+ cells (Fig. 4B). There was a positive correlation between the expression of ViTAL cytotoxic genes (GZMB, PRF1, and GNLY) and the predicted relative abundance of CD8+ T cells and NK cells (ρ > 0.6 and P < 0.002) (Fig. 4C). These results suggest that the increase in cytolytic genes is likely because of increased recruitment of cytotoxic cells to the lesion. To confirm this conclusion, we assessed the percentage of granzyme B (GzmB)–expressing CD8+ T cells in patient lesions before treatment using flow cytometry and observed a significant positive correlation between the expression of CD8+ T cells and GzmB (ρ = 0.75 and P = 0.003) (Fig. 4D). With one exception, the patients for whom treatment failed had the highest percentage of GzmB+ CD8+ T cells in their lesions (Fig. 4D), which was consistent with our RNA-seq results. CCL3 and CCL4 were up-regulated in lesions from patients who remained uncured, suggesting these genes might be responsible for increased recruitment of cytotoxic cells to the site of the infection. Consistent with this notion, we observed a positive Spearman’s correlation between the expression of CCL3, CCL4, and CCR5 (the receptor for CCL3 and CCL4) and the predicted relative abundance of CD8+ T cells and NK cells (ρ > 0.6 and P < 0.003) (Fig. 4C). These findings reinforce the role of CD8+ T cells as mediators of pathology during CL infection and implicate the recruitment of cytolytic CD8+ T cells in treatment failure.

Fig. 4 Elevated expression of cytotoxic genes associates with higher predicted abundance of CD8+ effector T cells and NK cells in lesions from patients unresponsive to treatment.

ImmQuant and the DMAP database were used to predict the relative abundance of cell subsets in biopsies based on RNA-seq data. (A) Relative cell abundance scores (ImmQuant scores) of CD4+ effector T cells (CD4+CD62L+CD45RA), CD8+ effector memory T cells (CD8+CD62LCD45RA), CD8+ effector T cells (CD8+CD62LCD45RA+), and NK cells (CD56+CD16+CD3) between samples (n = 21), where each point corresponds to an individual patient sample. Mean and SD are shown. (B) ImmQuant scores between the patients who did not respond to treatment versus patients who were cured. Mann-Whitney U tests were used for statistical analysis. (C) Correlation matrix between the ImmQuant scores of CD4+ and CD8+ effector cells and NK cells between samples and the expression of GZMB, PRF1, GNLY, CCL3, CCL4, and CCR5. (D) Cells isolated from an independent set of lesions of patients infected L. braziliensis were stained directly ex vivo for flow cytometric analysis before treatment with pentavalent antimony. Each point represents a patient, and the red ellipse clusters patients (95% confidence) for whom antimony treatment was unsuccessful. Spearman’s test was used for correlation statistical analysis, and ρ values are represented. ρ, Spearman’s rho correlation coefficient. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Parasite load correlates with gene expression variability

Lesions from patients infected with L. braziliensis have been reported to contain very few parasites (35), and consistent with this, we previously reported that L. braziliensis transcripts could only be identified in a subset of patients by RNA-seq (28). Similarly, in this study, we detected parasite transcripts in only 13 of 21 samples by dual RNA-seq of host and parasite transcripts (Fig. 5A). To better quantify the parasites in lesions, we used qPCR (fig. S4) and found that 19 of the 21 lesion biopsies were positive by qPCR for L. braziliensis 18S ribosomal RNA (rRNA) (Fig. 5B). We found a strong positive correlation between the abundance of parasite transcripts detected by RNA-seq and absolute number of parasites determined by qPCR (ρ = 0.90 and P <0.0001) (Fig. 5C). We next asked whether parasite load was associated with the cytotoxic signature seen in lesions from patients whose lesions did not respond to treatment and found a strong correlation between parasite number measured by either qPCR or RNA-seq and both the expression of ViTALs (Fig. 5D and fig. S5) and the predicted relative abundance of CD8+ T cells and NK cells by ImmQuant (Fig. 5E). Highly expressed genes with limited variability between patients (cV < 0.65), such as MMP2, STAT1, CXCL10, and CD19, did not correlate with parasite load (ρ < 0.29 and P > 0.23) (fig. S6). These results suggest that parasite load contributes to the increased recruitment of cytotoxic cells to the site of infection and, consequently, the high abundance of transcripts from cytotoxic genes detected in lesions.

Fig. 5 Parasite load correlates with gene expression variability and increased abundance of CD8+ T cells and NK cells.

(A) Relative number of L. braziliensis transcripts detected in RNA-seq data from CL lesions (n = 21) and healthy skin (n = 7). The horizontal line marks the median. (B) Absolute number of parasites detected in 4-mm punch biopsies from lesions and healthy skin, as measured by qPCR. Filled points represent samples with detectable parasite transcripts based on RNA-seq analysis shown in (A). Spearman’s correlation between the absolute number of parasites quantified in each lesion biopsy by qPCR and (C) abundance of L. braziliensis transcripts by RNA-seq or (D) expression of ViTALs or (E) ImmQuant scores for CD8+ effector T cells and NK cells. Two samples with undetectable parasites were excluded from the correlation plots. *P < 0.05, **P < 0.01, and ***P < 0.001.

Parasite load predicts treatment outcome

Given the association between ViTALs and treatment failure and ViTALs and parasite load, we hypothesized that parasite load would predict treatment failure. We found that lesions from patients who were unresponsive to the first round of treatment had more parasites compared to the lesions from patients who were cured, as measured by qPCR (median of 34,370 parasites versus 17,805 parasites, respectively, P < 0.001) (Fig. 6A) or by RNA-seq [median of 5198 counts per million (CPM) versus 78 CPM, P < 0.0001] (Fig. 6B, left). We confirmed these results by applying the same analytic pipeline used in the current study to the 2016 RNA-seq dataset (Fig. 6B, right) (28). Last, Kaplan-Meier analysis revealed that patients with biopsies containing more than 32,000 parasites per 4-mm biopsy took nearly 2 months longer to heal than those with less than 32,000 parasites (P = 0.01; Fig. 6C). Receiver operating characteristic (ROC) analysis showed that parasite load in lesions from CL lesions by itself predicts treatment outcome with 96% accuracy (fig. S3). These data suggest that the parasite load may have a direct and proportional effect on the recruitment of cytotoxic cells, thereby affecting how patients respond to treatment.

Fig. 6 Parasite load predicts treatment outcome.

(A) Absolute number of parasites per 4-mm biopsy and (B) relative abundance of L. braziliensis transcripts in patients who were cured versus those for whom treatment was unsuccessful for the current dataset and the 2016 dataset (28). (C) Kaplan-Meier curve showing time to cure from initial clinical screening. Patients were grouped by those with less (blue) or more (red) than 3.2 × 104 parasites per biopsy, as measured by qPCR. **P < 0.01, log-rank (Mantel-Cox) test. **P < 0.01, ***P < 0.001, and ****P < 0.0001, Mann-Whitney U test.


Patients infected by L. braziliensis develop chronic ulcerated lesions that often fail to respond to drug treatment (3). By analyzing gene expression in lesions from these patients, we identified ViTALs that correlated with treatment failure, including several components of the cytolytic machinery (GZMB, PRF1, and GNLY). We and others previously found that cytolysis is a major component of cutaneous lesions (712), and from a combination of human and murine studies, we discovered that cytolysis leads to NLRP3 activation, IL-1β production, and subsequent severe inflammation (8). We found that blocking components of this pathway, such as NLRP3 or IL-1β, ameliorated severe disease in murine models (8). Although these studies implicated cytolysis and IL-1β production in the severity of CL, here, we directly demonstrated that disease outcome in patients in this case, as assessed by treatment failure, associates with the degree of ongoing cytolysis within lesions. Thus, these results substantially strengthen the idea that therapies directed at blocking the NLRP3 inflammasome and IL-1β in combination with antiparasite drugs might be advantageous for patients.

Transcriptome profiling in infectious diseases is a powerful approach to identify the pathways that lead to disease and can be used to understand the pathogenesis of CL (9, 28, 36). Here, we hypothesized that the genes whose expression was most variable between patients would be instrumental in dictating treatment outcome. Our results confirm our previous findings that B cell–associated transcripts are highly up-regulated in the lesions of patients infected with L. braziliensis (28) and demonstrate their variable expression between patients. Why immunoglobulin transcripts are up-regulated in patient lesions remains unclear but might reflect a role in limiting immune responses because IL-10 has been associated with parasite persistence (37) and immunoglobulin G–opsonized parasites induce IL-10 in macrophages (38).

We next looked for correlations between highly variable genes and clinical outcome, and by analyzing two independent sets of patient data, we were able to identify eight ViTALs that correlated with treatment failure in both cohorts. These included PRF1, GNLY, and GZMB, directly linking these three cytolytic genes with treatment outcome. We also identified two genes (UNC13A and APOBEC3A) that have not previously been associated with leishmaniasis but are related to genes that may influence cytolysis. For example, UNC13D regulates cytolytic granule secretion (32), whereas APOBEC3G improves cytolytic T cell recognition of target cells (31). Another gene, KIR2DL4, is expressed on CD8+ T cells and NK cells; thus, six of eight of the genes that correlated with treatment failure appear to be related to cytolysis. In addition, we found that CCL4 was the only chemokine that correlated with treatment failure. Because CCL4 binds to CCR5, which is expressed on CD8+ T cells, expression of this chemokine gene may also relate to CD8+ T cell recruitment to the lesions. Last, we also found that IFNG correlated with treatment failure, highlighting that treatment failure is not likely to be due to the absence of a protective immune response.

To determine whether the number of parasites influenced variability in gene expression, we quantitated parasites per 4-mm biopsy by qPCR and compared these values to RNA-seq reads mapped to the L. braziliensis transcriptome. Using ImmQuant to infer relative cell composition from our RNA-seq data, we found that lesions with a higher number of parasites also exhibited an increased abundance of CD8+ effector T cells and NK cells, along with higher expression of GZMB, PRF1, GNLY, and IL1B. Correspondingly, before treatment, the parasite load in the lesions from patients for whom treatment failed was higher than in patients for whom treatment was successful. Differences in parasite number could be due to differing host responses, although it did not appear to be due to a lack of IFNG gene expression. Alternatively, the difference in parasite number could also be due to variation in parasite strains or the inability of the drug to eliminate the larger parasite burden. However, the parasite burden varied relatively little across the lesion samples examined in this study, and we favor the alternative possibility that small differences in parasite burden maintain a positive inflammatory feedback loop (fig. S8). Alternatively, the inflammation promoted by increased cytotoxicity may lead to recruitment of more cells that can be infected, thus leading to parasite accumulation in the lesion. However, in previous experimental models where we observed CD8+ T cell–mediated cytotoxicity and increased disease, depletion of CD8+ T cells or blocking the inflammasome pathway limited disease without changing the parasite burden (21).

There are important limitations to this study, one of which was the inability to monitor gene expression or parasite burden over time. Such information would allow us to determine whether there were differences in the ability of pentavalent antimony treatment to kill parasites between patients and whether a reduction in cytolytic gene and IL-1β expression during treatment was associated with treatment success. Another limitation is that these studies only considered L. braziliensis patients, and in patients infected with different Leishmania species, cytolysis may be less important in promoting pathology. Last, although we have associated treatment failure with expression of several genes linked with cytolysis, our results do not show that cytolysis or the production of IL-1β actually inhibits treatment success. Studies from our experimental models indicate that this is likely, but only a clinical trial blocking IL-1β or the pathway leading to IL-1β will determine whether this is what is happening in patients.

In addition to identifying pathways leading to disease in patients infected with L. braziliensis, this study identified potential markers of treatment failure, which could have a significant impact on patient treatment. Pentavalent antimony is the drug of choice in Brazil for CL and is administered by daily intravenous injection for 20 days. In addition to having toxic side effects, this drug often fails to resolve disease, and those patients who fail the first round of treatment undergo a second or, if needed, a third round of treatment (46). Prospective identification of patients at high risk of failing antimony treatment would save patients from treatment with an ineffective drug and would expedite implementation of alternative treatments, such as liposomal amphotericin, to those who would benefit most from them. It remains to be determined whether expression of certain ViTAL genes will be useful in predicting treatment outcome for other anti-Leishmania drugs. However, because our prediction is that parasites promote increased cytolysis and inflammation, in turn leading to greater pathology, we predict that drugs that are more efficient at eliminating parasites might break the cycle of inflammation.

By exploiting RNA-seq data, we have substantially expanded our understanding of the pathogenesis of CL. We found that genes that vary in expression between patients identify pathways that can influence treatment outcome. These results directly demonstrate the significance of cytolysis in patients infected with L. braziliensis and link treatment outcome to higher expression of genes associated with cytolysis. These results also suggest that experimental murine infections that mimic the CD8+ T cell–mediated cytolysis seen in patients may be better models for testing drugs than infection of C57BL/6 mice with Leishmania (7, 8, 21). In addition to CL, there are several infectious diseases where CD8+ cytolytic T cells have been implicated in promoting increased pathology—including cerebral malaria (39, 40), trypanosomiasis (41), and coxsackievirus (42)—raising the possibility that variability in the expression of genes associated with cytolysis might predict disease progression in these diseases as well. Furthermore, we have identified several genes that had not been previously known to influence leishmaniasis, and studies are ongoing to evaluate their role in disease. Last, this study identifies parasite load as a driving factor in cytolytic gene expression and demonstrates that parasite load and cytolytic gene expression might have a potential as markers of anti-Leishmania drug treatment failure. Identification of markers of disease progression in leishmaniasis could facilitate a substantial improvement in patient care and should spur on future studies to identify new markers of disease progression in other clinical forms of leishmaniasis, a critical need for this widespread but neglected tropical disease.


Study design

Our goal was to determine whether analysis of CL lesions before treatment with pentavalent antimony would implicate genes associated with treatment failure that would define mechanisms of pathogenesis and identify potential markers of treatment failure. A characteristic skin lesion positive PCR for L. braziliensis and a positive intradermal skin test with Leishmania antigen were used to identify patients with CL. Exclusion criteria included previous anti-leishmanial treatment, individuals under 18 years old, pregnancy, or the presence of other comorbidities. Before treatment, a 4-mm punch biopsy was collected from patients with CL and endemic controls (table S1) in Corte de Pedra, Bahia, Brazil and stored in RNAlater (Thermo Fisher Scientific), or cells were collected from lesions for flow cytometry. Patients were given standard-of-care treatment (daily intravenous injections of pentavalent antimony; 20 mg/kg per day for 20 days). At day 90 after the start of treatment, patients were evaluated for lesion resolution. Cure was defined as re-epithelialization of lesions and the resolution of inflamed borders. Patients with active lesions at 90 days were defined as failing treatment and were given an additional round of chemotherapy. This study was conducted according to the principles specified in the Declaration of Helsinki and under local ethical guidelines (Ethical Committee of the Maternidade Climerio de Oliveira, Salvador, Bahia, Brazil and the University of Pennsylvania Institutional Review Board). Informed consent for the collection of samples and subsequent analysis was obtained.

Transcriptional profiling by RNA-seq

Lesions were homogenized with an MP tissue homogenizer (MP Biomedicals), and RNA was extracted using the RNeasy Plus Mini Kit (QIAGEN) according to the manufacturer’s instructions and used to prepare sequence-ready libraries using the Illumina TruSeq Total Transcriptome kit with Ribo-Zero Gold rRNA depletion (Illumina). Quality assessment and quantification of RNA preparations and libraries were carried out using an Agilent 4200 TapeStation and Qubit 3, respectively. Samples were sequenced on an Illumina NextSeq 500 to produce 75–base pair single end reads with a mean sequencing depth of 45 million reads per sample for healthy controls and 56 million reads per sample for patients with CL. Raw reads from this study, or from a previously published study (28), were mapped to the human reference transcriptome (Ensembl; Homo sapiens version 86) using Kallisto version 0.44.0 (43). For identification of parasite transcripts, raw data were filtered to remove human reads using KneadData (44), and remaining reads were mapped to the L. braziliensis reference transcriptome, MHOM/BR/75/M2904 version 2 (Ensembl) (fig. S7). Raw sequence data are available on the Gene Expression Omnibus (GEO; accession no. GSE127831).

All subsequent analyses were carried out using the statistical computing environment R version 3.5.1 (45) in RStudio version 1.1.456 and Bioconductor version 3.8 (46). Briefly, transcript quantification data were summarized to genes using the tximport package (47) and normalized using the trimmed mean of M values (TMM) method in edgeR (48). Genes with <1 CPM in n + 1 of the samples, where n is the size of the smallest group of replicates, were filtered out. Normalized filtered data were variance-stabilized using the voom function in limma (49), and differentially expressed genes were identified with linear modeling using limma (FDR ≤ 0.01; absolute logFC ≥ 1) after correcting for multiple testing using Benjamini-Hochberg. GO analysis was carried out using DAVID (50) and biological process terms. GSEA was carried out using GSEA software version 3.0 and the “C2” canonical pathways collection from MsigDB (50, 51). ImmQuant (34), together with the differentiation map (DMAP) cell phenotyping database (52), was used to predict the relative abundance of different cell types from RNA-seq data. CombiROC software was used to compute the best sensitivity and specificity for different gene combinations and plot the results in as ROC (33). Parasite load ROC analysis was performed using GraphPad Prism v7. The full RNA-seq data analysis pipeline is documented in data file S2 and is available as a fully reproducible dockerized code “capsule” archived on Code Ocean (

Quantification of parasite burden in CL lesion biopsies by qPCR

The same lesion RNA preparations used for RNA-seq were used to quantify parasite burden by qPCR. A standard curve was prepared from total RNA extracted from 107 L. braziliensis promastigotes recovered from axenic culture using the RNeasy Plus Mini Kit (QIAGEN), and complementary DNA (cDNA) was generated with the High-Capacity RNA-to-cDNA Kit (Applied Biosciences). qPCR was carried out on a ViiA 7 machine (Applied Biosciences) using Power SYBR Green Master Mix (Applied Biosciences) and primers targeting the L. braziliensis 18S ribosomal subunit (forward, 5′-TACTGGGGCGTCAGAG-3′ and reverse, 5′-GGGTGTCATCGTTTGC-3′) (53, 54) and human GAPDH (housekeeping gene: forward, 5′-GGTGTGAACCATGAGAAGTATGA-3′ and reverse, 5′-GAGTCCTTCCACGATACCAAAG-3′) (fig. S4). All reactions were carried out in duplicate, and data are presented as parasites per 4-mm biopsy.

Flow cytometry

An independent set of 13 biopsies were collected from patients with CL before treatment. Biopsies were treated with Liberase (Roche) for 90 min at 37°C and 5% CO2, dissociated, and passed through a 50-μm Medicon filter (BD Pharmingen). Cells were surface-stained and intracellularly stained with anti-CD8b PeCy5.5 (eBioscience) and anti-GzmB allophycocyanin (Invitrogen), with data collected on a FACSCanto (BD) and analyzed using FlowJo software (Tree Star).

Statistical analysis

cV was calculated using FinCal v0.6.3 (55). Differential gene expression analysis was performed with limma package (49). Permutational multivariate analysis of variance (PERMANOVA) was used to compare distances between gene expression profiles of CL lesions and healthy skin. ROC curves were made with plotROC package to calculate area under the curve for potential markers of treatment outcome (56). Mann-Whitney U tests and Spearman’s rho (ρ) correlations were performed in the R statistical environment or GraphPrism v7, and P < 0.05 was considered statistically significant.

Code Ocean

All code used to analyze RNA-seq data and clinical metadata from this study is available as a reproducible dockerized code capsule on Code Ocean.

Code Ocean link


Fig. S1. Top 100 up-regulated genes in CL lesions relative to healthy skin.

Fig. S2. Age, sex, lesion size, and delayed-type hypersensitivity response are not associated with treatment failure.

Fig. S3. CombiROC analysis identifies gene combinations for accurately predicting treatment outcome.

Fig. S4. Absolute quantification of L. braziliensis parasites in CL lesions by qPCR.

Fig. S5. Relative abundance of parasite transcripts correlates with expression of eight ViTALs.

Fig. S6. Non-ViTAL genes do not correlate with parasite load.

Fig. S7. Coverage plot for L. braziliensis read mapping.

Fig. S8. Model showing proposed role of parasite load in determining treatment outcome in CL.

Table S1. Demographic and clinical metadata from patients with CL.

Data file S1. GO enrichment analysis results for 250 ViTALs.

Data file S2. Supplementary code file.


Acknowledgments: We thank L. Guimarães and E. Lago for assistance with patient screening and sample collection. Funding: This work was supported by U01 AI088650 [International Collaborations in Infectious Disease Research (ICIDR)] to P.S., P50 AI030639 [Tropical Medicine Research Center (TMRC)] to E.M.C., and the Center for Host-Microbial Interactions (CHMI). Author contributions: C.F.A., D.P.B., and P.S. conceptualized the study. C.F.A., B.T.N., A.M.M., and D.P.B. participated in the methods and validated the results. C.F.A. analyzed the results. C.F.A., F.O.N., B.T.N., A.M.M., L.P.C., E.M.C., D.P.B., and P.S. participated in the investigation. D.P.B., L.P.C., E.M.C., and P.S. provided resources and acquired funding. C.F.A., D.P.B., and P.S. participated in data curation. C.F.A., D.P.B., and P.S. contributed in the writing of original draft, review, and editing. C.F.A., D.P.B., and P.S. supervised this study. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data associated with this study are present in the paper or the Supplementary Materials. Raw sequence data are available at the Sequence Read Archive (SRA; project no. PRJNA525604), and tabular count data are available on the GEO (accession no. GSE127831). All code used to analyze RNA-seq data and clinical metadata from this study have been archived on Zenodo (DOI: 10.5281/zenodo.3374884) and is available as a reproducible dockerized code “capsule” on Code Ocean. ( This repository contains quality control analysis of raw data and an RMarkdown document that integrates code with outputs, which was used to generate data file S2.
View Abstract

Navigate This Article