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

Materials and methods

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.

TABLE 1

Table 1. Clinical and demographic characteristics of the study sample.

FIGURE 1

Figure 1. 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). (Created with BioRender.com).

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 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₁₀ (p value; Rajkumar et al., 2015).

RT² PCR profiler arrays

RT² 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²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² Profiler PCR Array (Cat. no. PAHS-011Z) in combination with RT²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 fold-regulation 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² 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₂fold_change (log₂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₂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).

FIGURE 2

Figure 2. (A) A pareto plot with 80% of the fraction variance explained by the first 5 PCs using PCs for each of the sample and (B). Clusters mapped between PC1 and PC2 in a 2-dimentional PCA plot. The clusters contained 16 non reactional controls (NR), 18 T1R and 19 T2R samples.

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(PPE-family protein) are over expressed in BI-zero NR samples whereas ML1466 (50S ribosomal protein L27) and ML1180 (Putative ESAT-6-like 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₂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 whole data set even after splitting the control sample based on RJ classification (Supplementary Material-1; Supplementary Figure S2).

FIGURE 3

Figure 3. (A) Venn diagram representing the number of significantly upregulated and downregulated genes of Mycobacterium leprae across the study groups. (B,C) Volcano plots with log₂FC (derived from normalized gene expression data) showing the top 10 differentially expressed genes in T1R and T2R compared to NR cases. The overexpressed genes are those with the volcano plot p-values < 0.05 and Log₂FC values of ≥ 0.6 (mean of each of the reaction groups) and underexpressed genes are those with the volcano plot p-values < 0.05 and Log₂FC values of ≤ − 0.6 (for the mean of each of the reaction groups). Scales in B and C represent the Log₂FC values.

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

FIGURE 4

Figure 4. (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.

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 to NR, for GO cellular component terms for integral component of membrane [GO:0016021], intrinsic component of membrane [GO:0031224], the membrane part [GO:0044425] and membrane [GO:0016020] were significantly overrepresented (p-value 0.0004). In the GO biological processes, terms for transcription and regulation of transcription, RNA biogenesis, regulation of RNA biosynthetic processes and lipid metabolism were overrepresented among others noted in Table 2 and Figure 5A. For the GO Molecular Function, genes involved in DNA binding and hydrolase activity are noted. For genes that are overexpressed in T2R in comparison to NR, the GO cellular component terms for integral component of membrane [GO:0016021], intrinsic component of membrane [GO:0031224], the membrane part [GO:0044425] and membrane [GO:0016020] were significantly overrepresented (p-value 0.008). For the Biological Process terms, lipid metabolic processes, cellular component biogenesis, RNA processing, RNA metabolic process, ribonucleoprotein complex biogenesis and genes involved in signal transduction were significantly overexpressed along with others shown in Table 2 and Figure 5B.

TABLE 2

Table 2. The enriched GO terms in T1R and T2R groups across the sample.

FIGURE 5

Figure 5. 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).

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.

TABLE 3

Table 3. The log₂FC for DEGs among the antigen coding genes in immune epitope database (IEDB) database for Mycobacterium leprae.

FIGURE 6

Figure 6. (A) Differentially expressed genes (DEGs) among antigen coding genes in T1R. (B) DEGs among antigen coding genes in T2R.

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

FIGURE 7

Figure 7. Predicted linear B-Cell epitopes in the possible membrane protein—ML2388.

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.

We noted a significant upregulation of CXC chemokines, CXCL9, CXCL10, CXCL2, CXCL11, CD40 ligand (CD40LG), and interleukin IL17A in T1R. In T2R, CXC chemokines, CXCL10, CXCL11, CXCL9, CXCL2 and CD40 ligand (CD40LG) were upregulated (Scollard et al., 2011).

Multiplex qPCR assays

We also tested the expression levels of a selected set of human immune genes that were observed elsewhere (Teles et al., 2013; Geluk et al., 2014; Tió-Coma et al., 2019) to have associations with reactional states in leprosy using the multiplex qPCR system in our study groups [T1R (n = 16), T2R (n = 9) and NR (n = 16)]. Between the NR and the T1R groups, a significant difference in expression was observed for genes GNLY (Granulysin), CD8A, CXCL10, IL10, PRF1 (Perforin 1), CCL2, FCGR1B (Fc Gamma Receptor Ib), OAS1 (2′-5’-Oligoadenylate Synthetase 1), IFI44 (Interferon Induced Protein 44) and CTLA4 (cytotoxic T-lymphocyte-associated protein 4). Gene expression was higher (lower ∆Ct) for all the genes tested in patients with T1R than NR patients. No statistically significant differences in expression were noted between NR and T2R or T1R and T2R (Figure 8).

FIGURE 8

Figure 8. 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.

Discussion

Reactional states in leprosy pose a significant challenge to the global efforts to contain the disease. These immune exacerbations are 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 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 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 significantly upregulated genes. These encode conserved membrane proteins, proteins of the 30S and 50S ribosomal subunits, possible enoyl-CoA hydratases and other enzymes. In addition, we noted that 10 virulence genes were significantly over expressed. They are ML0243—Putative acyl-CoA synthetase, ML1014—RNA polymerase sigma factor sigB, ML1076—RNA polymerase sigma factor SigE, ML1128—Diaminopimelate decarboxylase, ML1633—Possible secreted hydrolase, ML1727—O-phosphoserine phosphohydrolase(Acebrón-García-de-Eulate et al., 2021), ML1925—Superoxide dismutase [Cu-Zn], ML1954—Pantothenate kinase(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).

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

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.

Acknowledgments

The authors would like to thank all the administration and laboratory staff at Schieffelin Institute of Health-Research and Leprosy Center who were involved in sample collection and performing the experiments and Leiden University Medical Centre, Netherlands for also performing the experiments. The authors would like to thank Genotypic Technologies Pvt. Ltd. India for the Microarray experiments. Special thanks to the administration of SIH-R&LC, Karigiri for the infrastructural support throughout the study and for the financial support from the Leprosy Research Initiative (LRI) and the Turing Foundation under LRI Grant number 704.16.57.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2023.1113318/full#supplementary-material

References

Acebrón-García-de-Eulate, M., Blundell, T. L., and Vedithi, S. C. (2021). Strategies for drug target identification in Mycobacterium leprae. Drug Discov. Today 26, 1569–1573. doi: 10.1016/j.drudis.2021.03.026

CrossRef Full Text | Google Scholar

Chaitanya, V. S., Lavania, M., Nigam, A., Turankar, R. P., Singh, I., Horo, I., et al. (2013). Cortisol and proinflammatory cytokine profiles in type 1 (reversal) reactions of leprosy. Immunol. Lett. 156, 159–167. doi: 10.1016/j.imlet.2013.10.008

CrossRef Full Text | Google Scholar

de Oliveira, A. L., Amadeu, T. P., de França Gomes, A. C., Menezes, V. M., da Costa Nery, J. A., Pinheiro, R. O., et al. (2013). Role of CD8+ T cells in triggering reversal reaction in HIV/leprosy patients. Immunology 140, 47–60. doi: 10.1111/imm.12108

CrossRef Full Text | Google Scholar

Duthie, M. S., Goto, W., Ireton, G. C., Reece, S. T., Cardoso, L. P. V., Martelli, C. M. T., et al. (2007). Use of protein antigens for early serological diagnosis of leprosy. Clin. Vaccine Immunol. 14, 1400–1408. doi: 10.1128/CVI.00299-07

CrossRef Full Text | Google Scholar

Ferreira, H., Mendes, M. A., de Mattos Barbosa, M. G., de Oliveira, E. B., Sales, A. M., Moraes, M. O., et al. (2021). Potential role of CXCL10 in monitoring response to treatment in leprosy patients. Front. Immunol. 12:662307. doi: 10.3389/fimmu.2021.662307

CrossRef Full Text | Google Scholar

Geluk, A., van Meijgaarden, K. E., Wilson, L., Bobosha, K., van der Ploeg-van Schip, J. J., van den Eeden, S. J. F., et al. (2014). Longitudinal immune responses and gene expression profiles in type 1 leprosy reactions. J. Clin. Immunol. 34, 245–255. doi: 10.1007/s10875-013-9979-x

CrossRef Full Text | Google Scholar

Hasan, Z., Mahmood, A., Zafar, S., Khan, A. A., and Hussain, R. (2004). Leprosy patients with lepromatous disease have an up-regulated IL-8 response that is unlinked to TNF-alpha responses. Int. J. Lepr. Other Mycobact. Dis. 72, 35–44. doi: 10.1489/1544-581X(2004)072<0035:LPWLDH>2.0.CO;2

CrossRef Full Text | Google Scholar

Kahawita, I. P., and Lockwood, D. N. J. (2008). Towards understanding the pathology of erythema nodosum leprosum. Trans. R. Soc. Trop. Med. Hyg. 102, 329–337. doi: 10.1016/j.trstmh.2008.01.004

CrossRef Full Text | Google Scholar

Kwittken, J. (1968). Erythema nodosum leprosum. Int. J. Lepr. Other Mycobact. Dis. 36, 94–95

Google Scholar

Lastória, J. C., and de Abreu, M. A. M. M. (2014). Leprosy: a review of laboratory and therapeutic aspects—part 2. An. Bras. Dermatol. 89, 389–401. doi: 10.1590/abd1806-4841.20142460

CrossRef Full Text | Google Scholar

Leal-Calvo, T., Avanzi, C., Mendes, M. A., Benjak, A., Busso, P., Pinheiro, R. O., et al. (2021). A new paradigm for leprosy diagnosis based on host gene expression. PLoS Pathog. 17:e1009972. doi: 10.1371/journal.ppat.1009972

CrossRef Full Text | Google Scholar

Leal-Calvo, T., Martins, B. L., Bertoluci, D. F., Rosa, P. S., de Camargo, R. M., Germano, G. V., et al. (2021). Moraes MO. Large-scale gene expression signatures reveal a Microbicidal pattern of activation in mycobacterium leprae-infected monocyte-derived macrophages with low multiplicity of infection. Front. Immunol. 12:647832. doi: 10.3389/fimmu.2021.647832

CrossRef Full Text | Google Scholar

Lini, N., Shankernarayan, N. P., and Dharmalingam, K. (2009). Quantitative real-time PCR analysis of mycobacterium leprae DNA and mRNA in human biopsy material from leprosy and reactional cases. J. Med. Microbiol. 58, 753–759. doi: 10.1099/jmm.0.007252-0

CrossRef Full Text | Google Scholar

Manry, J., Nédélec, Y., Fava, V. M., Cobat, A., Orlova, M., Thuc, N. V., et al. (2017). Deciphering the genetic control of gene expression following mycobacterium leprae antigen stimulation. PLoS Genet. 13:e1006952. doi: 10.1371/journal.pgen.1006952

CrossRef Full Text | Google Scholar

Montoya, D. J., Andrade, P., Silva, B. J. A., Teles, R. M. B., Ma, F., Bryson, B., et al. (2019). Dual RNA-Seq of human leprosy lesions identifies bacterial determinants linked to host immune response. Cell Rep. 26, 3574–3585.e3. doi: 10.1016/j.celrep.2019.02.109

CrossRef Full Text | Google Scholar

Montoya, D., and Modlin, R. L. (2010). Learning from leprosy: insight into the human innate immune response. Adv. Immunol. 105, 1–24. doi: 10.1016/S0065-2776(10)05001-7

CrossRef Full Text | Google Scholar

Nery, J. A. C., Bernardes Filho, F., Quintanilha, J., Machado, A. M., Oliveira, S. S. C., and Sales, A. M. (2013). Understanding the type 1 reactional state for early diagnosis and treatment: a way to avoid disability in leprosy. An. Bras. Dermatol. 88, 787–792. doi: 10.1590/abd1806-4841.20132004

CrossRef Full Text | Google Scholar

Orlova, M., Cobat, A., Huong, N. T., Ba, N. N., van Thuc, N., Spencer, J., et al. (2013). Gene set signature of reversal reaction type I in leprosy patients. PLoS Genet. 9:e1003624. doi: 10.1371/journal.pgen.1003624

CrossRef Full Text | Google Scholar

Pandhi, D., and Chhabra, N. (2013). New insights in the pathogenesis of type 1 and type 2 lepra reaction. Indian J. Dermatol. Venereol. Leprol. 79, 739–749. doi: 10.4103/0378-6323.120719

CrossRef Full Text | Google Scholar

Rajkumar, A. P., Qvist, P., Lazarus, R., Lescai, F., Ju, J., Nyegaard, M., et al. (2015). Experimental validation of methods for differential gene expression analysis and sample pooling in RNA-seq. BMC Genomics 16:548. doi: 10.1186/s12864-015-1767-y

CrossRef Full Text | Google Scholar

Reece, S. T., Ireton, G., Mohamath, R., Guderian, J., Goto, W., Gelber, R., et al. (2006). ML0405 and ML2331 are antigens of mycobacterium leprae with potential for diagnosis of leprosy. Clin. Vaccine Immunol. 13, 333–340. doi: 10.1128/CVI.13.3.333-340.2006

CrossRef Full Text | Google Scholar

Richardus, J. H., Finlay, K. M., Croft, R. P., and Smith, W. C. (1996). Nerve function impairment in leprosy at diagnosis and at completion of MDT: a retrospective cohort study of 786 patients in Bangladesh. Lepr. Rev. 67, 297–305. doi: 10.5935/0305-7518.19960030

CrossRef Full Text | Google Scholar

Ridley, D. S., and Jopling, W. H. (1966). Classification of leprosy according to immunity. A five-group system. Int. J. Lepr. Other Mycobact. Dis. 34, 255–273

Google Scholar

Sapkota, B. R., Macdonald, M., Berrington, W. R., Misch, E. A., Ranjit, C., Siddiqui, M. R., et al. (2010). Association of TNF, MBL, and VDR polymorphisms with leprosy phenotypes. Hum. Immunol. 71, 992–998. doi: 10.1016/j.humimm.2010.07.001

CrossRef Full Text | Google Scholar

Scollard, D. M., Chaduvula, M. V., Martinez, A., Fowlkes, N., Nath, I., Stryjewska, B. M., et al. (2011). Increased CXC ligand 10 levels and gene expression in type 1 leprosy reactions▿. Clin. Vaccine Immunol. 18, 947–953. doi: 10.1128/CVI.00042-11

CrossRef Full Text | Google Scholar

Sehgal, V. N. (1987). Reactions in leprosy. Clinical aspects. Int. J. Dermatol. 26, 278–285. doi: 10.1111/j.1365-4362.1987.tb00188.x

CrossRef Full Text | Google Scholar

Sharma, R., Lavania, M., Chauhan, D. S., Katoch, K., Amresh, P., Rakhi, R., et al. (2009). Potential of a metabolic gene (accA3) of M. leprae as a marker for leprosy reactions. Indian J. Lepr. 81, 141–148

Google Scholar

Sherman, B. T., Hao, M., Qiu, J., Jiao, X., Baseler, M. W., Lane, H. C., et al. (2022). DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res. 50, W216–W221. doi: 10.1093/nar/gkac194

CrossRef Full Text | Google Scholar

Tang, A. S. O., Wong, Q. Y., Yeo, S. T., Ting, I. P. L., Lee, J. T. H., Fam, T. L., et al. (2021). Challenges in managing a lepromatous leprosy patient complicated with Melioidosis infection, Dapsone-induced Methemoglobinemia, hemolytic anemia, and Lepra reaction. Am J Case Rep 22:e931655. doi: 10.12659/AJCR.931655

CrossRef Full Text | Google Scholar

Teles, R. M. B., Graeber, T. G., Krutzik, S. R., Montoya, D., Schenk, M., Lee, D. J., et al. (2013). Type I interferon suppresses type II interferon-triggered human anti-mycobacterial responses. Science 339, 1448–1453. doi: 10.1126/science.1233665

CrossRef Full Text | Google Scholar

Teles, R. M. B., Lu, J., Tió-Coma, M., Goulart, I. M. B., Banu, S., Hagge, D., et al. (2019). Identification of a systemic interferon-γ inducible antimicrobial gene signature in leprosy patients undergoing reversal reaction. PLoS Negl. Trop. Dis. 13:e0007764. doi: 10.1371/journal.pntd.0007764

CrossRef Full Text | Google Scholar

Texereau, J., Chiche, J.-D., Taylor, W., Choukroun, G., Comba, B., and Mira, J.-P. (2005). The importance of toll-like receptor 2 polymorphisms in severe infections. Clin. Infect. Dis. 41, S408–S415. doi: 10.1086/431990

CrossRef Full Text | Google Scholar

Tió-Coma, M., Kiełbasa, S. M., van den Eeden, S. J. F., Mei, H., Roy, J. C., Wallinga, J., et al. (2021). Blood RNA signature RISK4LEP predicts leprosy years before clinical onset. EBioMedicine 68:103379. doi: 10.1016/j.ebiom.2021.103379

CrossRef Full Text | Google Scholar

Tió-Coma, M., van Hooij, A., Bobosha, K., van der Ploeg-van Schip, J. J., Banu, S., Khadge, S., et al. (2019). Whole blood RNA signatures in leprosy patients identify reversal reactions before clinical onset: a prospective, multicenter study. Sci. Rep. 9:17931. doi: 10.1038/s41598-019-54213-y

CrossRef Full Text | Google Scholar

van Hooij, A., Tjon Kon Fat, E. M., Richardus, R., van den Eeden, S. J. F., Wilson, L., de Dood, C. J., et al. (2016). Quantitative lateral flow strip assays as user-friendly tools to detect biomarker profiles for leprosy. Sci. Rep. 6:34260. doi: 10.1038/srep34260

CrossRef Full Text | Google Scholar

WHO Expert Committee on leprosy, World Health Organization . WHO Expert Committee on Leprosy: Eighth Report. Geneva: World Health Organization (2012).

Google Scholar

Williams, D. L., Pittman, T. L., Deshotel, M., Oby-Robinson, S., Smith, I., and Husson, R. (2007). Molecular basis of the defective heat stress response in mycobacterium leprae. J. Bacteriol. 189, 8818–8827. doi: 10.1128/JB.00601-07

CrossRef Full Text | Google Scholar

Williams, D. L., Torrero, M., Wheeler, P. R., Truman, R. W., Yoder, M., Morrison, N., et al. (2004). Biological implications of mycobacterium leprae gene expression during infection. J. Mol. Microbiol. Biotechnol. 8, 58–72. doi: 10.1159/000082081

CrossRef Full Text | Google Scholar

Yuan, Y.-H., Liu, J., You, Y.-G., Chen, X.-H., Yuan, L.-C., Wen, Y., et al. (2021). Transcriptomic analysis of mycobacterium leprae-stimulated response in peripheral blood mononuclear cells reveal potential biomarkers for early diagnosis of leprosy. Front. Cell. Infect. Microbiol. 11:714396. doi: 10.3389/fcimb.2021.714396

CrossRef Full Text | Google Scholar

Zhang, X., Cheng, Y., Zhang, Q., Wang, X., Lin, Y., Yang, C., et al. (2019). Meta-analysis identifies major Histocompatiblity complex loci in or near HLA-DRB1, HLA-DQA1, HLA-C as associated with leprosy in Chinese Han population. J. Invest. Dermatol. 139, 957–960. doi: 10.1016/j.jid.2018.09.029

CrossRef Full Text | Google Scholar