Ethics statement. This study (IRB code #1273) was approved and monitored by the Institutional Review Board for ethical conduct of research involving human subjects of the Centro Internacional de Entrenamiento e Investigaciones Médicas (CIDEIM) in accordance with national (Resolution 008430, República de Colombia, Ministry of Health, 1993) and international (Declaration of Helsinki and amendments, World Medical Association, Fortaleza, Brazil, October 2013) guidelines. All individuals voluntarily participated in the study and written informed consent was obtained for each participant.
Study design and subjects. This study was designed to identify host biomarkers and innate immune functions that participate in the response to antileishmanial treatment in CL patients. Transcriptional profiling of innate immune cells and lesion biopsies was conducted. Adult patients (18 to 60 years of age) with parasitological diagnosis of active CL with a time of evolution < 6 months, and without apparent immune deficiencies (negative HIV test, no evidence of immunological disorder or treatment with medication having immunomodulating effects), who received standard-of-care treatment with Glucantime (GLUC, 20 mg/kg/day for 20 days) or miltefosine (MLF, 1.8–2.5 mg/kg daily dose for 28 days) were included in this study. Treatment outcome was evaluated at week 13 following initiation of treatment for GLUC and at week 26 for MLF. Cure was defined as complete re-epithelization and absence of inflammatory signs for all lesions. TF was defined as incomplete re-epithelization and/or the presence of induration, raised borders, or redness in any lesion; reactivation of the original lesion(s); or the appearance of new lesions.
Skin lesions biopsies samples. Skin lesion punch biopsies were obtained before initiation of treatment. Biopsy punches of 3 mm were obtained under local anesthetic, taking into account the following ratio: 1/3 of healthy skin and 2/3 of the edge of the lesion (the indurated edge, which does not include necrotic tissue). Skin biopsies were immediately stored in 1 mL Allprotect® (Qiagen). Samples were equilibrated overnight at 4°C and then stored at -20°C until processing 47.
Isolation of monocytes, neutrophils and eosinophils from peripheral blood samples. Ninety mL of whole blood anticoagulated with EDTA was obtained from each patient. PBMCs and polymorphonuclear leukocytes (PMNs) were isolated by centrifugation over a Polymorphprep™ (Axis-Shield) gradient according to the manufacturer's instructions. CD14 + Mo were purified from PBMCs using the CD14 microbeads ultrapure kit (Milteny Biotec) coupled to magnetic cell sorting (MACS). CD16 + Nφ were purified from PMNs using CD16 microbeads human kit (Milteny Biotec) and Eφ were obtained by negative selection using the eosinophil isolation kit human (Milteny Biotec). Cells were washed with cold PBS, precipitated by centrifugation, and the cell pellet was resuspended in 100–200 µL of RNAlater (Ambion) and stored at -80°C for later use. Purity of isolated cell populations was evaluated by flow cytometry and light microscopy. Samples were used only when purity was > 98%.
Leishmania strains, typing and drug susceptibility testing. Leishmania isolates were obtained from all patients and propagated in Senekjie´s biphasic blood agar and immediately stored in liquid nitrogen until use. Strains were typed by immunoreactivity to monoclonal antibodies as previously described 48. Drug susceptibility of intracellular amastigotes was estimated by evaluation of % parasite survival in PMA-differentiated U-937 cells after exposure to pentavalent antimony (SbV) at a final concentration of 32 µg/mL, compared to control without drug exposure. Leishmania strains were defined as Sb-resistant when percent reduction of the parasite burden after drug exposure was < 78%. Susceptibility cutoff was defined based on a panel of well-characterized clinical isolates presenting with intrinsic resistance or susceptibility to SbV 49,50.
RNA isolation and cDNA synthesis. RNA isolation from purified cell populations (Mo, Nφ and Eφ) stored in RNAlater was performed using TRIzol™ (Invitrogen, USA), followed by RNA cleanup with RNeasy Mini Kit columns (Qiagen, USA). Isopropanol/water (1:1) was used for RNA precipitation. RNA isolation from lesion biopsies was performed by tissue disruption, homogenization and extraction using TRIzol™ reagent as previously described 47. RNA integrity was assessed using an Agilent 2100 bioanalyzer. For RNA-seq, poly(A)-enriched cDNA libraries were generated using the Illumina TruSeq v2 sample preparation kit (San Diego, CA) and checked for quality and quantity using bioanalyzer and quantitative PCR.
RNA-seq data generation, preprocessing, and quality trimming. Single and paired-end reads were obtained on an Illumina NovaSeq 6000 at the Genetic Resources Core Facility, Johns Hopkins Department of Genetic Medicine, Baltimore, MD; or on an Illumina HiSeq1000 at the Brain & Behavior Institute - Advanced Genomic Technologies Core (BBI-AGTC) at the University of Maryland, College Park, MD. Trimmomatic 51 was used to remove Illumina adapter sequences, discard reads shorter than 40 nucleotides, and trim any 4 nucleotide rolling window with a mean Phred quality score less than or equal to 20. Sequence quality metrics were assessed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Raw data and project information available via NCBI-dbGaP study ID # 38338; phs003545.v1
Mapping cDNA fragments and abundance estimation. Reads were aligned against the human (hg38 revision 100), L.V. panamensis (TriTrypDb release 36), and L.V. braziliensis (release 26) genomes with HISAT2 (2.1.0) 52 using the default parameters. The resulting accepted hits and mapped reads were sorted and indexed via SAMtools 53 and passed to HTSeq 54 for generating count tables.
Global data assessment, visualization and differential expression analysis. Biological replicates and batch effects were assessed and visualized using an R package, hpgltools (https://github.com/elsayed-lab/hpgltools), developed and maintained in the El-Sayed lab. Normalized data were visualized using log2 transformed counts per million (CPM) reads following filtering to remove low counts (defined as any gene with a sum less than twice the number of samples or when all samples had fewer than 2 counts). Samples in which fewer than 11,000 genes were observed in a non-zero genes plot, were removed. Following data filtration, visualizations were performed to observe the sample relationships; these included density plots, boxplots of depth, coefficient of variance, hierarchical clustering analyses based on Pearson’s correlation coefficient and Euclidean distance, variance partition analyses, and principal component analyses (PCA) before and after normalization. Several combinations of normalization and batch adjustment strategies were evaluated along with surrogate variable estimation via SVA 55. Samples were queried via cure/TF status, visit number, clinic, and cell type in order to calculate the surrogate variable (SV) loadings, and the F-statistic was calculated for each variable with respect to each SV.
Differential expression analyses were performed using a single pipeline which performed all pairwise comparisons using the Bioconductor packages: limma 56, edgeR 57, DESeq2 58, EBSeq 59, and an explicitly basic analysis (using only Log2CPM values). In each case (except EBSeq and the basic analysis), the surrogate variable estimates provided by SVA were used to adjust the statistical model in an attempt to address the batch/surrogate effects. Each contrast was evaluated in context of its agreement with other methods, but the interpretations were primarily informed by the DESeq2 results. Detailed information on the analytical pipeline and scripts is available in Supplementary Material. Genes with significant changes in abundance (|Log2 fold change| ≥ 1 and false discovery rate adjusted P values ≤ 0.05) were passed to gProfileR 60 and GSVA 61. Gene ontology analyses were supplemented with manual data curation. Network analyses were performed with STRING 11.0 62. Simultaneously, gene set variation analysis (GSVA) was performed to produce an enrichment score against the mSigDB 63 datasets (C2, C7, and H) on a per-sample basis. These scores were passed to limma to evaluate the difference in GSVA score distributions for each gene set in the samples. Results from limma were then filtered according to log2 fold change, adjusted P value, and maximum GSVA score mean.
Detection of virus sequences. Kraken 2, with a supplemented version of its viral database 64, was used to check each sample specifically for the Leishmania RNA virus (LRV), as well as any putative viral reads. Confirmation of LRV absence was performed by qRT-PCR as described 36.
Weighted gene co-expression network analysis (WGCNA). WGCNA co-expression networks were generated and examined using low-count filtered, SVA-adjusted, RPKM (Reads per kilobase per million) by average CDS length, log2-transformed counts as input. Pairwise Pearson correlations between each gene pair were calculated and transformed into a signed adjacency matrix using the minimum power that resulted in a scale-free R2 fit of 0.8. The resulting modules and associated eigengenes were produced via the default correlation matrix blockwise module detection methods from WGCNA. Module scores by sample were queried against the metadata factors (cure and TF) via Pearson's correlations and scored using the P value metrics provided by WGCNA. The eigengenes were extracted from modules with scores deemed significant, manually examined, and passed to gene set enrichment methods. Modules significantly associated with the outcome and containing more than 1,000 genes were excluded from the analyses.
Composite scores. Pre-Tx RPKM values were used to construct raw and normalized (Z-score) composite scores. Selection of genes that constitute each individual composite was based on the DE genes between cures and TF that were common among Nφ, Eφ, and Mo (signature genes). Raw scores were calculated per patient, as the sum of RPKM values of the signature genes in each cell type (Eq. 1). A normalized score was also computed based on the sum of normalized RPKM values (Eq. 2). For raw and normalized scores, “I” indexes genes and “k” patients.
For each cell type:
equation 1)
$${\left(raw score\right)}_{k}={\sum }_{i=1}^{n}{RPKM}_{ik}$$
equation 2)
(\({normalized score)}_{k}={\sum }_{i=1}^{n}\frac{{RPKM}_{ik}-{\overline{x}}_{i}}{{\sigma }_{\text{i}}}\)
where
\({\overline{x}}_{i}=\frac{1}{m}{\sum }_{k=1}^{m}{RPKM}_{ik}\) and
$${\sigma }_{i}=\sqrt{\frac{1}{m} {\sum }_{k=1}^{m}{{(RPKM}_{ik}- {\overline{x}}_{i}) }^{2}}$$
Receiver operating characteristic (ROC) curves were used to explore the predictive potential of raw scores and Z-scores to discriminate therapeutic cure and failure. Sensitivity and specificity parameters were calculated for scores. Youden’s J statistic was used to define cutoff values.
Machine learning models. A series of machine learning models were generated and examined via the caret R package 19 to create transcriptome-informed classifiers of patients likely to cure or fail treatment. Initial expression sets were selected to include all data from innate cells from patients recruited in Cali and Tumaco, or data only from Tumaco patients. The starting data were log2-transformed CPM values, normalized, filtered to exclude genes with CV < 0.1, centered, filtered to exclude genes with correlations ≥ 0.95. The most variable 3,000 genes were then selected. The remaining data was split into training (0.4) and testing (0.6) sets 10 times. The training datasets were used to create k-nearest neighbor, random forest, GLM, and gradient boost models with an arbitrarily chosen mix of bootstrap and CV sampling. The test partitions were evaluated for accuracy and sensitivity/specificity with respect to the known outcome of each patient.
Statistical analyses. For the exploration and description of the sociodemographic and clinical variables, univariate analyses were performed. Categorical variables were described with frequencies and percentages. Quantitative variables were described as means (± SD) or medians (IQR) according to the distribution of the data. For the comparison of qualitative variables, Fisher's exact test or the chi2 test was used according to the distribution of the data. Fisher´s test was also used for comparison of tables larger than 2x2 65. Quantitative variables were compared using t-test or U-Mann-Whitney tests. Normality was determined with qq plots and the Shapiro wilk test. In all analyses P-values < 0.05 were considered significant. Statistical analysis was performed using GraphPad Prism version 9 and R version 4.1.3.