Large-Scale Gene Expression Signatures Reveal a Microbicidal Pattern of Activation in Mycobacterium leprae-Infected Monocyte-Derived Macrophages With Low Multiplicity of Infection

Leprosy is a disease with a clinical spectrum of presentations that is also manifested in diverse histological features. At one pole, lepromatous lesions (L-pole) have phagocytic foamy macrophages heavily parasitized with freely multiplying intracellular Mycobacterium leprae. At the other pole, the presence of epithelioid giant cells and granulomatous formation in tuberculoid lesions (T-pole) lead to the control of M. leprae replication and the containment of its spread. The mechanism that triggers this polarization is unknown, but macrophages are central in this process. Over the past few years, leprosy has been studied using large scale techniques to shed light on the basic pathways that, upon infection, rewire the host cellular metabolism and gene expression. M. leprae is particularly peculiar as it invades Schwann cells in the nerves, reprogramming their gene expression leading to a stem-like cell phenotype. This modulatory behavior exerted by M. leprae is also observed in skin macrophages. Here, we used live M. leprae to infect (10:1 multiplicity of infection) monocyte-derived macrophages (MDMs) for 48 h and analyzed the whole gene expression profile using microarrays. In this model, we observe an intense upregulation of genes consistent with a cellular immune response, with enriched pathways including peptide and protein secretion, leukocyte activation, inflammation, and cellular divalent inorganic cation homeostasis. Among the most differentially expressed genes (DEGs) are CCL5/RANTES and CYP27B1, and several members of the metallothionein and metalloproteinase families. This is consistent with a proinflammatory state that would resemble macrophage rewiring toward granulomatous formation observed at the T-pole. Furthermore, a comparison with a dataset retrieved from the Gene Expression Omnibus of M. leprae-infected Schwann cells (MOI 100:1) showed that the patterns among the DEGs are highly distinct, as the Schwann cells under these conditions had a scavenging and phagocytic gene profile similar to M2-like macrophages, with enriched pathways rearrangements in the cytoskeleton, lipid and cholesterol metabolism and upregulated genes including MVK, MSMO1, and LACC1/FAMIN. In summary, macrophages may have a central role in defining the paradigmatic cellular (T-pole) vs. humoral (L-pole) responses and it is likely that the multiplicity of infection and genetic polymorphisms in key genes are gearing this polarization.

Introduction

Schwann cells in the peripheral nerves and macrophages in the skin are the major host cells for Mycobacterium leprae (ML) infection (1). These cells operate with high plasticity, induced by different environmental factors, and M. leprae has a unique ability to subvert and reprogram these host cells in order to establish a more favorable niche in which to replicate and spread. Huge transcriptomic variations may induce phenotypic modifications, as evidenced by transformations in Schwann to mesenchymal-like cells upon infection (2). Furthermore, increased glucose uptake, mitochondrial shutdown, and lipid biosynthesis resembling the Warburg effect are all phenomena induced upon M. leprae infection in these cells, although some of these are restricted to specific clinical forms (3).

Clinical presentation of leprosy is a spectrum encompassing a myriad of manifestations (1), where the tuberculoid pole (T-pole) is restrictive to bacillus growth leading to localized disease, while the lepromatous forms (L-pole) present a permissive and disseminated clinical form with high bacterial loads (1). M. leprae has suffered a reductive evolution resulting in low genetic variability, which suggests that the diversity of the disease phenotypes is attributable to the host responses (4). This landscape makes leprosy a unique model to understand the mechanisms involved in the immunopathogenesis of infectious diseases.

In the skin, macrophages are pivotal in the host-pathogen interaction, having important roles from proinflammatory and microbicidal activity to tissue remodeling and wound healing, which are features of the so-called M1 and M2 macrophages, respectively. Most of the skin macrophages are derived from monocytes that migrate and differentiate under inflammatory stimuli, referred to as monocyte-derived macrophages (MDMs) (5). Macrophages present huge functional plasticity according to the milieu in order to maintain skin homeostasis (6). Although less efficient in T cell activation than dendritic cells, they are vastly superior in their phagocytic ability (7). Additionally, monocytes engulf the bacillus through phagocytosis and produce cytokines helping to dictate the host-specific immune response at the lesion. The initial immune events involved in leprosy disease progression are probably triggered by the macrophage-M. leprae interaction. This hypothesis can be reinforced by the fact that key innate immune genes, pattern recognition receptors, and autophagic genes have been associated with disease outcome in the mouse model for bacilli replication (infected footpads of Balb/C lineage mice) that carries the NRAMP1 polymorphism. Furthermore, human genome-wide association studies and other genetic analysis has identified PRKN, LRRK2, NOD2, TLR1, and MRC1 as genes associated with disease outcome and are expressed in M. leprae-infected macrophages (8). In this regard, polarized macrophages are found according to the clinical form of leprosy. In tuberculoid lesions, there is a predominance of classically-activated macrophages, which are able to partially contain M. leprae replication by activation of cellular responses, vitamin D-dependent pathways, and granuloma formation consistent with the M1 profile (9). On the other hand, lepromatous lesions present permissive, scavenging, phagocytic, and foamy macrophages with anti-inflammatory profile harboring a large number of bacilli, associated with a poor microbicide activity, which are the phenotypes of M2 macrophages (9, 10).

Pathways such as apoptosis and autophagy, combined with lower levels of proinflammatory cytokines and antigen presentation molecules are inhibited by live, but not dead, M. leprae in macrophages and Schwann cells (11–13). This type of response is triggered by type I interferons (IFNs) through a STING/TBK/IRF3 pathway that inhibits IFN-gamma and other microbicidal mechanisms within infected macrophages in vitro (12).

In a recent re-analysis of leprosy public microarray datasets, mainly comparing skin lesions from different clinical forms, we confirmed previous genes and pathways corroborating the predominance of cellular immunity, leukocyte differentiation, vitamin D receptor (VDR)-mediated microbicidal responses, and granuloma formation at the T-pole, while the L-pole exhibited scavenger receptors and lipid metabolism genes (14). Nevertheless, we also pointed out new differentially expressed genes (DEGs) in leprosy related to skin development and keratinocyte differentiation (14).

Macrophages play a central role in orchestrating gene modulation within the nerve or skin microenvironment. Nevertheless, genome-wide expression patterns from M. leprae-infected macrophages have never been obtained. The gene expression signatures of these macrophages can provide target discoveries for therapies and prevention strategies to interfere with leprosy outcomes. Herein, we show the DEGs of macrophages related to immune and other responses after M. leprae infection. Our data indicate that M. leprae-infected MDMs submitted to a low multiplicity of infection (10:1) for 48 h overly express a cellular immunity profile. This suggests an M. leprae-induced program that might create a granuloma formation type of response in the tissues. Moreover, a comparison between the dataset generated in this study with that of a previously published dataset involving M. leprae-infected primary Schwann cells in vitro was performed. Through this comparison it was found that the expression profile obtained after a higher multiplicity of infection (100:1) for 48 h has a distinct activation pattern, with genes associated with lipid biosynthesis and consistent with a scavenging and phagocytic profile. Understanding the interaction between macrophage and M. leprae is highly relevant in deciphering key features that could further aid in the definition of the clinical forms.

Materials and Methods

Study Subjects

Volunteers were recruited from the staff of the Lauro de Souza Lima Institute (Secretary of Health of São Paulo State, Brazil) including health care, cleaning, and security workers, as well as students. Firstly, 29 healthy individuals (14 men and 15 women, from 20 to 30 years), free from cancer and infectious and autoimmune diseases, were enrolled for the collection of peripheral blood. A signed, written informed consent was obtained from all participants. The study was approved by the Local Ethics Committee (Protocol: 56169616.5.0000.5475).

Monocyte-Derived Macrophage (MDM) Differentiation

Fifty milliliters of peripheral blood were collected by venipuncture in tubes with anticoagulant (heparin). Peripheral blood mononuclear cells (PBMCs) were purified by density gradient (Histopaque 1077, Sigma Co., St. Louis, MO, USA) and monocytes were isolated by positive selection employing anti-CD14-coated magnetic microbeads (Miltenyi Biotec, Auburn, CA, USA) according to the manufacturer’s instructions. Purity greater than 95% was confirmed by flow cytometry using an anti-CD14 antibody (BD Biosciences, Franklin Lakes, NJ, USA).

Macrophages were differentiated by adherence, without using recombinant cytokines, in order to obtain a non-polarized macrophage with M0 phenotype, according described by Vogel and collaborators (15) with modifications. Briefly, monocytes were cultured in 24-well plates (5 × 10⁵ cells/well) in Iscove’s modified Dulbecco’s medium (IMDM, Gibco, Grand Island, NY, USA) supplemented with 10% of pooled human sera from the participants and 100 units/mL of penicillin and 100 μg/mL of streptomycin (Gibco; Thermo Fisher Scientific, Inc.) for 6 days at 37°C in 5% CO₂ humidified atmosphere for the differentiation into macrophages. Cultures were fed by replacing half part of the complete medium on the third day. The percentage of differentiated macrophages was evaluated during the standardization of the protocols considering the expression of CD68 and was greater than 95%.

M. leprae Purification and Macrophages Infection

M. leprae (Thai-53 strain, kindly provided by Dr. Yuji Yamamoto, NIH, Japan) was obtained from footpads of athymic nude mice, according to a previously described protocol (16). Briefly, each footpad was inoculated with 3 × 10⁶ acid-fast bacilli in 30 μL of saline solution. After 6-7 months, mice were euthanized and the footpads were collected (CEUA 06/2006). The protocol for bacilli recovery included footpad dissection, tissue isolation, and enzymatic digestion with 0.05% trypsin, followed by purification, quantification, and evaluation of ML viability using the Live/Dead BacLight Bacterial Viability Kit (Molecular Probes, Inc., Eugene, OR, USA). The M. leprae count was done after Ziehl-Neelsen staining (16, 17).

Previously differentiated MDMs were infected with M. leprae in the multiplicity of infection of 10 bacilli/cell (MOI = 10:1) and kept at 37°C in 5% CO₂ humidified atmosphere. After 48 h of infection, macrophages were collected in TRIzol Reagent (Thermo Fisher Scientific, Carlsbad, CA, USA). Uninfected MDMs incubated under the same conditions were used as controls.

RNA Isolation and Microarray

Total RNA was extracted by an in-house method using phenol:chloroform and isopropyl alcohol (18). Glycogen (Thermo Fisher Scientific) was added to improve RNA recovery. After centrifugation at 12,000 RPM for 5 minutes, the RNA pellet was washed with 70% ethanol, air dried, and resuspended in DEPC-treated water. RNA quality was evaluated on a Bioanalyzer 2100 Instrument (Agilent Technologies Inc., Palo Alto, CA, USA), by using the RNA 6000 Nano Kit (Agilent Technologies). For all samples, the RNA integrity number (RIN) was higher than 8.

For the transcriptome, RNA samples were then purified with the RNeasy MinElute Cleanup Kit (QIAGEN, Hilden, Germany). The reverse transcription synthesis and biotin labeling were performed with Epicentre TargetAmp Kit (Illumina, CA, USA), and transcriptomes were obtained by chip hybridization using the HumanHT-12 v4 BeadChip followed by scanning in iScan equipment (Illumina) according to manufacturer’s instructions.

Microarray Data Analysis

Raw.idat files were imported into the R v. 3.6.1 (BiocVersion 3.9.0) environment using limma v. 3.40.0. Quality control was carried out by investigating raw and normalized expression intensities across arrays with Tukey box plots. Background correction and between-array quantile normalization were performed using negative/positive control probes (19) from the manufacturer with limma::neqc function (Figure S1C) (20–22). Next, an ExpressionSet (Biobase v. 2.44.0) object was assembled to hold assay, gene, and phenotype data (23). Principal Component Analysis (PCA) was used to inspect dataset structure and biological and technical effects. PCA was conducted with FactoMineR v. 1.41 (24) and visualized with factoextra v. 1.0.5 and cowplot v. 1.0.0 (Figures S2A–F). Outlier samples were removed before attempting statistical inference based on the first three principal components (Figure S2). Differential expression analysis was performed by fitting gene-wise linear models with moderated standard errors by the empirical Bayes method (22, 25). The final model included the independent variables: chip (categorical with seven levels), sex (categorical, two levels), treatment (categorical, two levels), and individuals as a random effect. Genes were mapped to Entrezid and HGNC official symbols using the illuminaHumanv4.db v. 1.26.0 annotation. Duplicated Entrezid were removed by keeping the one with the largest average across all arrays. Finally, nominal P-values were inspected with histograms and adjusted for multiple testing with the Benjamini-Hochberg method to control the false-discovery rate (FDR) (26). Genes were considered differentially expressed (DE) if FDR ≤ 10% and absolute fold-change ≥ 1.5 (|log₂FC| ≈ 0.58) with alternative thresholds indicated wherever used. A volcano plot was drawn to illustrate the DE results with ggplot2 v. 3.3.0 (27). Exploratory hierarchical clustering was constructed with the top 92 DEGs using pheatmap v. 1.0.12 (28), with Euclidean distance (samples) plus average agglomeration and Pearson correlation (genes). Raw and normalized data are available in the Gene Expression Omnibus (GEO) accession GSE162416. Also, data and computer source code are readily available in Zenodo (https://dx.doi.org/10.5281/zenodo.4401968).

Over-Representation Analysis (ORA) and Gene Set Enrichment Analysis (GSEA)

Gene ontology (GO) biological processes (BP) were evaluated separately for up- and downregulated DE lists using clusterProfiler v. 3.12.0 (29), fgsea v. 1.10.0 (30) and org.Hs.eg.db v.3.8.2 annotations. The universe set contained all Entrezid genes (n = 21207) used in DE analyses. The minimum gene set size was 5 and the FDR cutoff was set at 10%. Enriched BP were visualized with dot plots or heat plots. For GSEA (Figure S3), the gene list was constructed using limma’s estimated log₂FC. GO BP and Reactome GSEA were estimated with 5000 permutations and with gene sets containing at least five genes (29, 31, 32).

Comparison to Infected Schwann Cell Dataset

Microarray dataset GSE35423 was processed as described elsewhere (12, 14). DEGs from Schwann cells infected with M. leprae for 48 hours were filtered by FDR ≤ 10% and |log₂FC| ≥ 0.26 (20% difference). Genes differentially expressed common to both datasets were visually compared with a dot plot along with log₂FC and confidence intervals from original results (ggplot2 v. 3.3.0). Since the number of common DEG was large, only the top 50 DEG (sorted by decreasing log₂FC) with same and opposite modulation signs were drawn. Table S6 contains the full results. UpSetR v.1.4.0 was used to visualize the intersection between DEG according to the dataset and modulation sign.

Gene Set Variation Analysis (GSVA)

Pathway activity was estimated using GSVA with custom gene sets (33). GSVA is an unsupervised non-parametric alternative to ORA and GSEA as it does not depend on prior selection or ranking of the genes from group comparisons. The score produced by GSVA can be interpreted as the coordinated activation of genes from the gene set, summarizing the expression profile within individual samples. The gene sets used herein were compiled from literature and functional annotation databases. The ‘granulomatosis’ gene set was assembled with genes sourced from the Human Phenotype Ontology (HP:0002955), Gene Ontology (GO:0002432), and DisGeneNet (granulomatous diseases with score ≥ 0.1) totaling 30 genes (Table S7). The macrophage polarization signatures are mainly from (34), where ‘M1’ (n genes = 25), ‘M2’ (n=20), ‘M2a’ (n=12), ‘M2b’ (n=9) and ‘M2c’ (n=12). ‘Pro-M1’ and ‘Pro-M2’ gene sets were built from literature (35, 36). Autophagy genes were retrieved from Reactome v.75 accessions: ‘Macroautophagy’ (R-HSA-1632852.8, n=137), ‘Chaperone Mediated Autophagy’ (R-HSA-9613829.3, n=22), and ‘Late endosomal microautophagy’ (R-HSA-9615710.3, n=34). Wilcoxon signed-rank test was used to test the differences between mean ranks between infected and mock macrophages signatures. Spearman’s rank correlation coefficient was calculated alongside 95% confidence intervals using DescTools R package v.0.99.40 (37). Principal component analysis (PCA) was applied to further explore the macroautophagy genes in the dataset. Mean-centered and variance standardized expression matrix with the 137 genes from macroautophagy gene set was subjected to PCA computed with FactoMineR v.2.4 (24). PCA scatter plot and contribution bar plots were graphed with factoextra v.1.0.6.

Results

Differentially Expressed Genes in MDMs Infected With M. leprae

Here, we included 29 healthy volunteers (15 females and 14 males) with an age range of 20-30 years. PBMCs were collected and used to obtain monocytes for differentiation to MDMs, which were subsequently infected with M. leprae at 10:1 MOI for 48 h. Microarray quality control was used to discard aberrant arrays, and outlier samples were removed (Figure S2).

We designed our study to evaluate early changes in the macrophage-M. leprae interaction. Thus, human MDMs were analyzed after 48 h of infection with M. leprae, which resulted in 325 unique upregulated (FDR ≤ 10% and log₂FC ≥ 0.58) and 117 downregulated (FDR ≤ 10% and log₂FC ≤ -0.58) genes. Genes with the largest effect size (|log₂FC| > 1.5) are annotated in the volcano plot of Figure 1A, of which some members of the metallothionein family can be observed, like MT1G, MT1E and MT1P. Figure 1B shows a heatmap with hierarchical clustering of all samples and the top 92 genes with an FDR ≤ 10% and |log₂FC| ≥ 1. The heatmap demonstrates a cluster pattern that clearly distinguishes infected from non-infected samples (Figure 1B), at the same time highlighting the high heterogeneity among individuals for some genes, such as ORM1, MT1G, MT1H, MMP12, CXCL5, COL22A1, and GAL. Some of the gene families identified such as metallothioneins, chemokines, and interleukins are involved in inflammation and autophagy, and among the induced genes, we found that MT1G, MMP7, TNFAIP6, CYP27B1, and CCL5/RANTES had the largest fold-changes (Figure 1A and Table S1). Conversely, among the most repressed were SELENOP, TLR7, CDCP1, HPSE, and GNG2 (Figure 1A and Table S1).

FIGURE 1

Figure 1 (A) Volcano plot showing the DEGs from monocyte-derived macrophages infected with live M. leprae (MOI 10:1) for 48 h. Blue dots represent genes with an FDR ≤ 10% and |log₂FC| ≥ 0.58. Gene symbols are given for those with an FDR ≤ 10% and |log₂FC| ≥ 1.5. (B) Heatmap and unsupervised hierarchical clustering of genes with an FDR ≤ 10% and |log₂FC| ≥ 1 (n = 92). Samples were clustered based on Euclidean distance and genes with Pearson correlation coefficient, both with average agglomeration. Color key displays expression values in standard deviation units away from the mean (i.e., scaled and centered row-wise). FDR, false discovery rate; FC, fold change.

To further understand which biological processes the DEGs were involved in, we performed an over-representation analysis (ORA) with Gene Ontology Biological Process annotation. ORA revealed several biological processes modulated by M. leprae infection, whereby the expression patterns of the MDM genes were consistent with resistant responses and a reprogramming profile that would induce M. leprae killing and infection control. Strikingly, the upregulated pathways had a greater number of DEGs and more robust FDR. Among these predominantly induced pathways, there were redundancies, which are expected since several of these genes participate in multiple pathways. Nevertheless, there is a clear activation of ‘regulation of inflammatory response’, ‘regulation of leukocyte activation’, ‘cytokine secretion’, ‘regulation of T cell activation’, ‘response to IFN-gamma’, as well as some others (Figure 2A and Table S2). On the contrary, downregulated genes pertain to biological processes such as membrane lipid metabolic process, sphingolipid metabolic process, and locomotory behavior (Figure 2B and Table S3).

FIGURE 2

Figure 2 Dot plot showing the top significant GO biological processes enriched from ORA of genes (A) upregulated (n = 35) or (B) repressed (n = 30) by M. leprae infection. Gene ratio is the fraction of genes belonging to an ontology over the total number of modulated genes. Circle size shows the number of modulated genes per biological process. FDR, false discovery rate.

Commonly Modulated Genes in Macrophages and Schwann Cells Upon Infection With Live-M. leprae

As a comparison, we decided to use a public dataset from M. leprae infected (100:1) Schwann cells where approximately 30-100 intracellular bacilli per cell were observed (2). In this case, a shift toward a de-differentiation phenotype was noticed. By comparing the lists of DEGs from this study and the results from Schwann cells infected with live M. leprae (100:1) for 48 h, we identified 920 DEGs with FDR ≤ 10% and a difference in mean expression of at least 20% (i.e., |log₂FC| ≥ 0.26). Of these, 229 (24.89%) were upregulated in both datasets, while 179 (19.45%) were jointly downregulated (Figure 3A and Table S6). On the other hand, 512 genes (55.65%) were regulated in opposite directions in the two experiments (Figure 3A). We next graphed the top 50 genes with the greatest |log₂FC| for both concordant and discordant DEG, as this can be used to further enlighten both similarities and differences in how the pathogen interacts with each host cell as well as the effect of a higher/lower multiplicity of infection (Table S6). Indeed, a higher multiplicity of infection (Schwann dataset) leads to activation of the type I IFN pathways (OAS1, TRIM7, TNFSF10) and subversion of energetic metabolism where increased glucose uptake is redirected from glycolysis/mitochondrial respiration to lipid metabolism (MVK, DHCR7, HMGCS1, LDLR, MSMO1)  (3, 12). In the same manner, some members of TNF signaling, inflammasome pathway, and regulators of NFκB and JNK pathways were upregulated in MDMs but repressed in Schwann cells, such as TNF, CXCL5, FAS, MAPK13, TRAF1, NKFB2, AIM2, and PANX1 (Figure 3C and Table S6). Furthermore, in Schwann cells, we observed a consistent upregulation of genes involved with prostaglandin biosynthesis, lipogenesis, mitochondrial metabolism, and negative regulation of immune ADORA1, LDLR, NOV, PPARA, SERPINF1, FFAR4, and ARG2). The patterns were quite distinct from the MDMs where Schwann cells showed a more pronounced expression of genes SELENOP, GAL, RGS18, AIF1, DEF6, ANGPTL6, DHCR7, MVK, and MSMO, some involved in cholesterol biosynthesis, along with LACC1/FAMIN, which is also genetically associated with leprosy (Table S6).

FIGURE 3

Figure 3 Common DEGs between this study and Schwann cells infected with M. leprae (MOI 1:100 for 48 h) [GSE35423] (12). (A) The number of commonly DEG by dataset and modulation sign passing |log₂FC| ≥ 0.26 and FDR ≤ 10% in both datasets. Arrows indicate (↑) upregulation and (↓) downregulation. Top 50 DEG common to both datasets modulated in the same manner (B) or oppositely (C). Points represent log₂FC (unstandardized) from each dataset original differential expression analysis alongside error lines indicating nominal 95% confidence intervals. All genes shown have an FDR ≤ 10% and a difference in mean expression of at least 56% (≈ 0.65 log₂FC).

Infected MDMs Express Mainly M1 Polarization Genes That Are Correlated With Granuloma Formation

The somewhat contrasting gene expression profiles identified before suggested that the infected MDMs could be expressing M1 markers with autophagic and phagocytic profiles. To examine this, we calculated gene scores representative of multiple macrophage polarization phenotypes, autophagy, and granulomatosis, which are features often observed in paucibacillary leprosy patients. Figure 4A shows the scores for each signature indicating a transcriptional activity increase in genes involved with M1 macrophage polarization, as well as decreased or similar activity for M2, M2a, and Pro-M2. It seems that although genes responsible for inducing M2 phenotype are unaffected, the M2b signature shows activation upon infection, which could indicate either a mixed MDM specialization phenotype or distinctive polarization patterns between individuals, like the clinical disease presentation. As for the autophagy signatures, macroautophagy is specifically active in infected MDMs, whereas chaperone-mediated autophagy appears the opposite (Figure 4A). Finally, we tested the correlation between granulomatosis gene scores with M1, M2, and macroautophagy scores. There was a moderated positive correlation between M1 polarization signature with granulomatosis, and no correlation with M2 score, which together corroborated our hypothesis that a pro-inflammatory macrophage profile seems predominant in this model with low MOI (Figures 4B, C). The granulomatosis signature does not correlate with macroautophagy, which could indicate that these processes may not be simultaneously active at the transcriptional level (Figure 4D). Finally, we wanted to see which macroautophagy genes were driving the signature activity signal. Principal component analysis (PCA) with the 137 genes from the macroautophagy signature revealed that the second principal component (PC2) separated MDMs according to the infection treatment, explaining around 12% of the variability (Figure 4E). As expected, there is considerable between-individual variation and noise captured mainly by the first and second principal components (PCs). We then explored PC1, and only the “microarray chip” variable partially correlated with this PC (data not shown). Some uninfected MDMs had basal higher expression (placed vertically higher across PC2) as seen in blue dots over the horizontal dashed line marking zero (Figure 4E). Next, we extracted the top 25 genes most correlating to PC2. The gene TOMM40 alone contributed with explaining almost 4% of the variability of that axis (Figure 4F), where more than 25 genes contribute with more variability than would be expected if contributions were uniform (dashed vertical line in Figure 4F). Infected MDMs separated from their mock pair by great distance vertically (Figure 4F) are cells with the most difference in expression for the genes correlating with PC2, which varies among individuals. Most of the top genes with higher contributions to PC2 showed upregulation after infection with live M. leprae, with some exceptions that were significantly downregulated (Figure 4G).

FIGURE 4

Figure 4 (A) Gene scores for macrophage polarization profiles, granulomatosis and autophagy. Tukey box plots display first, second (median), and third quartiles with whiskers extending ± 1.5 × interquartile range (IQR). Each point illustrates a MDM culture (n=50) from one human donor in mock and infected conditions (paired within donor). Nominal P-values displayed are from Wilcoxon signed-rank test. (B–D) Scatter plots illustrating the correlation between granulomatosis, M1, M2, and macroautophagy gene scores. Spearman’s rank correlation coefficient is shown with nominal 95% confidence intervals calculated only with the infected condition (n=25). Slope and intercept for drawing the lines were estimated using robust linear regression from the MASS v.7.3-53.1 R package. (E) Scatter plot with the two first principal components from PCA on the subset of 137 macroautophagy-related genes. Lines connect samples from the same donor. (F) Percentage contribution of the top 50 genes most correlated to PC2. (G) Log₂FC estimates from differential expression analysis for the top 50 macroautophagy-related genes from PCA and their nominal 95% confidence intervals. Light-shaded red indicates gene with FDR > 10%.

Discussion

Here, we evaluated the gene expression patterns of M. leprae-infected MDMs, at a 10:1 MOI for 48 h, using a large-scale technique. The data were consistent with an M1-like differentiation profile, with several enriched pathways associated with increased microbicidal activity, T cell activation, and IFN-gamma cytokine secretion. Among the most DE genes were those found in the vitamin D processing pathway (CYP27B1, VDR, IL1B), CC-chemokine ligands (CCL5, CCL4, CCL3, CCL4L2, CCL2), receptors involved in cellular migration (GREM1, MCOLN2, MDK, CCL3), and inflammation (CCL5/RANTES). It seems that a combination of MDM differentiation with M. leprae infection at a low multiplicity of infection (10:1) induces a gene expression program in agreement with a protective response. One of these genes, CCL5/RANTES, is a key chemoattractant for monocytes in the skin suggesting a pro-inflammatory feedback loop during M. leprae infection toward M1-like macrophages, which is consistent with the higher expression in PB leprosy (38, 39). Curiously, the chemokine-clustered genomic region has been associated with leprosy indicating that genetic variations within this gene could contribute to clinical form polarization (40). We observed an upregulation of metallothionein gene’s expression, like MT1E and MT1G. These genes produce proteins that regulate metal availability, while zinc homeostasis is crucial to activation of transcriptional factors and reactive oxygen species and have been associated with antimicrobial immune responses (41).

Some metalloproteinases, such as MMP7, ADAMDEC1 and MMP12, were also highly differentially expressed. These are involved with tissue/matrix remodeling and host defense. It is interesting that the macrophages in this model appear to have a regulatory activation of gene expression involved in tissue remodeling. Indeed, MMP2 and MMP9 expression and activity were higher among T-pole patients (42). This pattern is dependable with the reshaping of cellular morphology and the adjacent tissue during M1-like macrophage differentiation observed here.

The upregulation of mitochondrial genes in M. leprae-infected-macrophages, as compared to Schwann cells, is also consistent with the aforementioned alterations observed in infected MDMs toward epithelioid transformation programming. In this same direction, NDUFAF6, which encodes a protein that participates in the mitochondrial respiratory chain complex I (NADH:ubiquinone oxidoreductase) assembly; TACO1, which is involved in the translational activation of mitochondria cytochrome c; and TOMM40, the encoded protein of which produces a pore to channel protein precursors into mitochondria were also observed. We also detected an over-expression in MDMs of metabolite transporters across the inner mitochondrial membrane (SLC25A12).

Other sets of induced genes in the MDMs were involved in tubulin assembly and cytoskeleton remodeling, such as TUBB6. Some of these genes are involved in autophagy and the mTOR pathway (43), which is increased during macrophage differentiation in the presence of a low M. leprae MOI or dead mycobacteria (44). Curiously, Yang and colleagues showed that a 20:1 MOI leads to higher CD163 and an anti-inflammatory profile with higher IL-10 and lower HLA expression, although dead M. leprae induced proinflammatory cytokines (45). In clinical samples these patterns are observed in low bacterial index patients (T-pole) and type 1 reactional patients (RR) where vitamin D-mediated microbicidal clearance is observed. Furthermore, IFN-gamma and IL-15 have been demonstrated to gear the polarization of M1-like macrophages, exhibiting a phenotype that is probably orchestrating the milieu inducing granuloma formation and autophagy (44, 46). Thus, there is a clear antagonistic pattern that, on the one hand, we observed granuloma formation and M1-like polarization for macrophages (10:1 MOI), whereas lipid biogenesis, wound healing and a phagocytic signature were observed for Schwann cells (100:1 MOI). This M2-like profile occurred in both Schwann cells and macrophages when higher concentrations of M. leprae were used to infect the cells (11, 12, 47). Indeed, molecular and biochemical subversion induced by M. leprae infection leads to a type I IFN response reducing autophagy (12, 48) and turning on the Warburg-like effect (3). This is expected since this downregulation was detected as associated with diminished capacity to produce ATP through the respiratory chain (3) and the redirection of glycolysis to lipid biosynthesis. Furthermore, we observed a prominent upregulation toward higher MOI (100:1) when genes associated with cholesterol and fatty acid metabolism (DHCR7, MVK, and MSMO, and LACC1/FAMIN were analyzed, which are pathways involved in lepromatous leprosy immunopathogenesis.

Genes associated with leprosy outcome are also involved with ulcerative colitis pathogenesis and other diseases, such as Crohn’s disease, Parkinson’s and Alzheimer’s (49). Curiously, some of the genes highlighted here, such as TOMM40 and TUBB6, were associated not only with neurodegenerative diseases, like Alzheimer’s, but also with ulcerative colitis (50, 51). These proteins interact with LRRK2, which is a gene that has been independently associated with leprosy, and Parkinson’s disease (49, 52). Another important gene upregulated in M. leprae-Schwann cells was LACC1. The protein which LACC1/FAMIN encodes is involved in fatty acid oxidation and bacterial clearance. Interestingly, the gene was associated with leprosy in several studies (53, 54).

The global gene expression can also be analyzed in a genotype-phenotype correlation perspective. In leprosy, whole blood cells stimulated with M. leprae sonicates indicate quantitative trait loci (eQTL) associated with transcript levels when samples were compared before and after stimulation. The data revealed SNPs controlling immunoinflammatory responses such as type II IFNs and bacterial/pathogen recognition (55). Indeed, there are SNPs regulating expression of genes such as LACC1, which indicates that it is likely that a combination of M. leprae-induced activation with genetic polymorphisms may define the commitment toward the T- or L-pole. Thus, together these could contribute to identifying potential pharmacologic targets as adjuvants to personalize treatment for each specific clinical form.

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 below: https://www.ncbi.nlm.nih.gov/geo/, GSE162416.

Ethics Statement

The studies involving human participants were reviewed and approved by Local Ethics Committee, Lauro de Souza Lima Institute. The patients/participants provided their written informed consent to participate in this study. The animal study was reviewed and approved by Local Ethics Committee, Lauro de Souza Lima Institute.

Author Contributions

AP, MM, and VB contributed to the design and implementation of the research. TL-C, AL, and MM to the writing of the manuscript. BM, DB, RC, and GG performed the experiments and processed the experimental data. TL-C performed the analysis. PR and DB contributed to M. leprae maintenance in athymic nude mice and purification. TL-C, MM, and AP to the analysis of the results. All authors provided critical feedback and helped shape the research, analysis and manuscript. All authors contributed to the article and approved the submitted version.

Funding

This study was funded in part by the Fundação de Amparo do Estado do Rio de Janeiro (FAPERJ, E_09/2019 - Cientista do Nosso Estado – 2019 and E_34/2014 - PENSA RIO – Apoio ao Estudo de Temas Relevantes e Estratégicos para o RJ - 2014) and the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) [grant #313657/2018-1 and 4000170/2017-2]. Fundação de Amparo à Pesquisa do Estado de São Paulo - FAPESP, Grant no. 2015/01744-9.

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.

Acknowledgments

TL-C was supported by a scholarship from the National Council for Scientific and Technological Development (Conselho Nacional de Desenvolvimento Científico e Tecnológico, CNPq) from 2018 to 2022.

Supplementary Material

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

References

1. Scollard DM, Adams LB, Gillis TP, Krahenbuhl JL, Truman W, Williams DL. The Continuing Challenges of Leprosy The Continuing Challenges of Leprosy. Clin Microbiol Rev (2006) 19:338–81. doi: 10.1128/CMR.19.2.338

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Masaki T, Qu J, Cholewa-Waclaw J, Burr K, Raaum R, Rambukkana A. Reprogramming adult Schwann cells to stem cell-like cells by leprosy bacilli promotes dissemination of infection. Cell (2013) 152:51–67. doi: 10.1016/j.cell.2012.12.014

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Medeiros RCA, Girardi K do C de V, Cardoso FKL, Mietto B de S, Pinto TG de T, Gomez LS, et al. Subversion of Schwann Cell Glucose Metabolism by Mycobacterium leprae. J Biol Chem (2016) 291:21375–87. doi: 10.1074/jbc.M116.725283

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Cole ST, Eiglmeier K, Parkhill J, James KD, Thomson NR, Wheeler PR, et al. Massive gene decay in the leprosy bacillus. Nature (2001) 409:1007–11. doi: 10.1038/35059006

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Delavary BM, van der Veer WM, van Egmond M, Niessen FB, Beelen RHJ. Macrophages in skin injury and repair. Immunobiology (2011) 216:753–62. doi: 10.1016/j.imbio.2011.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Mantovani A, Biswas SK, Galdiero MR, Sica A, Locati M. Macrophage plasticity and polarization in tissue repair and remodelling. J Pathol (2013) 229:176–85. doi: 10.1002/path.4133

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Kashem SW, Haniffa M, Kaplan DH. Antigen-Presenting Cells in the Skin. Annu Rev Immunol (2017) 35:469–99. doi: 10.1146/annurev-immunol-051116-052215

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Cardoso CC, Pereira AC, de Sales Marques C, Moraes MO. Leprosy susceptibility: genetic variations regulate innate and adaptive immunity, and disease outcome. Future Microbiol (2011) 6:533–49. doi: 10.2217/fmb.11.39

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Montoya D, Cruz D, Teles RMB, Lee DJ, Ochoa MT, Krutzik SR, et al. Divergence of Macrophage Phagocytic and Antimicrobial Programs in Leprosy. Cell Host Microbe (2009) 6:343–53. doi: 10.1016/j.chom.2009.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Inkeles MS, Teles RMB, Pouldar D, Andrade PR, Madigan CA, Lopez D, et al. Cell-type deconvolution with immune pathways identifies gene networks of host defense and immunopathology in leprosy. JCI Insight (2016) 1:e88843. doi: 10.1172/jci.insight.88843

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Teles RMB, Graeber TG, Krutzik SR, Montoya D, Schenk M, Lee DJ, et al. Type I Interferon Suppresses Type II Interferon-Triggered Human Anti-Mycobacterial Responses. Science (2013) 339:1448–53. doi: 10.1126/science.1233665

PubMed Abstract | CrossRef Full Text | Google Scholar

12. de Toledo-Pinto TG, Ferreira ABR, Ribeiro-Alves M, Rodrigues LS, Batista-Silva LR, Silva BJ de A, et al. STING-Dependent 2′-5′ Oligoadenylate Synthetase–Like Production Is Required for Intracellular Mycobacterium leprae Survival. J Infect Dis (2016) 214:311–20. doi: 10.1093/infdis/jiw144

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Ma Y, Pei Q, Zhang L, Lu J, Shui T, Chen J, et al. Live Mycobacterium leprae inhibits autophagy and apoptosis of infected macrophages and prevents engulfment of host cell by phagocytes. Am J Transl Res (2018) 10:2929–39.

PubMed Abstract | Google Scholar

14. Leal-Calvo T, Moraes MO. Reanalysis and integration of public microarray datasets reveals novel host genes modulated in leprosy. Mol Genet Genomics (2020) 295:1355–68. doi: 10.1007/s00438-020-01705-6

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Vogel DYS, Glim JE, Stavenuiter AWD, Breur M, Heijnen P, Amor S, et al. Human macrophage polarization in vitro: maturation and activation methods compared. Immunobiology (2014) 219:695–703. doi: 10.1016/j.imbio.2014.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Trombone APF, Pedrini SCB, Diório SM, Belone A de FF, Fachin LRV, do Nascimento DC, et al. Optimized protocols for Mycobacterium leprae strain management: frozen stock preservation and maintenance in athymic nude mice. J Vis Exp (2014) (85):50620. doi: 10.3791/50620

CrossRef Full Text | Google Scholar

17. Shepard CC. Acid-fast bacilli in nasal excretions in leprosy, and results of inoculation of mice. Am J Hyg (1960) 71:147–57. doi: 10.1093/oxfordjournals.aje.a120098

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Chomczynski P. A reagent for the single-step simultaneous isolation of RNA, DNA and proteins from cell and tissue samples. Biotechniques (1993) 15:532–534, 536–537.

PubMed Abstract | Google Scholar

19. Shi W, Oshlack A, Smyth GK. Optimizing the noise versus bias trade-off for Illumina whole genome expression BeadChips. Nucleic Acids Res (2010) 38:e204–4. doi: 10.1093/nar/gkq871

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Smyth GK, Speed T. Normalization of cDNA microarray data. Methods (2003) 31:265–73. doi: 10.1016/S1046-2023(03)00155-5

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Phipson B, Lee S, Majewski IJ, Alexander WS, Smyth GK. Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression. Ann Appl Stat (2016) 10:946–63. doi: 10.1214/16-AOAS920

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Huber W, Carey VJ, Gentleman R, Anders S, Carlson M, Carvalho BS, et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat Methods (2015) 12:115–21. doi: 10.1038/nmeth.3252

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Lê S, Josse J, Husson F. FactoMineR : An R Package for Multivariate Analysis. J Stat Softw (2008) 25:1–18. doi: 10.18637/jss.v025.i01

CrossRef Full Text | Google Scholar

25. Smyth GK. Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments. Stat Appl Genet Mol Biol (2004) 3:1–26. doi: 10.2202/1544-6115.1027

CrossRef Full Text | Google Scholar

26. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing. J R Stat Soc Series B (1995) 57(1):289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

27. Wickham H. ggplot2: Elegant Graphics for Data Analysis. New York: Springer-Verlag (2016).

Google Scholar

28. Kolde R. pheatmap: Pretty Heatmaps. R package version 1.0.12 (2019). Available at: https://cran.r-project.org/package=pheatmap.

Google Scholar

29. Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters. OMICS: A J Integr Biol (2012) 16:284–7. doi: 10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

30. Korotkevich G, Sukhov V, Sergushichev A. Fast gene set enrichment analysis. bioRxiv [Preprint] (2016). Available at: https://www.biorxiv.org/content/10.1101/060012v3 (Accessed April 10, 2021).

Google Scholar

31. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U States America (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

CrossRef Full Text | Google Scholar

32. Mootha VK, Lindgren CM, Eriksson K-F, Subramanian A, Sihag S, Lehar J, et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet (2003) 34:267–73. doi: 10.1038/ng1180

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

34. Benoit M, Desnues B, Mege J-L. Macrophage polarization in bacterial infections. J Immunol (2008) 181:3733–9. doi: 10.4049/jimmunol.181.6.3733

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Murray PJ. Macrophage Polarization. Annu Rev Physiol (2017) 79:541–66. doi: 10.1146/annurev-physiol-022516-034339

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Ni Y, Zhuge F, Nagashimada M, Ota T. Novel Action of Carotenoids on Non-Alcoholic Fatty Liver Disease: Macrophage Polarization and Liver Homeostasis. Nutrients (2016) 8:391. doi: 10.3390/nu8070391

CrossRef Full Text | Google Scholar

37. Signorell A. DescTools: Tools for Descriptive Statistics. R package version 0.99.40. (2021). Available at: https://cran.r-project.org/package=DescTools.

Google Scholar

38. Belone A de FF, Rosa PS, Trombone APF, Fachin LRV, Guidella CC, Ura S, et al. Genome-wide screening of mRNA expression in leprosy patients. Front Genet (2015) 6:334. doi: 10.3389/fgene.2015.00334

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zhang D-F, Wang D, Li Y-Y, Yao Y-G. Integrative analyses of leprosy susceptibility genes indicate a common autoimmune profile. J Dermatol Sci (2016) 82:18–27. doi: 10.1016/j.jdermsci.2016.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Jamieson SE, Miller EN, Black GF, Peacock CS, Cordell HJ, Howson JMM, et al. Evidence for a cluster of genes on chromosome 17q11-q21 controlling susceptibility to tuberculosis and leprosy in Brazilians. Genes Immun (2004) 5:46–57. doi: 10.1038/sj.gene.6364029

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Subramanian Vignesh K, Deepe. GS. Metallothioneins: Emerging Modulators in Immunity and Infection. Int J Mol Sci (2017) 18:2197. doi: 10.3390/ijms18102197

CrossRef Full Text | Google Scholar

42. Teles RMB, Teles RB, Amadeu TP, Moura DF, Mendonca-Lima L, Ferreira H, et al. High Matrix Metalloproteinase Production Correlates with Immune Activation and Leukocyte Migration in Leprosy Reactional Lesions. Infect Immun (2010) 78:1012–21. doi: 10.1128/IAI.00896-09

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Jha AK, Huang SC-C, Sergushichev A, Lampropoulou V, Ivanova Y, Loginicheva E, et al. Network Integration of Parallel Metabolic and Transcriptional Data Reveals Metabolic Modules that Regulate Macrophage Polarization. Immunity (2015) 42:419–30. doi: 10.1016/j.immuni.2015.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Silva BJ de A, Barbosa MG de M, Andrade PR, Ferreira H, Nery JA da C, Côrte-Real S, et al. Autophagy Is an Innate Mechanism Associated with Leprosy Polarization. PloS Pathog (2017) 13:e1006103. doi: 10.1371/journal.ppat.1006103

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Yang D, Shui T, Miranda JW, Gilson DJ, Song Z, Chen J, et al. Mycobacterium leprae-Infected Macrophages Preferentially Primed Regulatory T Cell Responses and Was Associated with Lepromatous Leprosy. PloS Neglect Trop Dis (2016) 10:e0004335. doi: 10.1371/journal.pntd.0004335

CrossRef Full Text | Google Scholar

46. Chun Y, Kim J. Autophagy: An Essential Degradation Program for Cellular Homeostasis and Life. Cells (2018) 7(12):278. doi: 10.3390/cells7120278

CrossRef Full Text | Google Scholar

47. Guerreiro LTA, Robottom-Ferreira AB, Ribeiro-Alves M, Toledo-Pinto TG, Rosa Brito T, Rosa PS, et al. Gene Expression Profiling Specifies Chemokine, Mitochondrial and Lipid Metabolism Signatures in Leprosy. PloS One (2013) 8(6):e64748. doi: 10.1371/journal.pone.0064748

PubMed Abstract | CrossRef Full Text | Google Scholar

48. de Toledo-Pinto TG, Batista-Silva LR, Medeiros RCA, Lara FA, Moraes MO. Type I Interferons, Autophagy and Host Metabolism in Leprosy. Front Immunol (2018) 9:806. doi: 10.3389/fimmu.2018.00806

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Fava VM, Dallmann-Sauer M, Schurr E. Genetics of leprosy: today and beyond. Hum Genet (2020) 139:835–46. doi: 10.1007/s00439-019-02087-5

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Law BMH, Spain VA, Leinster VHL, Chia R, Beilina A, Cho HJ, et al. A direct interaction between leucine-rich repeat kinase 2 and specific β-tubulin isoforms regulates tubulin acetylation. J Biol Chem (2014) 289:895–908. doi: 10.1074/jbc.M113.507913

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Guo Y, Xu W, Li J-Q, Ou Y-N, Shen X-N, Huang Y-Y, et al. Genome-wide association study of hippocampal atrophy rate in non-demented elders. Aging (Albany NY) (2019) 11:10468–84. doi: 10.18632/aging.102470

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Fava VM, Xu YZ, Lettre G, Van Thuc N, Orlova M, Thai VH, et al. Pleiotropic effects for Parkin and LRRK2 in leprosy type-1 reactions and Parkinson’s disease. Proc Natl Acad Sci (2019) 116:15616–24. doi: 10.1073/pnas.1901805116

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Cader MZ, Boroviak K, Zhang Q, Assadi G, Kempster SL, Sewell GW, et al. C13orf31 (FAMIN) is a central regulator of immunometabolic function. Nat Immunol (2016) 17:1046–56. doi: 10.1038/ni.3532

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Wang D, Fan Y, Malhi M, Bi R, Wu Y, Xu M, et al. Missense Variants in HIF1A and LACC1 Contribute to Leprosy Risk in Han Chinese. Am J Hum Genet (2018) 102:794–805. doi: 10.1016/j.ajhg.2018.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar