A new paradigm for leprosy diagnosis based on host gene expression

Abstract

Transcriptional profiling is a powerful tool to investigate and detect human diseases. In this study, we used bulk RNA-sequencing (RNA-Seq) to compare the transcriptomes in skin lesions of leprosy patients or controls affected by other dermal conditions such as granuloma annulare, a confounder for paucibacillary leprosy. We identified five genes capable of accurately distinguishing multibacillary and paucibacillary leprosy from other skin conditions. Indoleamine 2,3-dioxygenase 1 (IDO1) expression alone was highly discriminatory, followed by TLR10, BLK, CD38, and SLAMF7, whereas the HS3ST2 and CD40LG mRNA separated multi- and paucibacillary leprosy. Finally, from the main differentially expressed genes (DEG) and enriched pathways, we conclude that paucibacillary disease is characterized by epithelioid transformation and granuloma formation, with an exacerbated cellular immune response, while multibacillary leprosy features epithelial-mesenchymal transition with phagocytic and lipid biogenesis patterns in the skin. These findings will help catalyze the development of better diagnostic tools and potential host-based therapeutic interventions. Finally, our data may help elucidate host-pathogen interplay driving disease clinical manifestations.

Fig 1. Differentially expressed genes from RNA-seq in leprosy vs. non-leprosy.

(A) Volcano plot depicting DEG from leprosy vs. non-leprosy, where violet dashed line marks |log₂FC| = 1. For clarity, gene symbols are shown only for the largest log₂FC. (B) Heatmap with hierarchical clustering of samples based on expression of the DEG from leprosy vs. non-leprosy comparison. Color scale ranges from lower expression (blue) to higher expression (red). Library size is given in millions. LIB, logarithmic index of bacilli. (C) Biological processes from GO enriched for up-regulated DEG from leprosy vs. non-leprosy comparison. FDR, false discovery rate; NL, non-leprosy; GA, granuloma annulare; non-leprosy: GA + healthy individuals.

Fig 2. Technical and biological validation for selected DEG discovered from RNA sequencing.

(A) Tukey boxplots with RT-qPCR normalized (2–3 reference genes) log₂ expression values (A.U) according to clinical and histopathological diagnosis. ODD samples are colored according to M. leprae 16S rRNA qPCR status as positive (blue) or negative (green). (B) log₂FC from MB-ODD and PB-ODD comparisons estimated from Bayesian linear mixed models and their 95% credible intervals. (C) Tukey boxplot highlighting IDO1 RT-qPCR normalized log₂ expression values by final diagnosis grouped into ODD category. Missing values are omitted.

Fig 3. Hierarchical clustering of RT-qPCR replicated DEG and ROC analysis.

(A) Hierarchical clustering with scaled and centered normalized log₂ RT-qPCR expression values (arbitrary units) and annotated according to group and specific diagnosis. Dendrogram tree was cut arbitrarily and cluster analysis is for hypothesis generating purposes only. Two samples had more than 13 missing expression values and were removed from A. (B) Principal component analysis (PCA) with 15 genes measured by RT-qPCR and using log₂ normalized scaled data. For PCA only, missing values were imputed by the gene arithmetic mean. NA, not amplified, i.e., Cp > 40. In this regard, there were two outliers (psoriasis and erythema), which are samples with high numbers of NA values and that were imputed using the gene arithmetic mean. (C) Receiver operating characteristic analysis for genes with largest AUC (97% confidence intervals) from RT-qPCR replication samples (complete data are shown in S6 Table). See also S1 Appendix.

Fig 4. Gene candidates identified with the penalized logistic regression (LASSO) model as the most important to distinguish PB and MB leprosy lesions.

(A) Coefficients (log odds) from the top 10 most selected genes (i.e., non-zero) across 10,000 bootstrap samples using the microarray from Belone et al. as training dataset. (B) Frequency of non-zero coefficients across all bootstrap samples. (C) Misclassification error distribution estimated from 4-fold cross-validation (k-) across 10,000 bootstrap samples, with median error of 3.70% (±5.4% median absolute deviation). (D) Number of genes kept across all resamples. Predicted probability from the final model performance on this study test RNA-seq (E) and Montoya et al. RNA-seq (F). Normalized log₂ gene expression (z-score) of the two most frequently selected variables for distinguishing MB from PB samples in the (G) microarray training dataset and (H) this study test RNA-seq. PB, paucibacillary leprosy; MB, multibacillary leprosy. Tukey box plots with 1st, 2nd and 3rd quartiles ± 1.5 × inter quartile range (IQR) whiskers. See also S1 Fig.

Fig 5. Differentially expressed genes from multibacillary (MB) vs. paucibacillary (PB) leprosy lesions.

(A) Volcano plot showing DEG from the MB vs. PB comparison, where blue points are DE with |log₂FC| ≥1 and FDR < 0.1. (B) Scatter plots with the 161 DEG common between this study and Belone et al. (24) microarray for the same comparison. Red and green dashed lines indicate log₂FC of -1 and 1, respectively. Blue points are genes with the same modulation signal and red indicates discordancy. Rho, Spearman’s rank correlation coefficient. CCC, Lin’s concordance correlation coefficient. Venn diagram on the right displays the number of DEG in each study according to FDR < 0.01. (C) Biological processes from GO enriched from up-regulated and (D) down-regulated DEG. FDR, false discovery rate.

Fig 6. Strongest correlations between keratinocyte and EMT-related genes in leprosy lesions.

(A) Heat plot with Spearman’s rho correlation coefficient of the strongest correlations after multiple testing adjustment with at least one gene-pair passing FDR ≤ 0.0001 and rho ≤ -0.8. Correlations with FDR > 0.1 are filled with white. Row colored squares identify gene annotations. Scatter plots of average log₂ expression calculated with keratinocyte/epidermal development-related genes previously documented as down-regulated in leprosy skin against dedifferentiation-related genes using either (B) this study RNA-seq dataset or (C) Belone et al. microarray (GSE74481). Lines were drawn based on intercept and beta parameters estimated from robust linear regression for all samples (black line) or separately for PB (blue line), and MB (red line). Spearman’s rho coefficient along with 95% nominal confidence intervals are shown inside scatter plots calculated from all samples. See also S2 and S3 Figs.

Fig 7. Hypothetical hourglass model contextualizing the observed findings for leprosy clinical outcomes.

The host-pathogen interaction in the skin leads to opposing leprosy clinical forms. Upon infection, M. leprae induces baseline metabolic alterations such as an increase in glucose uptake, modulation of lipid biosynthesis, reduction of mitochondrial metabolism, and upregulation of IDO-1 and type I IFN. Eventually, progression towards an unspecified inflammatory state can be observed where three ways could be anticipated: I) self-healing; II) progression towards the tuberculoid pole; or III) progression to lepromatous pole. These outcomes are driven by specific environmental and host genetic factors. It is expected that lower (or shorter) M. leprae exposure, food shortage, BCG vaccination, and polymorphisms in genes controlling autophagy/granuloma formation (NOD2, LRRK2, PRKN) all contribute to developing leprosy per se. Excessive inflammation is one phenotype observed, that is also seen in other granulomatous diseases (e.g., cutaneous sarcoidosis, granuloma annulare), especially in paucibacillary lesions. On the other pole, epithelial-mesenchymal transition and local immunosuppression are present due to a probably higher (and/or longer) M. leprae exposure, combined with host single-nucleotide polymorphisms (SNPs) at key genes, like lipid biogenesis (APOE) and central metabolism (HIF1A, LACC1/FAMIN), culminating in disease progression.