Host transcriptomic signature as alternative test-of-cure in visceral leishmaniasis patients co-infected with HIV.

BACKGROUND: Visceral leishmaniasis (VL) treatment in HIV patients very often fails and is followed by high relapse and case-fatality rates. Hence, treatment efficacy assessment is imperative but based on invasive organ aspiration for parasite detection. In the search of a less-invasive alternative and because the host immune response is pivotal for treatment outcome in immunocompromised VL patients, we studied changes in the whole blood transcriptional profile of VL-HIV patients during treatment.
METHODS: Embedded in a clinical trial in Northwest Ethiopia, RNA-Seq was performed on whole blood samples of 28 VL-HIV patients before and after completion of a 29-day treatment regimen of AmBisome or AmBisome/miltefosine. Pathway analyses were combined with a machine learning approach to establish a clinically-useful 4-gene set.
FINDINGS: Distinct signatures of differentially expressed genes between D0 and D29 were identified for patients who failed treatment and were successfully treated. Pathway analyses in the latter highlighted a downregulation of genes associated with host cellular activity and immunity, and upregulation of antimicrobial peptide activity in phagolysosomes. No signs of disease remission nor pathway enrichment were observed in treatment failure patients. Next, we identified a 4-gene pre-post signature (PRSS33, IL10, SLFN14, HRH4) that could accurately discriminate treatment outcome at end of treatment (D29), displaying an average area-under-the-ROC-curve of 0.95 (CI: 0.75-1.00).
INTERPRETATION: A simple blood-based signature thus holds significant promise to facilitate treatment efficacy monitoring and provide an alternative test-of-cure to guide patient management in VL-HIV patients. FUNDING: Project funding was provided by the AfricoLeish project, supported by the European Union Seventh Framework Programme (EU FP7).

Evidence before this study

Human immunodeficiency virus (HIV) infection has been identified as a significant challenge facing visceral leishmaniasis (VL) control. VL treatment in HIV patients very often fails and results in extended treatment time followed by high relapse and case-fatality rates. In addition to limited and toxic treatment options, the treatment outcome assessment thus becomes imperative in VL-HIV patients to guide decisions on treatment extension, treatment adjustment or secondary prophylaxis initiation. However, with no alternatives to assess treatment efficacy, repeated invasive and painful aspiration from infected organs for microscopical detection of the parasite remains the only approach. Chronic patients thus undergo repeated tissue aspirates or empirical optimization of treatment regimens. Hence, the development of a less-invasive alternative to assess treatment efficacy represents an urgent and important unmet clinical need to achieve VL control.

Added value of this study

Our study is the first to investigate whole transcriptome changes in the severely neglected group of visceral leishmaniasis patients with a concurrent HIV infection, where the host immune response is largely uncharted but believed to be pivotal for treatment outcome. In a translational effort, the findings were biologically interpreted and systemically reduced to a 4-gene biomarker that could accurately discriminate treatment outcome.

Implications of the all the available evidence

This signature holds significant promise to facilitate treatment efficacy monitoring in R&D and could provide an alternative test-of-cure to guide patient management in VL-HIV patients.

Alt-text: Unlabelled box

Introduction

Visceral leishmaniasis (VL) is a potentially fatal, yet neglected, vector-borne disseminated infection caused by protozoans of the Leishmania donovani spp. complex [1]. Typical symptoms include fever spikes, substantial weight loss, splenomegaly and alterations of haematopoiesis. With a global estimate of 90,000 cases annually, Ethiopia together with Brazil, India, Kenya, Somalia, South Sudan and Sudan host more than 90% of all VL cases [2].

Compared to varying cure rates of around 90–95% in VL patients, treatment of patients with a concurrent Human Immunodeficiency Virus-1 (HIV) infection (referred to herein as ‘VL-HIV patients’) in endemic regions of East-Africa, Brazil and India frequently fails. Treatment failure results in extended treatments and case-fatality rates up to 25% [3,4]. This is particularly true for East Africa where antileishmanial drugs show lower efficacy rates and HIV prevalence rates of 10–20% are reported amongst VL patients [3]. Even if apparent parasitological clearance at end of treatment and viral suppression with ART is achieved, up to 60% of VL-HIV patients will develop recurrent relapse, typically within 3–6 months after initial cure (compared to 1–5% in immunocompetent VL patients) [3,5]. Hence, treatment outcome assessment is imperative in VL-HIV patients to guide decisions on treatment extension, treatment adjustment or secondary prophylaxis initiation.

To date, a repeated invasive and painful aspiration for microscopical detection of the parasite from infected organs (spleen, bone marrow or lymph nodes) remains the only approach for test-of-cure. Although spleen aspiration shows the highest sensitivity [6], it has a life-threatening risk of splenic haemorrhage that renders it unsuitable in patients with severe thrombocytopenia. In addition, these techniques require a great level of expertise, training of personnel and appropriate facilities where blood transfusion and management of intraabdominal bleeding is possible. Due to these reasons, chronic patients often undergo repeated tissue aspirates or empirical optimization of treatment regimens. Hence, the development of a less-invasive alternative to assess treatment efficacy represents an urgent and important unmet clinical need.

Molecular techniques for parasite detection seem promising, but could be less suitable as the parasitic load in blood decreases steeply after two days of treatment and gives no information on the host's immunological recovery [7]. In immunocompromised individuals in particular, host immune response restoration has been shown to be pivotal in the efficacy of VL treatment [8]. Therefore, transcriptomic signatures in peripheral blood may reflect immunological responses underpinning clearance or persistence of parasites. In recent years, gene signatures derived from blood transcriptomic profiling have shown great promise in treatment monitoring for a number of infectious diseases [9], [10], [11]. With a 5-gene signature, robust prediction of treatment failure in tuberculosis patients could be achieved after 1 or 4 weeks of therapy [11]. Similarly, Liu et al. identified and validated a 10-gene signature that predicted Ebola treatment outcome with an accuracy of 85% to 92% [9]. Yet, previous studies were often confined to single timepoint measurements and purely statistical approaches that could be complicated by patient-to-patient variation and little biological relevance of selected genes, all jeopardizing their generalizability.

To date, two small-scale studies in VL patients described distinct expression profiles in respectively the blood and the lymph nodes before and after treatment with amphotericin B [12] or sodium stibogluconate [13]. Likewise, Gardinassi et al. also defined distinct immunological signatures in the blood for active and cured L. infantum infected patients [14]. Although performed in immunocompetent VL patients, these findings support the pursuit of a blood-based test-of-cure. None of the previous studies, however, sought to translate such knowledge into a clinically useful signature to guide the challenging clinical management of VL-HIV patients and facilitate the evaluation in clinical trials of urgently needed novel therapeutics.

Here, to minimize the impact of heterogeneity in gene expression levels amongst individual VL-HIV patients we studied relative changes in the gene expression profile during VL treatment. By combining biological insight with a stringent machine learning approach, we have identified a relevant and robust 4-gene blood signature that accurately determined treatment outcome in VL-HIV patients.

Materials and methods

Study population and design

In 2014–16, a total of 30 patients recruited at the Leishmaniasis Research and Treatment Centre (Gondar, Ethiopia) in a randomized clinical trial (RCT) on combination treatment (NCT02011958, sponsored by the Drugs for Neglected Diseases Initiative) were included in this study for additional whole blood sampling [15]. In short, AmBisome monotherapy and a combination treatment of AmBisome and miltefosine were administrated. AmBisome monotherapy dosage was 40 mg/kg total dose, IV infusion of 5 mg/kg on day 1–5, 10, 17 and 24. In the combination treatment AmBisome dosage was 30 mg/kg total dose, IV infusion of 5 mg/kg on day 1,3,5,7,9 and 11 and miltefosine every day for 28 days (50 mg if patient weight < 25 kg, 100 mg if patient weight is >25 kg). Treatment success was defined as absence of parasites in spleen aspirate at D29. Treatment failure at D29 was defined as presence of parasites at the D29 assessment, or no clinical response to treatment requiring rescue medication on or before D29. Treatment was extended or rescue medication was provided after D29 based on parasitological examination and clinical improvement. Once the patients had a negative parasitology result, they started a follow-up period of one year (up to D390) to assess long-term relapse-free survival and safety. Rescue therapy was also given to all patients who had a confirmed VL relapse during follow-up (for more details, see [16]). After exclusion of 2 patients due to a technical error in the RNA isolation process, 28 VL-HIV patients were included in the herein described analyses (Fig. 1).

Fig. 1

Flow charts of study participants and study design. (a) Flow chart of inclusion and stratification by treatment outcomes of study participants (b) Lay-out of the study design and generation of 4-gene signature. LTFU: lost to follow-up, D: Day of treatment, VL-HIV: visceral leishmaniasis-HIV co-infection, DGE: differential gene expression.

Ethics approval and consent to participate

The study protocol was approved by the University of Gondar Institutional Review Board, the Ethiopian National Research Ethics Review Committee, the Médecins Sans Frontiers Ethics Review Board, the London School of Hygiene and Tropical Medicine Research Ethics Committee, the Antwerp University Hospital Ethics Committee, the Prince Leopold Institute of Tropical Medicine Institutional Review Board and for this specific study by the University of York Biology Ethics Committee. All patients provided informed consent and the study was carried out in accordance with international guidelines (Helsinki declaration, Good Clinical Practices and local regulations).

RNA isolation and sequencing

Whole blood was collected and stabilized in 2.5 mL PAXgene Blood RNA tubes (PreAnalytiX GmbH, Hombrechtikon, Switzerland) and stored at −80 °C. Isolation and purification of total RNA was performed using the PAXgene blood RNA kit (PreAnalytiX GmbH, Hombrechtikon, Switzerland) according to the manufacturer's instructions. Next, messenger RNA was enriched by depleting ribosomal and globulin RNA (Globin-zero gold removal kit, Illumina Inc, San Diego, USA). RNA concentrations were measured with NanoDrop 1000 Spectophotometer (NanoDrop Technologies, USA) and RNA integrity using an Agilent 2100 Bioanalyzer (Nano kit, Agilent, CA, US). Subsequently, cDNA amplification, adaptor ligation and indexing were carried out on 1 µg of total RNA by using TruSeq stranded mRNA library preparation (Illumina Inc, San Diego, USA). Libraries were sequenced on an Illumina NextSeq500 instrument (single-end, 75 bp) using 1.2 pM and 1.89% PhiX with a total of 4 runs and an average coverage of 19.3 million reads per sample. 98% of the trimmed passed-filter reads mapped against the human genome.

Mechanistic approach

Longitudinal and inclusive DGE analyses (FDR corrected p-value ≤ 0.05 and a ≥ 1.5 absolute fold difference) between the timepoints D0 and D29 were performed for both the treatment success and failure group, based on a paired generalized linear model with CLC Genomic Workbench software V12 (Qiagen Bioinformatics). Next, the pre-ranked feature within Gene Set Enrichment Analysis (GSEA) software v.3 (BROAD institute, California, USA) was used to determine enriched gene sets in the treatment groups. The well-defined hallmark, curated canonical and gene ontology datasets from the Molecular Signature Database (MsigDB) were used for the enrichment analyses. In addition, previously described whole blood transcription modules (BTMs) were evaluated [17,18]. Gene sets were considered significant if their false discovery rate (FDR) was less than 25%. Default parameters were used with 1000 permutations and the conservative classic enrichment statistics was used for calculation of the enrichment score (ES). Gene sets that did have not a minimum of 10 genes shared with the pre-ranked gene list were excluded from the analyses, and normalized enrichment scores were depicted. Leading edge analyses were run on the GSEA results to identify core group genes.

Non-redundant biological terms were visualized in a functionally grouped network by means of the ClueGo plug-in v2.5.4 for Cytoscape v.3.7.1 [19,20]. Enrichment/depletion analyses were carried out with a two-sided hypergeometric test for datasets from the EBI-UniProt GO, KEGG, REACTOME and WikiPathways databases. All genes were recognised and 45 (11.3%) were not functionally annotated. The applied kappa statistics for connectivity was kept at the default value of 0.4. GO term fusion was selected to reduce redundancy. In addition, only terms with a p-value of <0.100 were shown, as well as a GO tree interval level of 3 to 8 (medium to detailed network specificity). The autoannotate plug-in was used for cluster delineation and naming.

Discriminatory approach

A cross-sectional and stringent machine learning approach was designed. The count table was normalized for the ‘per sample’ sequencing depth using the DESEQ2 package in R. Genes that had a read count of zero in 80% of the samples were removed. The D29/D0 expression values were used as input for the machine learning approach to correct for patient specific expression patterns. We adopted a random forest (RF) classification with leave-one-out cross-validation. Within each cross-validation loop, a feature selection was first applied on the remaining samples, only selecting the differentially expressed genes between treatment cure and failure based on a more stringent DGE analyses (FDR corrected p-value ≤ 0.01 and a ≥ 2 absolute fold difference) to target a limited gene set of highly discriminatory genes. These DEGs were then used to construct the RF classifier. The number of trees, with a maximum depth of 4 and equal minimal number of samples in each leaf, was limited to 100. The relative importancies of the genes in the forest decision making were stored if the left-out sample was correctly predicted. After summing the classifier importancies for each gene from every successful loop, the top 20 genes were depicted. The final 4-gene RF classifier was constructed similarly, but to reduce overfitting on the smaller training sets during bootstrapping (see ‘Statiscal and visual analyses’), we applied a maximum tree depth of 3 and a minimal number of 3 samples in each leaf. All these analyses were carried out using the sklearn package in Python v3.7.0 [21].

Reverse transcription quantitative polymerase chain reaction (RT-qPCR)

Predesigned PrimeTime®qPCR assays (Integrated DNA Technologies (IDT), Iowa, USA) were used for human gene expression analysis of HRH4 (Hs.PT.58.4403324), IL10 (Hs.PT.58.2807216), PRSS33 (Hs.PT.58.3087623), SLFN14 (Hs.PT.58.993231) and the reference gene TBP (Hs.PT.58v.39858774) (File S2). The recommended assay with the highest algorithm score was selected for each gene. After read counts were normalized via DESEQ2, TBP was selected as reference gene based on lowest maximal deviation between samples, lowest coefficient of variation, lowest LOG2FC change, commercial availability and literature review. Previous studies have validated TBP as a reliable reference gene for qPCR [22,23] and it was detected in all study samples.

Samples with sufficient left-over volume to reach an equal amount of 0.5 μg in every RNA extract were included (Fig. S2). Both, cDNA synthesis with iScript cDNA synthesis Kit (BioRad laboratories, California, USA) and qPCR were performed on the QuantStudio 5 real-time PCR system (Applied Biosystems, ThermoFisher, Waltham, MA USA). The PrimeTime gene expression master mix with low passive reference dye (IDT) was used. To reduce technical variation, all samples from 1 patient were run in duplicate for the four genes and reference gene on one plate simultaneously. As the efficiencies were almost all ranging between 90–110%, we used the comparative Ct-value method (ΔΔCt) for relative quantification between the timepoints. Individual LOG2FC between D29 and D0 were calculated based on CPM values or Ct-values for the RNA sequencing or RT-qPCR results, respectively.

Statistical and visual analyses

For patient characteristics (Table 1), continuous data was represented as medians with interquartile ranges and categorical data as numbers and frequency. Statistical significance between treatment success and failure patients was determined with the non-parametric Mann-Whitney U test for continuous variables and two-sided Fisher's exact test for categorical variables. P-values < 0.05 were considered significant. The ROC curve of the final 4-gene classifier was generated using the sklearn package in Python. The AUC, sensitivity and specificity were calculated over 1000 bootstrap iterations. For each iteration, the dataset was split in four parts (sklearn, StratifiedKFold), of which three parts were used for training and one part for testing. Principal component analysis (PCA) scatter plots, heatmaps (Euclidean distance, complete linkage), venn diagrams and volcano plots were created with the CLC Genomic Workbench software V.12 (Qiagen Bioinformatics). Violin plots were created with Graphpad Prism v8.0.1 (GraphPad Software, San Diego, USA).

Table 1

Patient characteristics before and after treatment, stratified by treatment outcome.

	Total (n = 28)	Success (n = 12, 43%)	Failure (n = 16, 57%)	p-value
Demographic characteristics
Age, median (1–99%)	34 (27–45)	33 (28–45)	34 (27–44)	0.815
Men, n (%)	28 (100)	12 (100)	16 (100)	1.000
Clinical characteristics at D0
Treatment regime, n (%)				0.491
AmBisome	9 (32.1)	3 (25)	6 (37.5)
AmBisome + Miltefosine	19 (67.9)	9 (75)	10 (62.5)
Primary VL, n (%)	14 (50)	5 (41.7)	9 (56.3)	0.704
Site of parasite detection, n (%)				0.141
Spleen (+1 to +2)	1 (4)	1 (8)	0 (0)
Spleen (+3 to +6)	26 (92)	11 (92)	15 (94)
Bone marrow	1 (4)	0 (0)	1 (6)
On ART, n (%)	23 (82.1)	8 (66.7)	15 (93.8)	0.133
HIV viral load, n (%)				0.049
Undetectable	8 (29)	1 (8)	7 (44)
>20–10,000 copies/mL	7 (25)	3 (25)	4 (25)
>10,000 copies/mL	13 (46)	8 (67)	5 (31)
CD4 count, median (1–99%)	47 (11–159)	44.5 (11–159)	61.5 (15–152)	0.516
Concomitant diseases, n (%)	23 (82.1)	10 (83.3)	13 (81.3)	1.000
Platelet count (x103/μl), median (1–99%)	100 (47–195)	116 (47–181)	98 (55–195)	0.693
Clinical characteristics at D29
Site of parasite detection, n (%)				NA
Spleen (+1 to +2)	6 (21)	NA	6 (37.5)
Spleen (+3 to +6)	6 (21)	NA	6 (37,5)
Bone marrow	4 (14)	NA	4 (25)
CD4 count, median (1–99%)	155 (20–341)	153 (20–285)	171 (42–341)	0.807
Platelet count (x103/μl), median (1–99%)	221 (87–701)	206 (115–464)	235 (87–701)	0.710

VL: visceral leishmaniasis, ART; antiretroviral therapy, HIV: human immunodeficiency virus, D: day of treatment, NA: not applicable.

Results

Nested in a clinical trial, 30 VL-HIV patients underwent additional whole blood sampling. After exclusion of two patients due to technical error, a total of 28 VL-HIV patients were included in the analyses (Fig. 1a). Twelve (43%) patients reached parasitological cure after 1 cycle of treatment and 16 required extended or rescue treatment. At admission, half of the cases had a history of VL and 82% were on antiretroviral therapy (ART) (Table 1). In both treatment outcome groups, a similar CD4 reconstitution of around 110 cells/μL was observed after 29 days of VL treatment compared to baseline. Except for a higher percentage of viral controllers in the treatment failure group (possibly linked to the higher number of people on ART), no significant differences in patient characteristics nor treatment regimens could be observed between the treatment success and failure group. These findings indicate that HIV-related markers were not associated with VL treatment outcome and that despite a successful suppression of the viral load or CD4 T cell recovery, VL treatment failure still occurred and that vice versa, non-controlled HIV patients were able to clear the infection.

Distinct gene expression changes in treatment success and failure patients

We first used principal component analysis (PCA) to reduce the dimensionality of the data and determine whether patients could be clustered. However, after pooling all the samples, the first principal component only explained 11.9% of the total variance in the dataset and this was associated with a lower read count in one of the samples (Fig. S1).

To account for this high variability in transcript abundance between individual VL-HIV patients, possibly due to diverse HIV or VL histories, genetic backgrounds, nutritional status and a high presence of comorbidities (>80%, Table 1), we performed longitudinal paired analyses of intra-individual changes in gene expression levels between Day 0 (D0) and D29 (using the mechanistic approach described in Methods, Fig. 1b). Differential gene expression (DGE) analyses between D0 and D29 resulted in distinct transcriptomic patterns by treatment outcome, with 397 and 194 unique DEGs for the treatment success and failure group, respectively (Fig. 2a, File S1 for complete lists). Only 142 were common DEGs between the groups, of which only one (CTD-2545M3.6) encoding an uncharacterised protein was regulated in the opposite direction in the groups. An almost equal distribution of up- (n = 218/397, 55%) and down- (n = 179/397, 45%) regulated genes was observed during treatment success, while in the treatment failure group less genes significantly changed in expression level and their expression was mainly downregulated (n = 166/194, 86%) (Fig. 2a). This observation was also reflected in corresponding unsupervised hierarchical clustering of the expression data (Fig. 2b&c). In the treatment success group, all samples were distinctly clustered by timepoint and showed a reversed pattern at D29 (Fig. 2b). In the treatment failure group, no primary clustering by timepoint could be observed with a more scattered pattern (Fig. 2c), indicating no mutual delineated impact of VL treatment on host transcriptomic profile. Altogether, these results suggest that significant changes take place in the abundance of blood transcripts during successful VL treatment and that the patterns of gene expression depend on the treatment efficacy.

Fig. 2

Distinct patterns of differentially expressed genes (DEGs) in treatment success and failure patients. (a) Overlap of differentially expressed genes (DEGs) between D0 and D29 for treatment success and failure patients. Venn diagram showing the number, uniqueness and directionality of the DEGs (>1.5 absolute fold change with FDR p-value <0.05 – see mechanistic approach in materials and methods) between D0 and D29 for treatment success (orange, n = 12) and failure cases (red, n = 16) (b) Unsupervised clustering of differentially expressed genes (DEGs) in treatment success and (c) failure patients. Heat maps showing the z-scores (bottom scale) or extent of up (red) or downregulation (blue) of the DEGs (Y-axis) at D0 (green samples) and D29 (purple samples) for all individual treatment success patients (n = 12) and all individual treatment failure patients (n = 16). Unsupervised hierarchal clustering of the samples was applied and based on Euclidean distance and complete linkage. D: Day of treatment.

Functional annotation of unique gene transcripts underlying successful VL treatment

To better understand the ontology and function of the 397 unique genes that showed significant alterations in their expression levels during successful VL treatment, we assessed a total of 6161 published gene sets across 4 different databases for significant enrichment. In total, 77 biological GO terms, 12 canonical pathways, 8 hallmark gene sets and 5 B.M. were enriched (Fig. 3). Despite more than 200 upregulated genes in the treatment success cases, enriched gene sets across the four different databases showed a robust downregulation of pathways by D29, in particular those linked to cell cycle, apoptosis, proteolysis and adaptive immune response processes. The top 8 most influential genes that were part of the core enrichment group in more than 35 gene sets, included different alpha subunits of the proteasome alpha (PSMA6, PSMA2, PSMA3, PSMA4), one subunit of the proteasome 26 s (PASMD14), cyclin-dependant kinase 1 (CDK1), interferon-γ (IFNG) and interleukin 6 (IL6). In contrast, no significant enrichments could be observed within the 194 unique DEGs from the failure group, consolidating the observation of no clear delineated impact of VL treatment on host pathways amongst patients with treatment failure after one month of treatment.

Fig. 3

Enriched gene sets in the treatment success group (n = 12). Normalized enrichment scores for hallmarks, GO terms and canonical pathways from the MsigDB database and previously described blood transcriptional modules are depicted in blue (downregulated) or red (upregulated) scale. The size column represents the number of genes from the respective gene set found in the expression dataset. NES: normalized enrichment score.

To obtain insight in the correlations between enriched terms and define overarching pathways or processes, a functionally grouped network of enriched categories was generated that showed 2 distinct clusters (Fig. 4). The main cluster was downregulated and consisted of cell cycle features (green), including (stem) cell differentiation, cell division, antigen presentation, cell proliferation and related catabolic or proteolytic processes. The second cluster was upregulated and covered features of phagolysosomal activity to kill or stop the growth of fungal and microbial organisms (purple), including G-protein coupled receptor signalling and vacuolar activity with primary (azurophilic) granules loaded with defesins and serine proteases.

Fig. 4

Functionally grouped network analyses of enriched gene sets in the treatment success group (n = 12). Functionally related clusters of enriched gene sets are represented by different colours and the label of the most significant term per cluster is shown. The node size represents the term enrichment significance and the connectivity by grey lines (kappa statistics 0.4). Similar terms were fused to reduce redundancy. Bottom tables showed the directionality of significant terms and the frequency of genes from the respective terms detected in the expression dataset. The exact number of genes detected is indicated behind each bar.

These findings were indicative of a clear remission of active disease in successfully treated patients. The data indicated a distinct and profound downregulation of the host cellular activity and immunity and suggested clearance of the parasite by antimicrobial peptides in phagolysosomes.

Gene set to accurately distinguish treatment outcome at D29

After having established that the whole blood transcriptome markedly changed during the course of treatment and reflected relevant biological processes, we investigated whether this transcriptomic change could be used to accurately determine the treatment outcome at D29. To answer this question, a machine learning approach was adopted that compared D0/D29 transcript abundance between treatment success and failure patients (using the discriminatory approach described in the Methods, Fig. 1b). To avoid overfitting, the random forest (RF) classifiers were evaluated by leave-one-out cross-validation. In total, 70% of the constructed classifiers could correctly predict the outcome of the left-out patient. The top 20 genes with the highest summed feature importance in the tree decision making over all correctly classifying classifiers were kept and included several pseudogenes and novel transcripts (Table 2, first column). To increase the biological relevance and robustness of selected genes, these 20 genes were subsequently compared with the list of DEGs from the mechanistic approach (D0 vs D29) for overlap (Table 2). Four genes were found to be shared with the treatment success group, but none with the failure group.

Table 2

Comparison of differential gene expression approaches and identification of the 4-gene set.

Discriminatory Approach		Mechanistic Approach	Gene details
DEG (success vs failure)		DEG (D0 vs D29)
Gene symbol	Relative importance in RF classifiers	Shared gene (treatment group)	Protein	Function	Cell specificity	Associated enriched gene sets	Regulation at D29 (LOG2FC)
MAN1C1	0.158	No
PRSS33	0.138	Yes (success)	Serine protease 33	Amidolytic activity	Predominantly expressed in macrophages and peripheral leukocytes	Proteolysis; protein kinase signalling; Ubiquitin/proteasome pathway	Upregulated (3.31)
AP001830.1	0.105	No
FAM153B	0.088	No
BCAS1	0.068	No
NA	0.061	No
AC007611.1	0.055	No
HRH4	0.051	Yes (success)	Histamine receptor H4	Binds to histamine in peripheral tissue and mediates inflammation	Predominantly expressed on haematopoietic cells	Purinergic (G-protein coupled) receptor activity; inflammatory response	Upregulated (2.6)
NAV1	0.049	No
VLDLR	0.042	No
IL10	0.039	Yes (success)	Interleukin 10	Cytokine with pleiotropic effects on immunoregulation and inflammation. Downregulates Th1 cytokines and antigen presentation	Primarily produced by monocytes and lymphocytes	negative regulation of antigen presentation; negative regulation of IFN-y production; Type 2 immune response; NF-κB signalling; defence response to protozoa/bacteria	Downregulated (−2.45)
KRT7	0.033	No
AC135068.2	0.031	No
AL390198.2	0.023	No
H1F0	0.018	No
LINC00205	0.017	No
GEMIN7-AS1	0.011	No
SLFN14	0.01	Yes (success)	Schlafen family member 14	Important role in platelet formation and function (endoribonuclease)	Abundant in immune cells, located in the nucleus	Regulation of DNA/protein metabolic and catabolic processes	Upregulated (1.74)
KLHL13	0.002	No
AC093159.1	<0.001	No

DEG: differential gene expression, RF: random forest, LOG2FC: log2(fold change), Th1: T-helper 1, IFN-γ: interferon-γ, NF-κB: nuclear factor kappa-light-chain-enhancer of activated B-cells.

Next, a new RF classifier was built to assess whether the combination of these 4 genes alone was able to accurately determine the treatment outcome. Bootstrapping this classifer a 1000 times (each iteration one quarter of dataset was used for testing) showed a mean area-under-the-ROC-curve (AUC) of 0.95 (95% CI: 0.75–1.00) - with a mean sensitivity of 84% (95% CI: 83–85) and specificity of 85% (95% CI: 83–86) - to correctly distinguish treatment failure from success cases (Fig. 5). We validated the findings by reverse transcription quantitative polymerase chain reaction (RT-qPCR) and showed similar distribution patterns in log2(foldchange) (LOG2FC) changes for all 4 genes (Fig. S2 and File S2). Due to the lack of follow-up samples in the 16 treatment failure patients, we could not investigate whether these patients also expressed the 4-gene signature when reaching parasitological cure after extended treatment (Fig. 1a). We did, however, observe that three out of the four genes were also differentially expressed in the 13 VL-HIV patients with long-term cure (cf. no relapse) at 6 months after treatment end (Fig. S3a), of which eight were treatment failure cases at D29. In agreement, the four genes could not be detected in treatment failure patients at D29 with low parasite burden (1+,2+) in splenic aspirates, indicating a true cure signature that reflected disease remission after an apparent complete resolution of the parasite.

Fig. 5

Receiver operator curves (ROC) for the final random forest classifier based on the 4-gene set. (a) Showing the average AUC value of the 4-gene random forest classifier, calculated over 1000 bootstrap replicates. (b) Bar chart showing the sum of relative importancies in the classifier over 1000 bootstrap replicates.

With regard to long-term cure prediction, distinct patterns at D29 could be observed in gene expression levels of successfully treated patients who relapsed in the following year (n = 5, 42%) and those who did not (Fig. 1). Yet, our 4-gene set was mostly found in the common DEGs between long-term cured and relapsed patients, except for SLFN14, indicating little value in relapse prediction (Fig. S3b). The relapsed patients overlapped for 61% with the failure group at D29, compared to only 30% with the long-term cured patients. These findings argue for the investigation of a predictive signature at D29 to predict relapse development in the next 6–12 months, but patient numbers were too small to perform robust analyses.

Discussion

In this study, we investigated the whole blood transcriptome of 28 well-characterized Ethiopian VL-HIV patients before and after 1 month of VL treatment. We observed a distinctive pattern of disease remission in successfully cured patients and the complete lack of it in treatment failure cases. Subsequently, we identified a 4-gene signature able to discriminate treatment success at D29 with a sensitivity of 84% and specificity of 85%. Application of this signature as a low-invasive tool to assess treatment efficacy in VL-HIV patients could have significant value in guiding patient management and future R&D.

Due to a possible homoeostatic mechanism to control persistent infection-induced inflammation in active VL, elevated levels of the regulatory cytokine IL-10 that significantly decreased following successful treatment, have been reported in numerous animal and human studies [24,25]. This marked downregulation of IL10 transcription in successfully treated patients was strongly recognized in our 4-gene signature. This is, however, in contrast with the two previous transcriptomic studies in VL patients that did not find a significant regulation of IL10 [12,13]. While we cannot rule out that our signature is HIV-specific, our longitudinal approach that accounts for natural variation in gene transcription between patients and the use of RNAseq instead of microarray, could have contributed to detecting these associations. Likewise, the other three genes (HRH4, PRSS33 and SLFN14) have not been reported in the previous two studies.

In addition, we could not confirm the IFNG gene as the major regulator gene as reported in previous transcriptome studies of VL patients [12]. Although the IFNG gene and pathways were also strongly downregulated in treated VL cases compared to active cases, we rather observed a prominent upregulation of vacuolar activity with primary (azurophilic) granules loaded with defesins and serine proteases. These genes and terms suggested an enhanced phagocytosis, neutrophil involvement and complement activation; and may reflect the proteolytic degradation responsible for parasite clearance. This process was reflected in the 4-gene signature by the marked upregulation of the amidolytic activity of serine protease 33. Its predominant expression in peripheral macrophages and highest discriminatory importance in the classifier supports future studies into the role of PRSS33 in phagolyosomal degradation of Leishmania parasites.

It has been postulated that the interactions of Leishmania spp. with eosinophils and mast cells influence the macrophage response to infection and the development of an adaptive immune response [26]. An important mediator in this process is the histamine receptor H4, also designated as the ‘immune system histamine receptor’ due to its chemotactic properties (eosinophil migration and mast call chemotaxis) and involvement in dendritic cell activation and T-cell differentiation. Because histamines have been described to have an assisting role in vitro clearance of Leishmania infection [26], the observed upregulation of HRH4 in successfully treated patients could indicate clearance of parasites. Nevertheless, due to its importance in the 4-gene signature, the role of histamine, mast cells and its receptors on treatment efficacy in VL and VL-HIV patients merits further research. In particular because the H4 receptor is also being explored as a promising drug target in many chronic inflammatory disorders [27].

Defects in the SLFN14 gene of the signature are associated with thrombocytopenia and excessive bleeding [28]. Likewise, the platelet count is known to be altered during VL infection due to splenic sequestration and often leads to bleeding tendency. Hence, the upregulation in SLFN14 expression may reflect the described return of platelet counts to normal levels after successful cure. Yet, no clear correlation could be observed between SLFN14 expression and recovery of platelet count after therapy in our study (Figure S4). It is of note that the gene was absent in the DEG list at D210 in long-term cured patients (Figure S4). Therefore, its importance should be validated as it could have been influenced by a small group of patients with severe thrombocytopenia. Our top 20 list from the discriminatory approach consisted of many pseudogenes and novel transcripts that also warrant further research, as they are not yet fully understood but are increasingly acknowledged as key contributors to immune responses [29].

Unfortunately, we could not study differences in treatment regimens due to sample size restrictions. Therefore, it remains unknown whether the broad antiproliferative activities of miltefosine contributed to the strong predominance of the downregulated cell cycle cluster in the network analyses. This predominance is however believed to reflect a cessation of the massive proliferation and differentiation of T and B cell clonotypes during the acute adaptive immune response, as similar results were observed in studies of immunocompetent VL patients in India and Brazil treated with amphotericin B (inciting cell wall disruption) and pentavalent antimonials (broad macromolecule inhibitor), respectively [12,14]. This would indicate it is not HIV-specific nor affected by the administered drugs. Nevertheless, a true mechanistic evaluation of affected organs with single cell evaluations should be studied to assess the causal and time relationship with drug modulation and parasite reduction. It is also to be noted that we could not observe enriched gene sets in the treatment failure group but noticed a significant decline in expression of many IgG-related genes (almost 40% of all DEGs). Many of these genes were also shared at D29 with the seemingly cured patients who relapsed within the next year (35% of common DEGs). This finding could suggest the initial presence of hypergammaglobulinemia at diagnosis as a risk factor for treatment failure or relapse, and merits further research.

Despite the high inter-patient variation in gene levels, we were able to construct a good classifier. We believe our methodology applied some good general practices that could be incorporated in future studies. Previous studies targeting limited transcriptomic gene sets for diagnosis or treatment outcome prediction in VL or even infectious diseases in general were mostly confined to cross-sectional comparisons at a single timepoint [9,10,30]. Such studies can be complicated by patient-to-patient variation in transcript abundance and are highly dependant on the stability of the used housekeeping gene for quantification in future validation studies with RT-qPCR. Therefore, we targeted intra-individual and relative changes in transcript abundance that are likely to enhance its generalizability but require both pre and post sampling. Additionally, purely statistical approaches were often adopted to select the best discriminatory genes. Such approaches can have the inherent danger of identifying genes with little biological relevance or interpretation, jeopardizing their validation in other patient cohorts [31]. Correspondingly, our top 20 list from the discriminatory approach consisted of many pseudogenes and functionally uncharacterized transcripts. By overlaying the machine learning results with mechanistic insight analyses, we believe the 4 genes represent relevant but distinct pathways that reflect underlying biological changes during treatment. When monitoring of treatment outcome is envisioned, the specificity of these pathways or genes to Leishmania infection is also less critical, compared to signature development for diagnosis or disease development prediction.

A limitation of this study is that the 4-gene set identified here is derived from data collected from a single, rather limited cohort of 28 male VL-HIV patients in a single geographic area in East Africa. Although it will be necessary to further validate this biomarker in larger cohorts of VL-HIV patients, and to include patients from other regions/countries, this should not detract from the value of such a biomarker. Even if restricted to use in East Africa, it would still be of high interest as the burden of VL-HIV is significant in this area. It will be of future interest to evaluate performance of this biomarker as the chemotherapeutic landscape alters, though the combination of AmBisome/miltefosine is the most effective treatment of VL-HIV patients in East Africa [15] and our data suggested that the identification of the 4-gene signature was not merely a reflection of drug choice.

In the near future, a composite endpoint consisting of clinical signs, antigen detection methods and a host 4-gene biomarker could be a powerful, less-invasive tool for research and drug development (R&D) in leishmaniasis, both for immunocompetent and VL-HIV patients. In particular, future studies should assess the predictive value of the signature at earlier timepoints and evaluate the methodology required to develop a more point-of-care detection system. This may guide early case management, treatment recommendation and could facilitate early futility analyses and dose-finding studies of novel compounds. It is also to be noted that its value in similar parasitic infections could be studied, as rather general disease remission signatures were detected. Our results also encourage future investigations in VL relapse prediction by blood-based transcriptomic signatures. Yet other processes are likely to be predictive in parasite recrudescence and to be reflected in such a signature. However, the specificity of the signature to Leishmania will also become more relevant. For example, an extensively validated blood-based 3-gene signature recently satisfied the WHO criteria for a non-sputum-based triage test for tuberculosis across heterogeneous culture-confirmed datasets [30].

In conclusion, we have identified an easily accessible 4-gene blood signature with high discriminatory value to assess treatment efficacy at the end of treatment in VL-HIV patients. Application of this signature as a low-invasive tool to assess treatment efficacy in VL-HIV patients could have significant value in guiding patient management and future R&D.

Declaration of Competing Interest

The authors declare that there are no competing interests.