Mycobacterium leprae and host immune transcriptomic signatures for reactional states in leprosy

Background Mycobacterium leprae transcriptomic and human host immune gene expression signatures that demonstrate a plausible association with type I (T1R) and type II reactions (T2R) aid in early diagnosis, prevention of nerve damage and consequent demyelinating neuropathy in leprosy. The aim of the study is to identify M. leprae and host-associated gene-expression signatures that are associated with reactional states in leprosy. Methods The differentially expressed genes from the whole transcriptome of M. leprae were determined using genome-wide hybridization arrays with RNA extracted from skin biopsies of 20 T1R, 20 T2R and 20 non reactional controls (NR). Additionally, human immune gene-expressions were profiled using RT2-PCR profiler arrays and real-time qPCRs. Results The RNA quality was optimal in 16 NR, 18 T1R and 19 T2R samples. Whole transcriptome expression array of these samples revealed significant upregulation of the genes that encode integral and intrinsic membrane proteins, hydrolases and oxidoreductases. In T1R lesional skin biopsy specimens, the top 10 significantly upregulated genes are ML2064, ML1271, ML1960, ML1220, ML2498, ML1996, ML2388, ML0429, ML2030 and ML0224 in comparison to NR. In T2R, genes ML2498, ML1526, ML0394, ML1960, ML2388, ML0429, ML0281, ML1847, ML1618 and ML1271 were significantly upregulated. We noted ML2664 was significantly upregulated in T1R and repressed in T2R. Conversely, we have not noted any genes upregulated in T2R and repressed in T1R. In both T1R and T2R, ML2388 was significantly upregulated. This gene encodes a probable membrane protein and epitope prediction using Bepipred-2.0 revealed a distinct B-cell epitope. Overexpression of ML2388 was noted consistently across the reaction samples. From the host immune gene expression profiles, genes for CXCL9, CXCL10, CXCL2, CD40LG, IL17A and CXCL11 were upregulated in T1R when compared to the NR. In T2R, CXCL10, CXCL11, CXCL9, CXCL2 and CD40LG were upregulated when compared to the NR group. Conclusion A gene set signature involving bacterial genes ML2388, ML2664, and host immune genes CXCL10 and IL-17A can be transcriptomic markers for reactional states in leprosy.


Introduction
Mycobacterium leprae (M. leprae), the causative bacillus for leprosy, continues to infect endemic populations in tropical countries, with approximately 200,000 new cases of leprosy emerging each year globally. Mycobacterium leprae infects the skin and the peripheral nerves causing skin lesions with loss of sensation resulting from demyelinating neuropathy as the bacilli infect the Schwann cells of the axonal myelin in the peripheral neurons (Richardus et al., 1996). Nerve damage in leprosy is mediated by M. leprae infection of the Schwann cells as well as exacerbated immune responses in the human host. Leprosy is manifested with a complex host immunological profile that classifies the disease into a cell-mediated immunity (CMI)-dominated tuberculoid pole (TT) and the humoral immune (HI) response-regulated lepromatous pole (LL; WHO Expert Committee on leprosy, World Health Organization, 2012). Both these poles are separated by three borderline intermediary groups that gradient from CMI towards the HI. These include the borderline-tuberculoid (BT), mid-borderline (BB) and marginal lepromatous forms (BL; Ridley and Jopling, 1966).
About 30-40% of leprosy infected individuals in the borderline forms and rarely in the polar states manifest delayed-type hypersensitivity reactions, the type 1 reaction also known as reversal reaction and the type 2 reaction known as Erythema Nodusum Leprosum (ENL; Kahawita and Lockwood, 2008;Montoya and Modlin, 2010). These inflammatory responses can occur before, during and after the treatment with multidrug therapy (MDT) and are managed by immunomodulatory drugs in high doses that often contribute to morbidity. Reactional states are a significant cause of nerve damage and associated disability in leprosy. Early detection of reactional episodes can facilitate prophylactic treatment interventions that minimize the risk of nerve damage (Kwittken, 1968;Sehgal, 1987;Nery et al., 2013).
Predictive genomic, transcriptomic and host immune biomarkers can play a critical role in detecting subclinical nerve damage and determining factors that trigger reactional states in leprosy. In leprosy endemic tropical countries, it is often challenging to characterize an individual's immune background due to varied antigenic exposure (Yuan et al., 2021). Thus, attributing specific immune responses (cytokine and antibody quantities or human immune gene expression signatures) alone to the onset of reactional states or leprosy per se may offer limited applicability in developing effective diagnostics for these M. leprae specific immune exacerbations in leprosy (Leal-Calvo et al., 2021).
A correlative gene expression signature that originates from both M. leprae and human host immune system provides comprehensive predictive and prognostic information for determining the onset of these inflammatory responses (Teles et al., 2013;Manry et al., 2017;Montoya et al., 2019;Leal-Calvo et al., 2021). In this study, we conducted a cross sectional analysis to quantify relative abundance of M. leprae and host immune gene transcripts in localized skin lesions of leprosy cases with type I and type II reactions. A gene expression signature that demonstrates significant association with reactional states has been determined. Follow up studies are warranted to validate these expression signatures in a longitudinal cohort (Tió-Coma et al., 2021).

Sample size
A total of 60 newly diagnosed untreated leprosy cases were recruited at the outpatient department of Schieffelin Institute of Health Research and Leprosy Centre in Karigiri, India. Following institutional ethical clearance, informed and written consent for participation was obtained from each subject prior to recruitment in the study following the ethical guidelines as laid down by the Indian Council of Medical Research. The sample was stratified as 20 with Type 1 Reaction, 20 with Type 2 reaction and 20 without any reaction. Post clinical examination, 5 mm x 5 mm excisional skin biopsies from skin lesions were collected of subjects for all the study groups. Clinical details of the sample were provided in Table 1. The study design and experiments were depicted in Figure 1.

Mycobacterium Leprae whole transcriptome hybridization arrays
Total RNA was extracted following the Trizol protocol (Qiagen RNeasy Lipi Tissue kit-Cat#74804) and bacterial RNA was enriched in the samples. The quality of RNA was estimated using BioAnalyzer 2,100 (Agilent Technologies) followed by labelling, reverse transcription, amplification and hybridization to the arrays. Based on the BioAnalyzer reports for RNA quality estimations,16 NR, 18 T1R and 19 T2R samples were found suitable for whole transcriptome hybridizations.
A 2x400K gene expression array (whole-genome tiling array) was designed with the probes having 60-mer oligonucleotides tiling every 10 bp of the genome sequence of M. leprae (NC_011896.1) by Genotypic Technology Pvt. Ltd. (Bengaluru, India). The array comprised 420,288 features which include probes and Agilent controls. The samples for gene expression were labelled using the Agilent Quick-Amp labelling Kit (p/n5190-0442). 500 ng each of total RNA was reverse transcribed at 40°C using oligo dT primer tagged to a T7 polymerase promoter and converted to double-stranded cDNA. Synthesized double-stranded cDNA were used as templates for cRNA generation. cRNA was generated by in vitro transcription and the dye Cy3 CTP (Agilent) was incorporated during this step. The cDNA synthesis and in vitro transcription steps were carried out at 40°C. Labelled cRNA was cleaned up using Qiagen RNeasy columns (Qiagen,Cat No: 74106) and quality was assessed for yields and specific activity using the Nanodrop ND-1000. The hybridized slides were scanned on a G2600D scanner (Agilent Technologies). The data thus acquired is analyzed using GeneSpring GX Version 12.1 software. Data were normalized and a fold difference in expression was noted from 359,922 probes which include sense and antisense orientations of 179,961 probes. The differentially expressing M. leprae genomic regions between type 1, type 2 reactions and non-reactional cases were noted. All the samples were performed in technical replicates to validate the observations and microarray data corresponding to 359,922 probes for each of the samples (Figure 1). The data as well as the array design was uploaded to the Gene Expression Omnibus (GEO) Repository of the National Center for Biotechnology Information (NCBI) with the accession numbers: GSE85948 and GPL22363.   The schematic representation for the cross-sectional analysis of Mycobacterium leprae transcriptome (A), human immune gene expression (B) and corresponding circulatory levels of cytokines, interleukins and chemokines (C) to derive a gene expression signature for reactional states of leprosy (D  The 75th percentile ranking was used to normalize the probe intensities. The fold difference in expression was noted by subtracting the gene intensities of reactional samples from that of non-reactional samples in each experiment using the geometric mean of the technical replicates. These fold changes were log-transformed to base 2 and volcano plots were generated to identify differentially expressed genes (DEGs). The fold change of ≥ 0.6 was considered as upregulated and ≤ − 0.6 was considered as down-regulated. The Benjamin Hochberg adjusted p values were represented as -log 10 (p value; Rajkumar et al., 2015).

RT 2 PCR profiler arrays
RT 2 PCR profiler arrays (Qiagen Inc., United States) for human inflammatory cytokines and receptors (PAHS-011Z) were used to quantitate expression levels of 96 human immune genes in the lesional skin RNA samples across the study group. Each catalogued RT 2 Profiler PCR Array contains a list of the human inflammatory cytokines and receptors genes as well as five housekeeping (reference) genes on the array. In addition, each array contains a panel of proprietary controls to monitor genomic DNA contamination (GDC) as well as the first strand synthesis (RTC) and real-time PCR efficiency (PPC). The list of genes was provided in Qiagen array Cat. no. PAHS-011Z. Total RNA was isolated from skin biopsy specimens using RNeasy kit (Qiagen Cat No: Cat. No./ID: 74104) according to the manufacturer's instructions. RNA quality was determined using a Nanodrop and was reverse transcribed using a QuantiTect Reverse Transcription Kit (Cat No: Cat. No./ID: 205311). The cDNA was used on the real-time RT 2 Profiler PCR Array (Cat. no. PAHS-011Z) in combination with RT 2 SYBR Green qPCR Mastermix (Qiagen Cat. no. 330529). Fold-change values greater than one indicates a positive-or an up-regulation, and the fold-regulation is equal to the fold-change. Fold-change values less than one indicate a negative or down-regulation, and the foldregulation is the negative inverse of the fold-change. The value of ps were calculated based on a Student's t-test of the replicate 2 (−Delta Ct) values for each gene in the control group and treatment groups. The data was analyzed using Qiagen GeneGlobe application for RT 2 PCR profiler arrays.

Multiplex qPCR
RNA extracted from lesional skin biopsy specimens (NR = 16; T1R = 16; T2R = 9) was used in the qPCR assays. The numbers were different from the size of the study groups as these are the samples in which RNA quality is optimal for cDNA preparation and qPCR. The genes GNLY, CD8A, CXCL10, IFI6, IL10, PRF1, CCL2, FCGR1B, OAS1, IFI44, and CTLA4 have been amplified using the conditions described elsewhere (Supplementary Material-1; Supplementary Table S3; Tió-Coma et al., 2021). The data of the qPCR were analyzed using Thermo Fisher Cloud and GraphPad Prism version 8.0.1. Samples were analyzed in duplicates and ΔCts were calculated using GAPDH as the reference gene. Mann-Whitney U-test was performed to determine if differences in gene expression is statistically significant.

Statistical analysis
For the transcriptome data, normalization of DEGs and the expression threshold of ± 0.6 for log 2 fold_change (log 2 FC) was determined using Agilent Gene Spring GX software (Agilent Inc.; Supplementary Table 1). The principal component analysis was performed using in-built prcomp () function in R and plots were generated using ggplot2 package. The functional GO term enrichment analysis was performed using DAVID database (Sherman et al., 2022) and software.

Results
Analysis of transcriptome-wide changes using principal component analysis Considering 1,600 genes and 45 rRNA transcripts whose intensities were noted from the hybridization arrays, the normalized data with log 2 FC values were subjected to principal component (PC) analysis to estimate the sample variance and reduce the dimensionality in the data. We first determined if the number of PCs are sufficient to explain the fraction of variance using a Pareto chart. A sequential reduction in variance across PCs was noted with PC 53 being zero indicating that the sample numbers are sufficient to explain variance ( Figure 2A). Further we visualized the clusters using ggfortify () package in R and plotted the clusters from the PCs ( Figure 2B).
Cluster 1 largely represents NR, Cluster 2 T1R, Cluster 3 T2R samples ( Figure 2B). Given the diverse Ridley Jopling (RJ) and bacteriological index (BI) classification across the NR samples, we did a differential gene expression analysis of the NR samples that have a BI of zero in comparison to those that have positive BI. We noted that genes ML0247(putative arsenate reductase), ML2269 (putative hydrolase), ML2296 (putative membrane protein) and ML1182(PPEfamily protein) are over expressed in BI-zero NR samples whereas ML1466 (50S ribosomal protein L27) and ML1180 (Putative ESAT-6like protein X) are upregulated in BI positive samples (Supplementary Material-1; Supplementary Figure S1).

Differentially expressed genes of Mycobacterium leprae across the reactional states
From the log 2 FC values, we noted transcripts corresponding to 132 genes of M. leprae for T1R and 117 genes in T2R that are significantly upregulated in comparison to NR. In both the reactional states, 70 genes were upregulated, and 38 genes were downregulated ( Figure 3A). Benjamin Hochberg adjusted value of ps were < 0.05 for all these associations. We identified only one gene (ML2664) that was upregulated in T1R and downregulated in T2R and no DEGs in converse (Supplementary Table 2). The top 10 upregulated and the lower 10 downregulated genes were labelled in the volcano plots in Figures 3B,C. The NR sample was further split based on the Ridley Joplin classification into two groups-the TT/BT group and the BL/ LL group. The TT/BT group was compared with T1R and BL/LL with T2R. We noted that the top 10 DEGs are retained the same as for the

Unsupervised clustering analysis for expression patterns
Hierarchical clustering analysis with the z scores of the fold changes gene-wise across the study groups revealed clusters with various enriched GO terms. Among the upregulated genes in T1R, we noted over representation of genes that encode integral membrane proteins, followed by cytosolic and ribosome bound protein coding genes ( Figure 4A). From the GO biological processes genes corresponding to proteins that mediate cell wall biosynthesis, lipid biosynthetic pathways, drug transport, fatty acid and amino acid metabolism, and the biotin and folic acid biosynthetic pathways are overrepresented ( Figure 4B). In the T2R, among the genes that had GO annotations for cellular component, those that encode integral Frontiers in Microbiology 07 frontiersin.org membrane components were higher in number ( Figure 4C) and from the GO biological processes ( Figure 4D), genes involved in translation, DNA recombination, fatty acid metabolism, protein transport and amino acid metabolism are present.

Functional term analysis from GO annotations
We further used the enriched GO terms from the DAVID database to ascertain the probability of gene co-occurrence in each GO term across the sample. Among the upregulated genes in T1R in comparison (A,B) represent heatmap of the genes that are upregulated in T1R in comparison to NR and (C,D) includes genes that are upregulated in T2R in comparison to NR. All the genes listed are significantly overexpressed (p < 0.05) and the color key represents the Log2FC values for each gene ranging from −2 to 2. The blank lines in GO terms are genes without GO annotations. In both the figures, we retained both T1R and T2R gene expressions to visualize expression changes across the sample types.

Expression profiles of known Mycobacterium leprae antigens from the IEDB database
From the 61 antigens recorded in Immune Epitope Database (IEDB) for M. leprae (identified by the search term Bacillus leprae), we selected 58 protein antigens which has Uniprot IDs (Table 3) and studied their expression across the sample. The significantly upregulated antigen coding genes in T1R and T2R are depicted in Figures 6A,B, respectively.

Consistent overexpression of ML2388 across the sample in reactional states of leprosy
Among all the significant DEGs noted across the sample, we have seen consistent over expression of ML2388 across the reactional states T1R and T2R in comparison to NR. This gene encodes a possible membrane protein and has GO term for cellular localization as the integral component of the membrane. We predicted the presence of linear B cell epitopes in this protein using Bepipred 2.0. The results were presented in Figure 7. Three linear epitopes were predicted with high confidence in the exposed and helical or coiled regions of the protein (Reece et al., 2006).

Differential expression of human immune genes
From the RT2 PCR Profiler arrays, the expression data was analyzed using Qiagen GeneGlobe pipeline using the ∆∆Ct method. The threshold cycle values across the study groups were normalized using the house keeping gene [GAPDH (Glyceraldehyde-3-phosphate dehydrogenase)] expression levels (Kwittken, 1968). The quality checks were performed by measuring the PCR array reproducibility, ∆Ct of the reverse transcription control and genomic DNA contamination. For the test sets (T1R and T2R), the PPC was 15 which indicates that the array reproducibility check has been passed, the transcription control has a ∆Ct of 5.28 and ∆Ct of genomic DNA contamination is 33 which indicates that the contamination is minimal. The cut off Ct value is set to 33. Fold-Change (2 (−∆∆Ct) ) is the normalized gene expression (2 (−∆Ct) ) in the Test Sample divided the normalized gene expression (2 (− ∆Ct) ) in the Control Sample. Fold-Regulation represents fold-change results in a biologically meaningful way.

Discussion
Reactional states in leprosy pose a significant challenge to the global efforts to contain the disease. These immune exacerbations are Frontiers in Microbiology 10 frontiersin.org a major cause of nerve function impairment and consequent disabilities in leprosy. Various host-related factors have been reported as risk factors for T1R and these include age, severity of the disease, and having a positive BI in the slit skin smears. Lepromatous leprosy with the BI of more that 4 + is identified as a risk factor for T2R in leprosy (Pandhi and Chhabra, 2013). Few studies have also implicated antigenic triggers in Type 1 reaction, leading to expansion of both cross-reactive and specific T-cells. Despite limited genotypic variability between various strains of M. leprae, the patterns of infection, pathogenicity and virulence largely differ across various clinical phenotypes. This suggests that success of infection and leprosy progression is largely dependent on host′s immune response and genetic complement (Texereau et al., 2005). Studies on M. leprae and host immune gene expression signatures can help identify potential biomarkers to predict reactional states in leprosy. In this study, we conducted a case-control association analysis of DEGs in M. leprae transcriptome in the lesional skin tissues of leprosy cases in T1R and T2R and compared them with those without any hypersensitivity reactions to decipher a gene expression signature for these reactional states. From the GEO datasets for M. leprae, while there are several studies on understanding host transcriptomic response (Teles et al., 2013;Manry et al., 2017;Zhang et al., 2019;Leal-Calvo et al., 2021) to M. leprae infection, we identified only one other study (Montoya et al., 2019) that profiled differential gene expression within M. leprae Over-represented GO terms in the T1R Group (p < 0.05) in comparison to NR (A) and in the T2R Group (p < 0.05) in comparison to NR group (B). The number of genes on x-axis indicates genes that have the same GO term in each of the classes (biological process, cellular localization, and molecular function).
Frontiers in Microbiology 11 frontiersin.org transcriptome to ascertain their implications on reactional states. Although, our experiments were conducted with total RNA, we limited our analysis to expression profiles in the protein coding genes. While identifying DEGs across study groups remained the main objective, we attempted to delineate the possible functional implications of DEGs on the reactional outcomes. The gene chip array from Agilent that we employed in the experiments is available at the GEO repository (GPL22363). From the upregulated genes in T1R (n = 132), genes ML2064, ML1271, ML1960, ML1220, ML2498, ML1996, ML2388, ML0429, ML2030 and ML0224 are the top 10 genes with highest Log2FC values. Most of these genes code for integral or intrinsic membrane proteins while others encode proteins of the 30S ribosomal subunit and enzymes involved in glycolysis and possible enoyl-CoA hydratases. Additionally, we noted 27 genes that encode integral and intrinsic components of the cell wall and plasma membrane (Supplementary Table S2). We noted 12 virulence genes among the upregulated which include ML0114-an ABC O-antigen transporter, ML1014-RNA polymerase sigma factor sigB, ML1076-RNA (Williams et al., 2004(Williams et al., , 2007 polymerase sigma factor SigE, ML1128-Diaminopimelate decarboxylase, ML1220-Biotin synthase (Lastória and de Abreu, 2014), ML1547-4′-phosphopantetheinyl transferase, ML1656-3-oxoacyl-[acyl-carrier-protein] synthase, ML1675-Uracil-DNA glycosylase, ML2124-Sensor-type histidine kinase, ML2307-Transcriptional regulator and ML2350-ATP-dependent efflux pump essential for phthiocerol dimycocerosates translocation (Sharma et al., 2009;Orlova et al., 2013). Overexpression of virulence genes were noted in several other studies in reactional states (Kahawita and Lockwood, 2008;Nery et al., 2013). We also noted over expression of heat shock proteins (hsp18) in T1R (Lini et al., 2009).
In T2R, genes ML2498, ML1526, ML0394, ML1960, ML2388, ML0429, ML0281, ML1847, ML1618 and ML1271 were the Top-10    Predicted linear B-Cell epitopes in the possible membrane protein-ML2388.  (Hasan et al., 2004), ML2307-Transcriptional regulator and ML2439-Sensory transduction protein RegX3 (Duthie et al., 2007;Sapkota et al., 2010;Tang et al., 2021). Among these we noted consistent over expression of ML2388. This gene encodes a putative membrane protein and possess three regions in the sequence that are predicted to be potential B-cell epitopes. To further understand the potential implications of previously characterized antigen coding genes in M. leprae with transcriptomic responses to reactions in leprosy, we measured their expression across the study groups and noted over-expression of integrase (ML0008), probable secreted protein (ML0757), phospotidylgycerol (ML1274) and 18Kda Heat shock Protein (ML1795) with T1R. In T2R, we noted overexpression of Methyltransf_25 domain-containing protein (ML0394).

Frontiers in
From the host genes, we noted that between the NR and the T1R groups, a significant difference in expression was observed for genes GNLY (Granulysin) responsible for transfer of granzymes (a family of serine proteases traditionally known for their role in promoting cytotoxicity of foreign, infected or neoplastic cells), CD8A (Cytotoxic T-lymphocytes) important for immune defense against intracellular pathogens, including viruses and bacteria, CXCL10 (Interferon gamma induced protein-10). Granulysin has been noted to express high in tuberculoid lesions (Geluk et al., 2014). Upregulation of CD8 antigen in reversal reactions has also been noted in leprosy/HIV co-infections (de Oliveira et al., 2013). The panel depicts immune gene expressions across the study groups as determined by qPCR. Statistically significant differences in mean of the ∆Ct for technical replicates of each of the sample was denoted by * for value of p of 0.05 and ** for value of p of 0.005.

Frontiers in Microbiology 14 frontiersin.org
Several studies revealed the upregulation of CXCL10 (IP10) in reactional states of leprosy (Chaitanya et al., 2013;van Hooij et al., 2016;Ferreira et al., 2021). Expression of IP-10 is seen in many Th1-type inflammatory diseases, where it is thought to play an important role in recruiting activated T cells into sites of tissue inflammation. Interleukin 10 (IL-10), a cytokine with potent anti-inflammatory properties, plays a central role in limiting host immune response to M. leprae, thereby preventing autoimmune inflammation and maintains tissue homeostasis. PRF1 (Perforin 1) is a glycoprotein responsible for pore formation in cell membranes of target cells and CCL2 chemokine which controls immunity by promoting regulatory T cell communication with dendritic cells in lymph nodes, are all upregulated in T1R. Further, FCGR1B (Fc Gamma Receptor Ib) which functions by binding to the Fc regions of immunglobulin, OAS1 (2′-5′-Oligoadenylate Synthetase 1) that helps in degradation of viral infections, IFI44 (Interferon Induced Protein 44) and CTLA4 (cytotoxic T-lymphocyte-associated protein 4) which acts as a halting mechanism and reduces the function of T cells are over expressed in T1R. Similar results were noted with RT2 profiler arrays in this study. CXC chemokines, CXCL9, CXCL10, CXCL2, CXCL11, CD40 ligand (CD40LG), and interleukin IL17A were upregulated in T1R. In T2R, CXC chemokines CXCL10, CXCL11, CXCL9, CXCL2 and CD40 ligand (CD40LG) were upregulated (Geluk et al., 2014;Teles et al., 2019;Tió-Coma et al., 2021).
In conclusion, we recommend bacterial genes ML2388, ML2664, and host immune genes CXCL10 and IL-17A as transcriptomic signatures for reactional states in leprosy. Our study profiled various upregulating and downregulating gene signatures from Mycobacterium leprae and human immune genes that demonstrated plausible association with reactional states in leprosy. Further studies are however required to validate the above identified gene expression signatures as predictive markers for leprosy reactions in a longitudinal cohort.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary material.

Ethics statement
The studies involving human participants were reviewed and approved by Karigiri Research Committee. The patients/participants provided their written informed consent to participate in this study.

Author contributions
MD, SV, and AG conceived the study. MT-C, DD, AV, and IH contributed to the data curation. MD, SV, and AV contributed to the formal analysis. SV and AG acquired funding. MD, SV, AV, DD, and MT-C investigated the study. SV, MD, and AG supervised the study. MD and SV wrote the original draft. AG, AV, and MT-C reviewed, and edited the manuscript. All authors reviewed, discussed, and agreed with the final manuscript.

Funding
This work received the financial support from the Leprosy Research Initiative (LRI) and the Turing Foundation under LRI Grant number 704.16.57.