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

  • 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).

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

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

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

Supplementary Materials

  • stm.sciencemag.org/cgi/content/full/11/519/eaax4204/DC1

    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.

  • The PDF file includes:

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

    [Download PDF]

    Other Supplementary Material for this manuscript includes the following:

    • Data file S1 (Microsoft Excel format). GO enrichment analysis results for 250 ViTALs.
    • Data file S2 (.pdf file). Supplementary code file.

Stay Connected to Science Translational Medicine

Navigate This Article